600 likes | 744 Views
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. 1000 Genomes – related work. Software Read alignment / assembly (Michael Stromberg)
E N D
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
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)
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
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)
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)
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
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
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
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)
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
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
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.
Variant calling tests in simulated data Aaron Quinlan (see poster at the Genome Meeting)
Thanks Chip Michael Derek Aaron Weichun
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
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
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
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
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)
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
Informatics challenges (cont’d) 4. SNP and short INDEL, and structural variation discovery 5. Data visualization 6. Data storage & management
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?
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
454 base quality values • the native 454 base caller assigns too low base quality values
PYROBAYES: determine base number New 454 base caller: data likelihoods priors posterior base number probability
PYROBAYES: base calls and quality values • call the most likely number of nucleotides • produce three base quality values: • QS (substitution) • QI (insertion) • QD (deletion)
PYROBAYES: Performance • better correlation between assigned and measured quality values • higher fraction of high-quality bases
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
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
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
Repeats at the fragment level “base masking” “fragment masking”
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
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)
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)
ABI/cap. 454/FLX 454/GS20 Illumina MOSAIK: co-assembling different read types
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
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
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
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
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)