Jointly aligning a group of DNA reads improves accuracy of identifying large deletions

Abstract

Performing sequence alignment to identify structural variants, such as large deletions, from genome sequencing data is a fundamental task, but current methods are far from perfect. The current practice is to independently align each DNA read to a reference genome. We show that the propensity of genomic rearrangements to accumulate in repeat-rich regions imposes severe ambiguities in these alignments, and consequently on the variant calls—with current read lengths, this affects more than one third of known large deletions in the C. Venter genome. We present a method to jointly align reads to a genome, whereby alignment ambiguity of one read can be disambiguated by other reads. We show this leads to a significant improvement in the accuracy of identifying large deletions (≥20 bases), while imposing minimal computational overhead and maintaining an overall running time that is at par with current tools. A software implementation is available as an open-source Python program called JRA at https://bitbucket.org/jointreadalignment/jra-src.

Variant calling tool evaluation for variable size indel calling from next generation whole genome and targeted sequencing data

“…HuRef data also allowed us to reproduce a realistic dataset that would capture the challenges of small indel and large indel calling in human genome. If the indels were to be inserted randomly across the genome it would underestimate the proportion of indels located in ambiguous regions, where the indels may be represented in different positions [39] (S1 and S2 Figs). Here, we reconstructed chromosome 1 of the HuRef genome, based on human reference genome hg19, by similar methods to [39,40] and inserted indels of a real human individual into the corresponding position in reference genome (S1 File).…”

“…If the indels were to be inserted randomly across the genome it would underestimate the proportion of indels located in ambiguous regions, where the indels may be represented in different positions [39] (S1 and S2 Figs). Here, we reconstructed chromosome 1 of the HuRef genome, based on human reference genome hg19, by similar methods to [39,40] and inserted indels of a real human individual into the corresponding position in reference genome (S1 File). In addition, we reconstructed chromosome 1 with two different haplotypes by randomly selecting variants from different size ranges and only inserting them into one of the haplotypes as heterozygous variants or into both haplotypes as homozygous variants.…”

“…A possible reason for this is that variants in repeat regions can be called by CHM1 cell line PacBio calls, but cannot be detected efficiently by Illumina short-read sequencing [47]. A previous study which used CHM1 cell line dataset as a truth set, also suggested that the truth set might be missing variants [39]. The results were separated into deletions (50bp -500bp and > 500bp) and insertions (50bp -500bp and > 500bp).…”

“…For the CHM1 cell line WGS datasets, we defined a TP as a tool-detected, filter-passed indel which has at least 20% overlap with a truth-indel via BEDtools [46] (S1 File). Criteria was kept loose because previous studies [39,47] have reported that the variant calling methodology for CHM1 cell line WGS datasets has its own bias, which leads to a low concordance between tool-detected indels and truth-indels (S1 File).…”

“…For 5 coverage sequencing data, the low read depth in indel breakpoints is the main reason that caused FPs. SR regions in genome may cause misalignment of reads and lead to further inaccurate indel calling, such as wrong positions or sizes due to the breakpoint ambiguity [39]. The breakpoint ambiguity of indels in SR regions may also cause incorrect allele frequency counting, thus leading to false genotype-level results ( S9 Fig).…”