Improving read alignment through the generation of alternative reference via iterative strategy

There is generally one standard reference sequence for each species. When extensive variations exist in other breeds of the species, it can lead to ambiguous alignment and inaccurate variant calling and, in turn, compromise the accuracy of downstream analysis. Here, with the help of the FPGA hardware platform, we present a method that generates an alternative reference via an iterative strategy to improve the read alignment for breeds that are genetically distant to the reference breed. Compared to the published reference genomes, by using the alternative reference sequences we built, the mapping rates of Chinese indigenous pigs and chickens were improved by 0.61–1.68% and 0.09–0.45%, respectively. These sequences also enable researchers to recover highly variable regions that could be missed using public reference sequences. We also determined that the optimal number of iterations needed to generate alternative reference sequences were seven and five for pigs and chickens, respectively. Our results show that, for genetically distant breeds, generating an alternative reference sequence can facilitate read alignment and variant calling and improve the accuracy of downstream analyses.

Whole-genome sequencing provides a comprehensive method to identify genomic variations 1,2 . As next-generation sequencing (NGS) has generated an ever-increasing volume of genomic data, it is not uncommon for current studies to involve hundreds or even thousands of individuals [3][4][5] . The accumulation of sequencing data provides great insight into biological problems and also improves the rigorousness and comprehensiveness of genomic analyses, especially for population genetics and association studies 6,7 . However, massive sequencing data processing and genome variant calling impose a heavy computational and storage burden and have thus become an emergent issue for large-scale studies. To eliminate computing bottlenecks, efforts have been made to develop more efficient algorithms or to utilize parallel computing 8 . Recently, heterogeneous computing with FPGA (Field-Programmable Gate Array) accelerators has shown significant potential to produce significant improvements in the computing efficiency of short-read alignment and variant calling while maintaining a very high level of consistency between the output and the original method [9][10][11] . For example, Menges F et al. proposed a new base-calling algorithm that was implemented in FPGA to achieve real-time performance 12 . Using FPGA, Arram J et al. accelerated the alignment of short reads 28 times faster than Bowtie2 running with 16 threads on dual Intel Xeon E5-2640 CPUs 13 . Acceleration via FPGA enables us to perform time consuming analyses that could not be done previously.
In addition to increasing computing efficiency, it is even more important to ensure that effective and accurate information is extracted from the sequencing data. Despite substantial genetic variations being found across breeds, only one complete reference genome is generally available for each species. For example, the reference genome for pigs came from a domesticated duroc pig 14,15 , while the reference genome for chickens is based on a wild red junglefowl 16,17 , which is the ancestor of the domestic chicken 18,19 . In most cases, the closest reference genome will be used for the read alignment and subsequent variant calling 20,21 . Unfortunately, in some cases, the closest reference genome is not very similar-especially for domesticated plants and animals, among which strong artificial and natural selection has led to extensive genetic differences between the domesticates and their Scientific Reports | (2020) 10:18712 | https://doi.org/10.1038/s41598-020-74526-7 www.nature.com/scientificreports/ wild counterparts or among different domesticated breeds. For example, a recent study on Chinese indigenous pigs identified more new SNPs (Single nucleotide polymorphisms) than those recorded in the dbSNP database 22 .
A high-quality and representative genome assembly can greatly facilitate studies. Currently, genome analyses at the population level often include individuals from multiple breeds. A genetic distance between the reference genome and the individual under investigation that is too large can lead to ambiguous alignment and inaccurate variant calling and, in turn, compromise the accuracy of the analysis. Although the costs of genome sequencing continue to fall, it remains unrealistic to build a reference genome for every breed of species since De Novo genome assembling is still technically cumbersome and very expensive. The generation of alternative reference sequences, therefore, is a cost-effective approach to satisfy the needs of this type of research. Cho et al. built a Korean consensus reference by incorporating common variants in the Korean population and found that a consensus reference can be beneficial for efficient variant detection 23 . By merging the alignment results, Okumur K et al. constructed alternative consensus references for Mycobacterium tuberculosis. An empirical evaluation showed that the use of a consensus reference significantly improved mapping efficacy and facilitated phylogenetic analysis 24 . Although researchers have realized that the standard reference genome might not always perform well for specific studies, there is still a lack of systematic research on which cases require an alternative reference sequence and how an alternative reference sequence should be generated.
In this study, with the help of the FPGA hardware platform, we performed an extensive evaluation of the generation of alternative reference sequences. Instead of using a one round substitution approach, we employed an iterative strategy, through which highly variable regions were recovered as the number of iterations increased. We showed that this process improves read alignment and variant calling not only for a genetically distant target breed but also for other breeds that are distant to the reference breed but close to the target breed. By considering a balance between sensitivity and computing costs, we also evaluated the optimal sequencing coverage and iterations that were required to generate an alternative reference sequence. Our results provide the first comprehensive study on the effective generation of an alternative reference sequence.

Results
Accurate variant calling by GTX-One. Since NA12878 is the gold standard publicly available variant set for variant caller benchmarking, to evaluate the performance of variant calling using the GTX-One platform, 30× (or 90 Gb) of whole genome sequencing data of NA12878 (H1) was used for variant calling using both the GTX-One and the GATK Best Practice (GPB) workflow. Based on the Genome in a Bottle Consortium (GIAB) gold standard callset, we defined the true positives (variants called with the same genotype as the gold standard callset, TP), false positives (variants called but not in the gold standard callset, FP), and false negatives (variants in the gold standard callset but not called, FN). We calculated the precision and sensitivity for both the GTX-One and GBP workflow using the following formulas: precision = TP/(TP + FP) and sensitivity = TP/(TP + FN). For SNPs, the GTX-One achieves a high precision of 99.54% and a high sensitivity of 99.36%. Overall, the performance of GTX-One is nearly identical to that of the GBP workflow, but it is much more efficient (Table 1 and  Supplementary Table S1).
Differences in the mapping rates among domestic breeds. Domesticated plants and animals are subject to directional selection. Moreover, the demographic effects of isolation and genetic drift also change the allele frequencies of a population. Over time, these factors work together to promote genetic divergence across breeds. This is why the standard reference genome sequences do not always perform well. To reveal the differences in the mapping efficiency for distinct breeds, we chose genetically different Chinese domestic pigs and European commercial pigs, as well as Chinese domestic chickens and commercial layer chickens, for comparison. Since the quality of the reference genome can also affect mapping efficiency, we repeated the mapping process for both the latest and the previous versions of the genome assemblies.
For the pig species, we included five breeds and two reference genome versions: Sscrofa10.2 and Sscrofa11.1. Duroc is the breed from which the pig reference genome was built. Similarly, we included six breeds and two reference genome versions, galGal5 and galGal6, for chickens. Red junglefowl is the breed from which the chicken reference genome was built. As listed in Table 2 and Supplementary Table S2, even when controlled for the version of the reference genome, the mapping rates still varied among breeds of the same species. For pig species, the mapping rates of different breeds showed larger variations compared to those of chickens. As expected, for the commercial breeds DU and LD, the mapping rates were higher than those of Chinese indigenous breeds. For chicken species, the overall mapping rates were relatively higher and less variable among different breeds; however, approximate 2% differences were still observed between the best and worst scenarios. In addition, the results show that the mapping rate was significantly influenced by the quality of the reference assembly (see Table 1. A comparison of the FPGA and CPU implementation of variant calling of the NA12878 data.  Table S3 for statistics on the reference assemblies), indicating the necessity to use high-quality reference sequences for genomic analysis. We next downloaded the genome sequencing data for one human (H2) and one chimpanzee (C1) as a control to evaluate the degree of mapping differences. The genome sequencing data for the human and the chimpanzee were both purposely aligned to the GRCh37 and GRCh38 reference genome. Surprisingly, we found that the mapping rate of the chimpanzee data to the human genome was as high as 97.80%, which is only 2% less than human data, as listed in Table 3. The difference in the mapping rate between humans and chimpanzees is even smaller than the within-species differences we observed in the pig species, indicating that the genomes of domestic animals changed significantly during the process of domestication. Thus, if only the standard reference genome sequences were used, the analysis results might be compromised, especially for studies including multiple domesticated breeds.
Optimal sequencing coverage for accurate genotype calls. At higher levels of genome coverage, the called variants afford a higher degree of confidence because each base is covered by a greater number of aligned reads. However, a higher coverage of sequencing means higher costs, while sequencing coverage that is too low often causes inaccurate genotype calls. Based on the above results, high-quality reference genome assemblies, Sscrofa11.1 and galGal6, were used for further analyses of pigs and chickens, respectively. Using WZS as an example, the statistics in Table 4 show that with an increase in sequencing coverage, both the number of variants called, and the sensitivity increased. The mapping rate, however, remained consistently high for all sequencing coverage rates, suggesting good and stable mapping quality. Statistics for the other pig and chicken breeds are listed in Supplementary Table S4 and Supplementary Table S5.
To determine the optimal genome sequencing coverage for accurate genotype calls, we plotted the variation counts (log-transformed) against the coverage. The curve was found to be best fitted by a logistic function ( Fig. 1, Supplementary Fig. S1 and Supplementary Fig. S2). We thus determined the threshold to be 0.0001 by selecting the slope of the tangent of the curve, which corresponds to 18.2× (rounded up to 19×) in pigs (Table 5) and to 15.4× (rounded up to 16×) in chickens (Supplementary Table S6). We applied the same method to the sensitivity and coverage ratio, and the results were similar. Thus, we report the results using variation counts unless otherwise specified.
Improving mapping by using alternative references generated by an iterative strategy. In short, the aim of generating alternative references is to increase the mapping rate. For highly divergent regions containing consecutive mismatches, the sequencing reads are not directly mappable, so substitutions to the reference base cannot be easily made. To solve this problem, we employed an iterative strategy for which the consecutive mismatches are substituted step by step, as illustrated in Fig. 2.
After the optimal sequencing coverage was determined, the raw reads were sampled at the optimal coverage for variant calling. The callset was then used as the input for the GATK FastaAlternateReferenceMaker function to generate a new reference sequence. In the first round of iterations, the publicly released reference genome sequence was used. In the following iterations, the updated reference and the callset from the updated reference were used as inputs for the next iteration. Since we had no prior knowledge of how many iterations would be enough, we performed 30 iterations for both species. To test how an alternative reference sequence improves mapping, in each round, the updated WZS alternative reference sequence was recorded, and the sequencing reads from BMX, SZL, and LD pig breeds were mapped against it. The same was done for the updated YJ alternative reference sequence for chickens, and the sequencing reads from the LS, LDH, and BLK chicken breeds were mapped against it. As shown in Table 6 and Supplementary Table S7, increased mapping rates were observed for the alternative references. The maximum increases in the mapping rate for pigs were 1.51% for WZS, a 1.16% increase for BMX, and 1.21% for SZL. For chicken breeds, the maximum increases in the mapping rate were 0.18% for YJ, 0.42% for ZJ, and 0.37% for LS. Statistics of the variant counts, mapping rates, and coverage ratios for the other breeds are detailed in Supplementary Tables S7 and S8. To exclude the possibility that the increased mapping rate is caused by chimeric sequences, we checked if the alternative reference sequence produced actually would become closer to the target breed. We aligned the progressive WZS alternative reference sequences against the public WGZ genome using Minimap2. As shown in table S9, the alternative reference genome became  www.nature.com/scientificreports/  www.nature.com/scientificreports/ progressively more similar to the WZS genome. Overall, using alternative reference sequences improved mapping for genetically distant breeds, which would help researchers analyze population level data more effectively. We noticed that the variant counts of WZS and BMX in Table 6 moved up and down during the iterations. Thus, we manually checked the output VCF files during each iteration. For the examples shown in Supplementary Table S10, we found that inconsistencies were caused by the switching over of heterozygous alleles in the individuals used to generate the alternative reference sequences, while the same loci in unrelated individuals, for read alignment, were only homozygous. This suggests that substitution to the reference base would only be meaningful for homozygous loci.
Since the iteration process is time consuming and the drops in variation counts with additional iterations (WZS 0 in Table 6) become less significant, to determine the optimal iteration number, we fitted the variation counts against the number of iterations. Again, the best fitting model was a logistic function (Fig. 3). Since the tangent to the curve approaches infinitely close to zero, we determined a threshold of -0.0001, which corresponds to seven iterations for pigs and five iterations for chickens (Table 7 and Supplementary Table S11).
For both BWA-MEM and GTX, we confirmed increased mapping rate upon the final WZS alternative reference sequences, using sequencing data of an unrelated WZS sample (Supplementary Table S12). We also discovered newly mapped sequences and novel high-quality variants that can only be identified using our alternative reference sequence (Supplementary Table S13 and S14). By aligning the WZS alternative reference sequence with the original pig reference sequence, we identified 138,545 highly variable regions (HVR) recovered by our iterative strategy, among which 44.76% were overlapped with genes, and 33.59% were overlapped with CDSs (detailed in Table 8). The GO (Gene Ontology) analysis results showed that the genes overlapping with HVR were mostly enriched in sensory perception (Fig. 4), which is consistent with previous reports that the genes related to sensory perception experience rapid evolution 25 .  Figure 2. Iterative substitutions of the reference sequence enable the mapping of readX, which could not be mapped before due to consecutive mismatches. In the case of allowing one base mismatch, in the first iteration, read1 was mapped to the original reference genome REF, and the genome Alternative-REF1 was generated by base replacement; in the second iteration, read2 was mapped to Alternative-REF1, and the genome Alternative-REF2 was generated by base replacement; after two iterations, readX was mapped to Alternative-REF2, and the genome Alternative-REF3 was generated by base replacement. REF indicates the reference genome, and Alternative-REF1/Alternative-REF2/Alternative-REF3 indicates the alternative reference sequence in different iterative rounds. www.nature.com/scientificreports/  www.nature.com/scientificreports/

Discussion
Due to the rapid development of sequencing technologies, the volume of sequencing increases rapidly, while the price per base pair continues to fall. However, it is always important to sequence with optimal coverage to reach a balance between the cost and accuracy of the analysis. For low-coverage sequencing, some regions of the genome might not be covered, or heterozygote sites might be misgenotyped as homozygous 26 . Ultra-high coverage sequencing will, on the other hand, greatly increase the costs for both sequencing and the data analysis. Here, the optimal sequencing coverage of accurate genotype calls for pigs and chickens was 19× and 16×, respectively. This will provide a useful reference for related studies. However, obtaining accurate variant information not only depends on sequencing coverage but also on good reference genomes. Since the standard reference genome might not satisfy the special needs of different studies, in some studies, researchers have attempted to look for alternative reference genomes. For example, Ai H et al. studied the domestication of Chinese pigs. In their work, the genome sequences for WZS, instead of the published Duroc assembly, were used as the reference to ensure a better alignment with Chinese pigs 22 . Incarnato et al. created an alternative reference sequence for the E14 genome based on the mm9 assembly. The sequencing reads mapped to the E14 genome increased by around 5% compared to those of the mouse mm9 genome 27 . In this study, by using alternative reference sequences, the mapping rate increased by 0.61-1.68% for Chinese pig breeds and by 0.09-0.45% for Chinese chicken breeds. Compared to the chickens that were domesticated from red junglefowl in southeast Asia, the domestication of the pig took place independently in two locations: East Anatolia and China 28 . Our result thus suggests that the generation of alternative reference sequences is more necessary for species with complex genetic backgrounds. Unlike previous methods that used single step substitution, we employed an iterative strategy of generation to produce alternative reference sequences. Although here the GTX-one performs the sequence mapping and variant calling processes, we confirmed that this iterative strategy can cooperate with other mappers/callers (Supplementary Table S15). The alternative reference sequences generated by this approach enable the substitutions of consecutive mismatches in highly variable regions. These regions, where sequencing reads are not directly mappable, could cause a complete loss of information and leave a false impression that the region is highly conserved.
We found that several highly variable regions overlap with the genes or even CDSs, indicating a possible overestimation of the conversation of coding sequences during domestication. Our method corrected this problem and recovered phylogenetically informative sites that are missed by using public reference sequences, which could improve the accuracy of downstream analyses, including GWAS and genetic diversity evaluations. The GO enrichment analysis of the genes in highly variable regions were found to be enriched in sensory www.nature.com/scientificreports/ perception pathways, suggesting that, even though previous studies reported that the genes involved in sensory perception are among the most rapidly evolving genes 25 , this phenomenon might still be underestimated.
Our results indicate that alternative reference sequences are not only effective for improving the mapping rates of the breeds used to generate the alternative reference sequences but are also effective for other genetically close populations. Interestingly, we found that the variation count reported in the second round of iterations was significantly lower than that obtained during the first round, indicating a substantial number of fixed substitutions across distinct breeds. These fixed substitutions are not informative for either the target breed or other genetically close breeds and can be avoided when using an alternative reference sequence.
In summary, our iterative strategy for generating alternative reference sequences facilities the read alignment for genetically distant breeds and will improve all variant-based downstream analyses, especially for population genetics analyses.

Materials and methods
Samples. The whole genome sequencing data for 6 pigs, 7 chickens, 2 humans, and 1 chimpanzee were used in this study. These data were either downloaded from the NCBI SRA database (https ://www.ncbi.nlm.nih.gov/) and DDBJ (https ://ddbj.nig.ac.jp/DRASe arch/) or sequenced by our lab. The details of the samples are shown in Table 9.
Read mapping and variant calling. GTX-One by the Genetalks company, a commercially available FPGA-based hardware accelerator platform, was used in this study for both the read alignment and variant calling. The alignment process for GTX-One is accelerated by the FPGA implementation of the parallel seedand-extend approach based on the Smith-Waterman algorithm, while the variant calling process is accelerated by the FPGA implementation of the Genome Analysis Toolkit (GATK 3.7) 29 HaplotypeCaller (Pair-HMM) (see Supplementary Text for details). Based on the Genome in a Bottle Consortium (GIAB) gold standard callset, we evaluated the performance of both the GTX and the BWA-GATK "Best Practices" workflow (GBP, BWA 0.7.17 + GATK 3.7) in terms of precision and sensitivity, using the following formulas: precision = TP/(TP + FP) and sensitivity = TP/(TP + FN), where TP represents true positive (variants called with the same genotype as the www.nature.com/scientificreports/ gold standard callset), FP represents false positive (variants called but not in the gold standard callset) and FN represents false negative (variants in the gold standard callset but not called).
Determining the optimal genome sequencing coverage. Raw reads were randomly sampled at 1× to 30× genome coverage with an increment of 1×. For each coverage, the sequencing data were aligned back to the reference genome, and variations were called using GTX-One. Variant counts, mapping rates (defined as the ratio of mapped reads to the total reads), and coverage ratios (defined as the proportion of the loci at a coverage greater than or equal to 1 compared to the reference genome) were recorded. Variant sets reported by the dataset to have excessive coverage (40× or higher) were used as gold standard. The BCFtools isec function 30 was used to compare the gold standard VCFs for each sequence's coverage, and we reported the sensitivity using the formula sensitivity = TP/(TP + FN), where TP represents a true positive (variants called with the same genotype as the gold standard callset), and FN represents a false negative (variants in the gold standard callset but not called). CurveExpert1.4 (https ://www.curve exper t.net) was used to fit the variant counts against the coverage. We determined the optimal genome sequencing coverage by selecting the tangent slope of the curve at a threshold of 0.0001 using an in-house python script.
Iterative strategy for the generation of alternative reference sequences. Given a specific variant callset (in VCF format), the GATK FastaAlternateReferenceMaker was used to replace the reference bases with the variations recorded in the callset to generate a new reference sequence. For each round, the variant callset was generated by GTX-One against the previous reference sequence. The reference sequence was continually updated during the iterative process. We selected WZS and YJ as examples for pigs and chickens, respectively, and, for both species, the iterative process continued for 30 rounds. Variant counts, mapping rates, and coverage ratios were recorded. Similarly, CurveExpert1.4 was used to fit the variant counts against the number of iterations to determine the optimal number of iterations. After the optimal number of iterations was determined, the final alternative reference sequence for the target breed was reported.
Whole genome alignments for identifying highly variable regions and genome similarities. LASTZ 31 was used to align the final alternative reference sequences and original reference sequences chromosome by chromosome. By parsing the alignment output of LASTZ, we reported highly variable regions that contain three or more consecutive mismatches. The gene-based annotations for highly variable regions were produced using the coordinate information in the genomic GTF file from the Ensembl repository. Functionbased annotations were then based on gene-based annotations using Gene ontologies via PANTHER (https ://panth erdb.org/). We also used Minimap2 32 to align the Duroc and the WZS alternative reference sequences from each iteration against the public WZS genome. The genomic similarities were calculated as the number of matched bases divided by the total genome size according to the Minimap2 output.