近年來由於次世代定序技術(NGS)的出現，使得我們可以同時對大量DNA進行定序，大幅降低定序成本，許多相關的研究也跟著大幅擴展，NGS數據中變異辨認幾乎是所有下游分析和解釋過程都依賴的關鍵步驟，但是對於準確基因變異體偵測仍然有許多問題及挑戰。因此我們先前提出過一個工具Explicit Alternative Genome Likelihood Evaluator(EAGLE)，透過使用顯示概率模型，處理變異辨認過程中固有的不確定性。明確評估測序數據與推定變異所隱含的替代基因組序列的匹配程度。但是不論是EAGLE的模型抑或是一般的變異辨認方法，或多或少都仰賴將讀序映射到參考基因組序列上的結果。而映射的方法是透過將讀序與參考基因組序列比對，根據比對的相似程度決定他們可能來自參考基因組序列上的位置。

由於我們比對的目標參考基因組序列是線性的，僅包含單個序列，因此無法捕獲基因組的多樣性。所以比對會出現一定程度的誤差，造成我們讀序可能被映射到錯誤的位置或是無法映射，這種現象我們稱之參考偏差。錯誤的讀序映射會導致假陰性或假陽性變異辨認，進一步影響到我們對於辨認基因體變異的準確性。而先前已經有許多研究討論消除參考偏差所帶來的影響，其中一種方法就是透過改進參考基因組序列，使其能夠包含基因組的多樣性，從而達到更好的映射結果。本文中我們提出一個類似概念的方法，並應用於我們先前提出的工具-EAGLE上，探討我們提出的方法對於降低參考偏差帶來的影響是否可行以及其潛力。

我們的方法可以分成五個部分，第一個部分受透過基因突變資料(VCF)中取得變異點資訊，結合參考基因組序列建立假設序列。第二個部分是為讀序檔案(FASTQ)建立一個能夠快速搜尋的索引結構。第三個部分是透過Burrows-Wheeler Aligner(BWA)以及我們上一部份建立的讀序索引結構在讀序檔案中找出所有能夠與我們建立的假設序列匹配的讀序。第四個部分則是檢查原本的映射結果檔案(BAM)中的堆積序列，堆積序列是在映射結果中被映射到相同區域的所有讀序，檢查是否存在原有的堆積序列中，對那些不存在堆積序列中的讀序建立相關的比對資訊。第五部分則是將我們找到的讀序以及比對資訊加入堆積序列中，改進我們EAGLE的計算結果。

在實驗中，我們模擬參考偏差發生的情景，透過我們的方法驗證是否能找到那些受到參考偏差的影響而遺失的讀序，結果顯示我們能夠找到當中大部分被遺失的讀序，但是在一些重複序列過多的區域可能有所限制。我們同時也在 dbSNP的變異資料及上進行測試，基本上也能夠取得相似的結果，惟在單核苷酸多態性的變異上效果較小。但大體而言，綜合我們所增加的時間以及成效，我們認為我們的方在在解決參考偏差的影響是相當有潛力的，但是應用到實際情況還是需要更進一步的驗證以及謹慎的解釋。

In recent years, due to the emergence of next-generation sequencing technology (NGS), we can sequence a large amount of DNA at the same time, greatly reducing the cost of sequencing, and expanding its applications. Genome variant calling from NGS data is a key step for many applications in biology and medicine, but there are still many problems and challenges for accurate variant calling. Most variant calling procedures start with the so called ``pile-up'', an alignment of reads to the reference genome around the position of a putative variant. However this alignment cannot be expected to be certain; because the human genome is repetitive (many distinct regions are be similar to each other) and the reference genome sequence differs from any individual. Moreover alignment only to the reference genome naturally tends to favor the reference sequence over others (so called reference bias).

Reference bias is a widely recognized problem in variant calling. One research direction is to represent the reference genome with a graph instead of a simple string to include some common variants. In a different approach, we previously addressed reference bias by proposing a method (EAGLE: Explicit Alternative Gene Likelihood Evaluator) which locally aligns reads not only to the reference sequence, but also to putative variant sequences; using an explicit probability model to combine evidence from all possible local alignments. The results were encouraging but fell short of completely eliminating reference bias, because that version of EAGLE still relies on read mapping to identify the reads which should be considered when evaluating a variant.

To completely eliminate reference bias in putative variant evaluation, any read similar to the variant sequence should be considered, even it does not match the reference sequence. Na\"ively matching all reads to all variants would be far too slow, but in principle the search indexes used to accelerate standard read mapping could be adapted to search reads against a variant sequence query. But one might ask 1) if this would be feasible in practice (in terms of computer memory and time) and 2) if it would really make a difference in the variant calling results. Last year we conducted a feasibility study which suggested that constructing a search index on reads is feasible for whole exome sequencing and that searching such a read index can indeed find relevant reads not present in the ``pile-up''. That study however was only a feasibility study and did not integrate the read index with the EAGLE variant likelihood computation.

Here we explore the read index idea more thoroughly and integrate it into the EAGLE software.
Our method can be divided into two preprocessing steps and then four steps for each putative variant. In the first pre-processing step we align the reads to the reference genome. In the second pre-processing step, we create a Burrows-Wheeler Aligner (BWA) index structure for the read sequence data (FASTQ format file) that can be quickly queried. Then for each variant in a list of putative variants (VCF format file), we first edit the reference genome to include the putative variant, then use that to query the read index structure, merge the matching reads with any reads mapping to the reference genome near the putative variant genome position, and finally apply the EAGLE probabilistic model computation to the merged set of reads to obtain an unbiased likelihood ratio between variant and reference sequence for that position.

We performed some experiments to evaluate the method. In one experiment we use real sequence read data but doctor the reference genome to produce a false reference genome sequence, thus simulating a situation in which we know the position of differences between individual and reference. Then we map reads to this false genome sequence to obtain a corresponding false pile-up. We tested if our read index can find reads supporting the actual genome sequence, but missing in the false pile-up. The experimental results show that we can find most of the missing reads, but have some difficulty in regions with repetitive sequences. We also tested real putative variants from the dbSNP dataset, observing qualitatively similar results. We show that for putative indel variants, the extra reads found by the read index have a large effect on the EAGLE likelihood ratio (and in the correct direction in cases where we know what the answer should be). The likelihood ratio of Single Nucleotide Variants (SNV)s, on the other hand is not greatly affected, a reasonable outcome since standard read mapping can tolerate single mismatchs, largely obviating the need to use a special read index. More comprehensive benchmarking is left for future work.

[1] Rasmus Nielsen, Joshua S Paul, Anders Albrechtsen, and Yun S Song. Genotypeand snp calling from next-generation sequencing data.Nature Reviews Genetics,12(6):443–451, 2011.
[2] Ryan Poplin, Valentin Ruano-Rubio, Mark A DePristo, Tim J Fennell, Mauricio OCarneiro, Geraldine A Van der Auwera, David E Kling, Laura D Gauthier, AmiLevy-Moonshine, David Roazen, et al. Scaling accurate genetic variant discoveryto tens of thousands of samples.BioRxiv, page 201178, 2018.
[3] Christopher T Saunders, Wendy SW Wong, Sajani Swamy, Jennifer Becq, Lisa JMurray, and R Keira Cheetham. Strelka: accurate somatic small-variant callingfrom sequenced tumor–normal sample pairs.Bioinformatics, 28(14):1811–1817,2012.
[4] Erik Garrison and Gabor Marth. Haplotype-based variant detection from short-read sequencing.arXiv preprint arXiv:1207.3907, 2012.
[5] Ryan Poplin, Pi-Chuan Chang, David Alexander, Scott Schwartz, ThomasColthurst, Alexander Ku, Dan Newburger, Jojo Dijamco, Nam Nguyen, Pegah TAfshar, et al. A universal snp and small-indel variant caller using deep neuralnetworks.Nature biotechnology, 36(10):983–987, 2018.
[6] Daniel C Koboldt. Best practices for variant calling in clinical sequencing.GenomeMedicine, 12(1):1–13, 2020.
[7] Xiaopeng Bian, Bin Zhu, Mingyi Wang, Ying Hu, Qingrong Chen, Cu Nguyen,Belynda Hicks, and Daoud Meerzaman. Comparing the performance of selected variant callers using synthetic data and genome segmentation.BMC bioinformat-ics, 19(1):1–11, 2018.
[8] Jiayun Chen, Xingsong Li, Hongbin Zhong, Yuhuan Meng, and Hongli Du. Sys-tematic comparison of germline variant calling pipelines cross multiple next-generation sequencers.Scientific reports, 9(1):1–13, 2019.
[9] Sohyun Hwang, Eiru Kim, Insuk Lee, and Edward M Marcotte. Systematic com-parison of variant calling pipelines using gold standard personal exome variants.Scientific reports, 5(1):1–8, 2015.
[10] Sen Zhao, Oleg Agafonov, Abdulrahman Azab, Tomasz Stokowy, and EivindHovig. Accuracy and efficiency of germline variant calling pipelines for humangenome data.Scientific reports, 10(1):1–12, 2020.
[11] Rui Tian, Malay K Basu, and Emidio Capriotti. Computational methods andresources for the interpretation of genomic variants in cancer.BMC genomics,16(8):1–19, 2015.
[12] Adam Cornish and Chittibabu Guda. A comparison of variant calling pipelinesusing genome in a bottle as a reference.BioMed research international, 2015, 2015.
[13] Jennifer K Sehn. Insertions and deletions (indels). InClinical genomics, pages129–150. Elsevier, 2015.
[14] Stephen B Montgomery, David L Goode, Erika Kvikstad, Cornelis A Albers,Zhengdong D Zhang, Xinmeng Jasmine Mu, Guruprasad Ananda, Bryan Howie,Konrad J Karczewski, Kevin S Smith, et al. The origin, evolution, and functional impact of short insertion–deletion variants identified in 179 human genomes.Genome research, 23(5):749–761, 2013.
[15] Tony Kuo, Martin C Frith, Jun Sese, and Paul Horton. Eagle: explicit alternativegenome likelihood evaluator.BMC medical genomics, 11(2):1–10, 2018.
[16] Vitor Sousa and Jody Hey. Understanding the origin of species with genome-scaledata: modelling gene flow.Nature Reviews Genetics, 14(6):404–414, 2013.
[17] Paolo Ferragina and Giovanni Manzini. Indexing compressed text.Journal of theACM (JACM), 52(4):552–581, 2005.
[18] Michael Burrows and David Wheeler. A block-sorting lossless data compressionalgorithm. InDigital SRC Research Report. Citeseer, 1994.
[19] Rui Martiniano, Erik Garrison, Eppie R Jones, Andrea Manica, and RichardDurbin. Removing reference bias and improving indel calling in ancient dna dataanalysis by mapping to a sequence variation graph.Genome biology, 21(1):1–18,2020.
[20] Nae-Chyun Chen, Brad Solomon, Taher Mun, Sheila Iyer, and Ben Langmead.Reference flow: reducing reference bias using multiple population genomes.Genome biology, 22(1):1–17, 2021.
[21] Torsten G ̈unther and Carl Nettelblad. The presence and impact of reference biason population genomic studies of prehistoric human populations.PLoS genetics,15(7):e1008302, 2019.
[22] Jessica Lau.Reference bias:Challenges and solutions.https://www.sevenbridges.com/reference-bias-challenges-and-solutions/.
[23] Erik Garrison, Jouni Sir ́en, Adam M Novak, Glenn Hickey, Jordan M Eizenga,Eric T Dawson, William Jones, Shilpa Garg, Charles Markello, Michael F Lin,et al. Variation graph toolkit improves read mapping by representing geneticvariation in the reference.Nature biotechnology, 36(9):875–879, 2018.
[24] Heng Li, Xiaowen Feng, and Chong Chu. The design and construction of referencepangenome graphs with minigraph.Genome biology, 21(1):1–19, 2020.