Read mapping and variant calling in human short-read DNA sequences

1. Read mapping and variant calling in human short-read DNA sequences Gabor T. Marth Boston College Biology Department 1000 Genomes Meeting Cold Spring Harbor Laboratory May 5-6. 2008

2. 1000 Genomes – related work • Software • Read alignment / assembly (Michael Stromberg) • SNP / short-INDEL calling (Gabor Marth) • Structural Variation calling (Chip Stewart) • Read simulator (Weichun Huang) • Benchmarking suite (Weichun Huang) • Read mapping based studies • Read accuracy / quality value analysis • Read simulations • Variant calling based study • SNP discovery: sample size / read coverage (Aaron Quinlan)

3. MOSAIK – sequence aligner/assembler Poster thumb Michael Strömberg (see poster at Genome Meeting) • maps reads to reference: short-hash based scan + Smith-Waterman alignment

4. MOSAIK – Features • produces gapped alignments essential for tolerating reads with insertion-deletion read errors and short insertion-deletion alleles • adapted to all available NGS platforms • can create “mixed” alignments of reads from different platforms (except SOLiD color-space reads)

5. MOSAIK – Resolving multiple map locations

6. MOSAIK – Performance Uses a lot of RAM for mammalian alignments – precomputed file based versions are available for RAM-limited users Run dissection (timeline figure from Michael)

7. MOSAIK – Accuracy

8. Erroricity – read accuracy / quality values Motivation: why does read accuracy matter? Why does quality value accuracy matter? we are using quality values to distinguish between sequencing error and true allelic difference Q-values should correspond well with actual sequencing error rate

9. Erroricity – study design & pipeline • Sampled 3 lanes each from 3 runs • Aligned reads with MosaikPE (up to 4 mismatches), keeping only consistently mapping pairs • Looked at read-specific, position-specific error rates • Compared SUB, IN, and DEL error rates • Looked at overall quality value vs. measured error rate, and position specific quality value vs. measured error rate • Compared the first and the second end-reads of read pairs • Compared RAW vs. CALIBRATED Q-values Derek Barnett

10. Read simulations (Weichun Huang, Aaron) • Input • Conceptual schema of read simulation • Representational biases (GC-driven and others…) [Chip] • Error and Q-value generation: 2D tables of read position, assigned Q-value  true Q-value, frequency • Speed / RAM / space • Data output • Software benchmarking system Weichun Huang, see poster at Genome Meeting

11. GigaBayes: SNP and short-INDEL calling • The new GigaBayes program: pop-gen and diploid priors, trio priors • Speed • Input / output behavior • Bayesian math focused on the individual genotype • How to deal with multiple reads from a single individual • Diploid individual • Multiple diploid individuals • Trio members • Prior frequency of an allele • Taking into account Q-value for allele • What is needed to call an allele? (# reads, Q, # people)

12. Variant calling (SNPs and short-INDELs) individuals fragments population reads G1 aacgtCaggct aacgtCaggct aacgtTaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct G2 aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtCaggct aacgtTaggct aacgtTaggct aacgtCaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct G3 aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtCaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtTaggct aacgtCaggct aacgtTaggct aacgtCaggct quality values priors sampling likelihood

13. Bayesian variant detection math Priors: (1) based on all the individuals from which reads are aligned. (2) Theta, P(AF=i), specific diploid genotype layout given AF=I Which of the two chromosomes a read represents? Calculated from multinomial (or binomial distribution) P(base call is an B | read template is T) – comes from quality values

14. SNP calling and genotyping P(SNP) = total probability of all non-monomorphic genotype combinations P(Gi) = marginal probability consequence: data from other individuals influence the genotype call of a given individual: include illustration using testProb program in GigaBayes package.

15. Variant calling tests in simulated data Aaron Quinlan (see poster at the Genome Meeting)

16. Variant calling – Estimated vs. population AF

17. Variant calling – AF (cont’d)

18. Variant calling – SNP discovery sensitivity

19. Variant calling – Genotype completeness

20. Variant calling – Genotype completeness

21. Summary / Conclusions

22. Thanks Chip Michael Derek Aaron Weichun

23. MOSAIK (Michael Stromberg) • MOSAIK is a reference-sequence guided aligner / assembler replace this with an animated figure illustrating mapping against reference and padding, by moving / stretching bases in the reads and in the reference sequence

24. MOSAIK – Features and characteristics • aligns reads to genome (higher RAM usage but also many desirable consequences) • offers several algorithmic levels that trade off between speed and accuracy • able to report “every” decent alternative map location for sequence reads, and distinguishes between uniquely and non-uniquely mapped reads • designed to work with all currently available technologies (Illumina, 454, AB, Helicos) and to include mixed read sets into a single anchored “assembly” • PE-aware • recently scaled up to mammalian alignments

25. Structural variation discovery • copy number variations (deletions & amplifications) can be detected from variations in the depth of read coverage • structural rearrangements (inversions and translocations) require paired-end read data

26. Software evaluation suite

27. GigaBayes

28. Read length and throughput Illumina/Solexa, AB/SOLiD short-read sequencers 1Gb (1-4 Gb in 25-50 bp reads) bases per machine run 100 Mb 454 pyrosequencer (20-100 Mb in 100-250 bp reads) 10 Mb ABI capillary sequencer 1Mb read length 10 bp 100 bp 1,000 bp

29. Current and future application areas Genome re-sequencing: somatic mutation detection, organismal SNP discovery, mutational profiling, structural variation discovery reference genome DEL SNP De novo genome sequencing • Short-read sequencing will be (at least) an alternative to micro-arrays for: • DNA-protein interaction analysis (CHiP-Seq) • novel transcript discovery • quantification of gene expression • epigenetic analysis (methylation profiling)

30. 3. Alignment of billions of reads Fundamental informatics challenges 1. Interpreting machine readouts – base calling, base error estimation 2. Dealing with non-uniqueness in the genome: resequenceability

31. Informatics challenges (cont’d) 4. SNP and short INDEL, and structural variation discovery 5. Data visualization 6. Data storage & management

32. Challenge 1. Base accuracy and base calling • machine read-outs are quite different • read length, read accuracy, and sequencing error profiles are variable (and change rapidly as machine hardware, chemistry, optics, and noise filtering improves) • what is the instrument-specific error profile? • are the base quality values satisfactory? • (1) are base quality values accurate? • (2) are most called bases high-quality?

33. 454 pyrosequencer error profile • multiple bases in a homo-polymeric run are incorporated in a single incorporation test  the number of bases must be determined from a single scalar signal the majority of errors are INDELs • error rates are nucleotide-dependent

34. 454 base quality values • the native 454 base caller assigns too low base quality values

35. PYROBAYES: determine base number New 454 base caller: data likelihoods priors posterior base number probability

36. PYROBAYES: base calls and quality values • call the most likely number of nucleotides • produce three base quality values: • QS (substitution) • QI (insertion) • QD (deletion)

37. PYROBAYES: Performance • better correlation between assigned and measured quality values • higher fraction of high-quality bases

38. Illumina/Solexa base accuracy • error rate grows as a function of base position within the read • a large fraction of the reads contains 1 or 2 errors

39. Illumina/Solexa base accuracy (cont’d) • Actual base accuracy for a fixed base quality value is a function of base position within the read (i.e. there is need for quality value calibration) • Most errors are substitutions  PHRED quality values work

40. RepeatMasker does not capture all micro-repeats, i.e. repeats at the scale of the read length • Near-perfect micro-repeats can be also a problem because we want to align reads even with a few sequencing errors and / or SNPs Challenge 2. Resequenceability • Reads from repeats cannot be uniquely mapped back to their true region of origin

41. Repeats at the fragment level “base masking” “fragment masking”

42. Fragment level repeat annotation • bases in repetitive fragments may be resequenced with reads representing other, unique fragments  fragment-level repeat annotations spare a higher fraction of the genome than base-level repeat masking

43. Find perfect and near-perfect micro-repeats • Hash based methods (fast but only work out to a couple of mismatches) • Exact methods (very slow but find every repeat copy) • Heuristic methods (fast but miss a fraction of the repeats)

44. Challenge 3. Read alignment and assembly • resequencing requires reference sequence-guided read alignment • to align billions of reads the aligner has to be fast and efficient • INDEL errors require gapped alignment • individually aligned reads must be “assembled” together • has to work for every read type (short, medium-length, and long reads) • must tolerate sequencing errors and SNPs • must work with both base-level and fragment-level repeat annotations • transcribed sequences require additional features e.g. splice-site aware alignment capability • most frequently used tools: BLAT (only pair-wise), SSAHA (pair-wise), MAQ (pair-wise and assembly), ELAND (pair-wise), MOSAIK (pair-wise and assembly, gapped)

45. ABI/cap. 454/FLX 454/GS20 Illumina MOSAIK: co-assembling different read types

46. Challenge 4. Polymorphism discovery • shallow and deep read coverage • most candidates will never be “checked” only very low error rates are acceptable • we updated PolyBayes to deal with new read types • made the new software (PBSHORT) much more efficient

47. Challenge 5. Data visualization • aid software development: integration of trace data viewing, fast navigation, zooming/panning • facilitate data validation (e.g. SNP validation): co-viewing of multiple read types, quality value displays • promote hypothesis generation: integration of annotation tracks

48. Challenge 6. Massive data volumes • two connected working groups to define standard data formats Short-read format working group ssrformat@ubc.ca (Asim Siddiqui, UBC) Assembly format working group Boston College http://assembly.bc.edu

49. Next-generation sequencing software Machine manufacturers’ sites plus third-party developers’ sites, e.g.: http://sourceforge.net/projects/maq/ http://bioinformatics.bc.edu/marthlab/Mosaik

50. Applications in various discovery projects 1. SNP discovery in shallow, single-read 454 coverage (Drosophila melanogaster) 2. Mutational profiling in deep 454 data (Pichia stipitis) (image from Nature Biotech.) 3. SNP and INDEL discovery in deep Illumina / Solexa short-read coverage (Caenorhabditis elegans)