SpeedSeq: Ultra-fast personal genome analysis and interpretation

SpeedSeq is an open-source genome analysis platform that accomplishes alignment, variant detection and functional annotation of a 50× human genome in 13 hours on a low-cost server, alleviating a bioinformatics bottleneck that typically demands weeks of computation with extensive hands-on expert involvement. SpeedSeq offers competitive or superior performance to current methods for detecting germline and somatic single nucleotide variants, indels, and structural variants, and includes novel functionality for streamlined interpretation.

Our implementations of FreeBayes and CNVnator are also faster than their original counterparts. SpeedSeq's parallelized FreeBayes produced SNV and indel calls in 1.4 hours, compared to 19.9 hours for the single-threaded version and 2.2 hours with naïve parallelization by chromosome. CNVnator took only 0.5 hours, compared with 13.5 hours using the original (version 0.3), a 27-fold speed improvement.
Our parallelized software implementations demonstrate little or no difference from the original versions. SAMBLASTER produces nearly identical output to Picard Tools duplicate marking 4 . Our parallelized FreeBayes output differs from the single-threaded standard version at 1,310/4,333,589 sites (0.03%), with only 43 of these variants validating against the GIAB truth set. In our implementation of CNVnator, there were minor differences in output due to the use of a lookup table and a slight loss of floating-point precision. Copy number variable regions differed by ~2,300 bp over the entire genome and often corresponded to ambiguous CNV calls with low-amplitude copy number deviations and low confidence scores.

Supplementary Note 2: Characteristics of miscalled variants
Non-unique 100-mers harbored higher percentages of false positive and false negative SNVs, and both SNVs and indels showed increased false positive rates in satellite repeats (Supplementary Table 1). We were not able to evaluate the role of segmental duplications and simple repeats in variant calling accuracy because the Genome in a Bottle truth set excluded these features, however, we expect performance to be similarly diminished in these regions.
Clusters of nearby variants can complicate variant detection due to misalignment, and indeed false positives were enriched in their proximity to indels, with 16.9% of false positive SNVs located within 10 bp of an indel compared to 0.06% of true positives. Upon manual inspection, many of these resulted from complex variants where flanking indels caused SpeedSeq to call phantom SNVs. A smaller portion of miscalls was due to different representations of the same variant, despite our adherence to the variant normalization protocols used by the Genome in a Bottle Consortium 5 . This reflects a persistent deficiency in tools to compare VCF files, and suggests that our metrics likely underestimate SpeedSeq's true performance.
To characterize the nature of miscalled SVs, we used their composite long-read/1KGP validation status to assess whether they were true (TP) or false positives (FP), and counted the number of SVs in which either of the breakends overlapped with annotated genomic features. We expressed this count as a percentage of the number of true positives (N=3,388) or false positives (N=3,088) (Supplementary Table 1). False positives overlapped more frequently with simple repeats (TP: 10.6%, FP: 29.4%) and segmental duplications (TP: 4.5%, FP: 23.7%). However, we note the important caveat that SVs in repetitive regions of the genome are also less likely to validate by PacBio or Illumina Moleculo long-reads, and may have been under-ascertained by 1KGP SV mapping efforts, and therefore some of this enrichment may be due to technical deficiencies in validation (Supplementary Table 1).

Supplementary Note 3: De novo variant detection
Identification of de novo variants from a family trio is a common WGS application and a natural extension of SpeedSeq's multi-sample variant detection functionality. To evaluate the performance of our tools, we defined a "trio" from the CEPH 1463 pedigree in which NA12878 represented the child, NA12878's biological mother (NA12892) represented one parent, and NA12877 (an unrelated individual) represented her father (Supplementary Fig. 4). Thus, the truth set of 288,409 de novo SNVs were defined by those variants that were present in the NA12878 GIAB 2.17 variant truth set but absent from NA12877 and NA12892, as determined by an independent variant set from Real Time Genomics 6 . Using the default minimum variant quality score of 1 and no additional filtering, SpeedSeq detected 98.3% of the "de novo mutations" with a false discovery rate of 5.9%. These putative mutation calls can be further refined within GEMINI (or using ad hoc scripts) by annotating known polymorphisms (e.g., dbSNP, 1KGP, NHLBI-ESP, ExAC), and by filtering for known sources of false positives. For example, simply requiring all three samples to have a minimum read-depth of 30 at putative mutation sites -a built-in filtering parameter in GEMINI -reduced FDR to 3.2%, with a minor effect on sensitivity (97.1%). Functionally relevant de novo mutations can be prioritized by filtering on diverse annotations and variant impact scores. Of course, reducing false positives to very low levelsas required for mutation rate studies -is a more difficult problem that requires extensive probabilistic filtering 7,8 .

Supplementary Note 4: Excluded regions
Despite the high quality of the human reference genome, artifacts remain in low-complexity regions and unannotated paralogous sequences that delay processing time and confound variant interpretation. Thus, we have excluded a static set of high-depth regions in the human genome from SNV, indel, and SV breakpoint calling modules of SpeedSeq. These regions exhibit an aberrant increase in sequencing coverage depth where reads from disparate parts of the genome accumulate, violating the diploid assumption of downstream variant calling algorithms. Most high-depth regions are caused by misassembled regions of the reference genome, where moderate or high copy repeats have been collapsed into a single (or few) copies, causing read pile-ups.
To identify these high-depth regions, we aligned reads from the CEPH 1463 pedigree to GRCh37 with BWA-MEM and measured aggregate coverage depth from all 17 family members plus one replicate. We then excluded the 15.6 Mb (0.6% of the genome) where the depth was greater than 2 × mode coverage + 3 standard deviations (Supplementary Fig. 4). While this static set of excluded regions was based on a single large family, the depth cutoff accommodates a 2-fold increase in copy number relative to the reference genome -corresponding to a homozygous duplication (4 copies) in a genomic region that is single copy in the reference -minimizing the bias toward fixed polymorphism in that family. In the future, we envision that excluded genomic regions will be defined by a much larger set of high-quality human genomes from diverse populations. For structural variation detection with LUMPY, we also exclude the mitochondrial genome, which is prone to false positive calls due to extremely high depth. These regions are not excluded from read-depth analysis using CNVnator, so it is possible (albeit technically complicated) to detect CNVs in these regions using SpeedSeq.
We tested for feature enrichment in repetitive elements from RepeatMasker and mappability from the CRG Alignability track 9,10 . As expected, the excluded regions are highly enriched in repetitive regions of the genome that are known to be poorly assembled in the reference genome including satellite repeats (54-fold enrichment), segmental duplications (9-fold), and sequences near centromeres (28-fold) and other assembly gaps (25-fold) (Supplementary Table 5