Linkage Disequilibrium Based SNP Genotype Calling from Short Sequencing Reads Ion Mandoiu Computer Science and Engineering Department University of Connecticut Joint work with S.
Download ReportTranscript Linkage Disequilibrium Based SNP Genotype Calling from Short Sequencing Reads Ion Mandoiu Computer Science and Engineering Department University of Connecticut Joint work with S.
Linkage Disequilibrium Based SNP Genotype Calling from Short Sequencing Reads Ion Mandoiu Computer Science and Engineering Department University of Connecticut Joint work with S. Dinakar, J. Duitama, Y. Hernández, J. Kennedy, and Y. Wu Ultra-High Throughput Sequencing Recent massively parallel sequencing technologies deliver orders of magnitude higher throughput compared to classic Sanger sequencing Roche/454 FLX Titanium 400bp reads 400Mb/10h run ABI SOLiD 2.0 25-35bp reads 3-4Gb/6 day run Illumina Genome Analyzer II 35-50bp reads 1.5Gb/2.5 day run Helicos HeliScope 25-55bp reads >1Gb/day Personal Genomes: The Future is Now! $100,000,000 C.Venter [email protected] $10,000,000 Cost $1,000,000 $100,000 J. Watson [email protected] NA18507 Illumina@36x SOLiD@12x $10,000 $1,000 $100 days weeks months Sequencing Time years Challenges for Genomic Medicine at Single-Base Resolution Medical sequencing focuses on genetic variation (SNPs, CNVs, genome rearrangements) Requires accurate determination of both alleles at variable loci This is limited by coverage depth due to random nature of shotgun sequencing For the Venter and Watson genomes (both sequenced at ~7.5x average coverage), comparison with SNP genotyping chips has shown only ~75% accuracy for sequencing based calls of heterozygous SNPs [Levy et al 07, Wheeler et al 08] [Wendl&Wilson 08] predict that 21x coverage is required for sequencing of normal tissue samples based on idealized theory that “neglects any heuristic inputs” What heuristic inputs help? How much can we gain from improved data analysis? Linkage Disequilibrium: Sources & Modeling HMM model of haplotype frequencies F1 F2 H1 H2 … Fn Hn Fi = founder haplotype at locus i, Hi = observed allele at locus i P(Fi), P(Fi | Fi-1) and P(Hi | Fi) estimated from reference panel such as Hapmap For given haplotype h with n SNPs, P(H=h|M) can be computed in O(nK2) using forward algorithm, where K=#founders Pipeline for LD-Based Genotype Calling Reference genome sequence Read sequences >gnl|ti|1779718824 name:EI1W3PE02ILQXT TCAGTGAGGGTTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTTGAGACAGAATTTTGCTCTT >gnl|ti|1779718824 name:EI1W3PE02ILQXT GTCGCCCAGGCTGGTGTGCAGTGGTGCAACCTCAGCTCACTGCAACCTCTGCCTCCAGGTTCAAGCAATT TCAGTGAGGGTTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTTGAGACAGAATTTTGCTCTT CTCTGCCTCAGCCTCCCAAGTAGCTGGGATTACAGGCGGGCGCCACCACGCCCAGCTAATTTTGTATTGT GTCGCCCAGGCTGGTGTGCAGTGGTGCAACCTCAGCTCACTGCAACCTCTGCCTCCAGGTTCAAGCAATT >gnl|ti|1779718824 name:EI1W3PE02ILQXT TAGTAAAGATGGGGTTTCACTACGTTGGCTGAGCTGTTCTCGAACTCCTGACCTCAAATGAC CTCTGCCTCAGCCTCCCAAGTAGCTGGGATTACAGGCGGGCGCCACCACGCCCAGCTAATTTTGTATTGT TCAGTGAGGGTTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTGTTTTTGAGACAGAATTTTGCTCTT >gnl|ti|1779718825 name:EI1W3PE02GTXK0 TAGTAAAGATGGGGTTTCACTACGTTGGCTGAGCTGTTCTCGAACTCCTGACCTCAAATGAC GTCGCCCAGGCTGGTGTGCAGTGGTGCAACCTCAGCTCACTGCAACCTCTGCCTCCAGGTTCAAGCAATT TCAGAATACCTGTTGCCCATTTTTATATGTTCCTTGGAGAAATGTCAATTCAGAGCTTTTGCTCAGCTTT >gnl|ti|1779718825 name:EI1W3PE02GTXK0 CTCTGCCTCAGCCTCCCAAGTAGCTGGGATTACAGGCGGGCGCCACCACGCCCAGCTAATTTTGTATTGT TAATATGTTTATTTGTTTTGCTGCTGTTGAGTTGTACAATGTTGGGGAAAACAGTCGCACAACACCCGGC TCAGAATACCTGTTGCCCATTTTTATATGTTCCTTGGAGAAATGTCAATTCAGAGCTTTTGCTCAGCTTT TAGTAAAGATGGGGTTTCACTACGTTGGCTGAGCTGTTCTCGAACTCCTGACCTCAAATGAC AGGTACTTTGAGTCTGGGGGAGACAAAGGAGTTAGAAAGAGAGAGAATAAGCACTTAAAAGGCGGGTCCA TAATATGTTTATTTGTTTTGCTGCTGTTGAGTTGTACAATGTTGGGGAAAACAGTCGCACAACACCCGGC >gnl|ti|1779718825 name:EI1W3PE02GTXK0 GGGGGCCCGAGCATCGGAGGGTTGCTCATGGCCCACAGTTGTCAGGCTCCACCTAATTAAATGGTTTACA AGGTACTTTGAGTCTGGGGGAGACAAAGGAGTTAGAAAGAGAGAGAATAAGCACTTAAAAGGCGGGTCCA TCAGAATACCTGTTGCCCATTTTTATATGTTCCTTGGAGAAATGTCAATTCAGAGCTTTTGCTCAGCTTT GGGGGCCCGAGCATCGGAGGGTTGCTCATGGCCCACAGTTGTCAGGCTCCACCTAATTAAATGGTTTACA TAATATGTTTATTTGTTTTGCTGCTGTTGAGTTGTACAATGTTGGGGAAAACAGTCGCACAACACCCGGC AGGTACTTTGAGTCTGGGGGAGACAAAGGAGTTAGAAAGAGAGAGAATAAGCACTTAAAAGGCGGGTCCA GGGGGCCCGAGCATCGGAGGGTTGCTCATGGCCCACAGTTGTCAGGCTCCACCTAATTAAATGGTTTACA >gi|88943037|ref|NT_113796.1|Hs1_111515 Homo sapiens chromosome 1 genomic contig, reference assembly GAATTCTGTGAAAGCCTGTAGCTATAAAAAAATGTTGAGCCATAAATACCATCAGAAATAACAAAGGGAG CTTTGAAGTATTCTGAGACTTGTAGGAAGGTGAAGTAAATATCTAATATAATTGTAACAAGTAGTGCTTG >gi|88943037|ref|NT_113796.1|Hs1_111515 Homo sapiens chromosome 1 genomic contig, reference GATTGTATGTTTTTGATTATTTTTTGTTAGGCTGTGATGGGCTCAAGTAATTGAAATTCCTGATGCAAGT assembly AATACAGATGGATTCAGGAGAGGTACTTCCAGGGGGTCAAGGGGAGAAATACCTGTTGGGGGTCAATGCC GAATTCTGTGAAAGCCTGTAGCTATAAAAAAATGTTGAGCCATAAATACCATCAGAAATAACAAAGGGAG CTCCTAATTCTGGAGTAGGGGCTAGGCTAGAATGGTAGAATGCTCAAAAGAATCCAGCGAAGAGGAATAT >gi|88943037|ref|NT_113796.1|Hs1_111515 Homo sapiens chromosome 1 genomic contig, reference CTTTGAAGTATTCTGAGACTTGTAGGAAGGTGAAGTAAATATCTAATATAATTGTAACAAGTAGTGCTTG TTCTGAGATAATAAATAGGACTGTCCCATATTGGAGGCCTTTTTGAACAGTTGTTGTATGGTGACCCTGA assembly GATTGTATGTTTTTGATTATTTTTTGTTAGGCTGTGATGGGCTCAAGTAATTGAAATTCCTGATGCAAGT AATGTACTTTCTCAGATACAGAACACCCTTGGTCAATTGAATACAGATCAATCACTTTAAGTAAGCTAAG GAATTCTGTGAAAGCCTGTAGCTATAAAAAAATGTTGAGCCATAAATACCATCAGAAATAACAAAGGGAG AATACAGATGGATTCAGGAGAGGTACTTCCAGGGGGTCAAGGGGAGAAATACCTGTTGGGGGTCAATGCC TCCTTACTAAATTGATGAGACTTAAACCCATGAAAACTTAACAGCTAAACTCCCTAGTCAACTGGTTTGA CTTTGAAGTATTCTGAGACTTGTAGGAAGGTGAAGTAAATATCTAATATAATTGTAACAAGTAGTGCTTG CTCCTAATTCTGGAGTAGGGGCTAGGCTAGAATGGTAGAATGCTCAAAAGAATCCAGCGAAGAGGAATAT ATCTACTTCTCCAGCAGCTGGGGGAAAAAAGGTGAGAGAAGCAGGATTGAAGCTGCTTCTTTGAATTTAC GATTGTATGTTTTTGATTATTTTTTGTTAGGCTGTGATGGGCTCAAGTAATTGAAATTCCTGATGCAAGT TTCTGAGATAATAAATAGGACTGTCCCATATTGGAGGCCTTTTTGAACAGTTGTTGTATGGTGACCCTGA Mapped reads Quality scores AATACAGATGGATTCAGGAGAGGTACTTCCAGGGGGTCAAGGGGAGAAATACCTGTTGGGGGTCAATGCC AATGTACTTTCTCAGATACAGAACACCCTTGGTCAATTGAATACAGATCAATCACTTTAAGTAAGCTAAG CTCCTAATTCTGGAGTAGGGGCTAGGCTAGAATGGTAGAATGCTCAAAAGAATCCAGCGAAGAGGAATAT TCCTTACTAAATTGATGAGACTTAAACCCATGAAAACTTAACAGCTAAACTCCCTAGTCAACTGGTTTGA TTCTGAGATAATAAATAGGACTGTCCCATATTGGAGGCCTTTTTGAACAGTTGTTGTATGGTGACCCTGA ATCTACTTCTCCAGCAGCTGGGGGAAAAAAGGTGAGAGAAGCAGGATTGAAGCTGCTTCTTTGAATTTAC AATGTACTTTCTCAGATACAGAACACCCTTGGTCAATTGAATACAGATCAATCACTTTAAGTAAGCTAAG TCCTTACTAAATTGATGAGACTTAAACCCATGAAAACTTAACAGCTAAACTCCCTAGTCAACTGGTTTGA ATCTACTTCTCCAGCAGCTGGGGGAAAAAAGGTGAGAGAAGCAGGATTGAAGCTGCTTCTTTGAATTTAC >gnl|ti|1779718824 name:EI1W3PE02ILQXT 28 28 28 28 26 28 28 40 34 14 44 36 23 13 2 27 42 35 21 7 27 42 35 21 6 28name:EI1W3PE02ILQXT 43 36 22 10 27 42 35 20 6 28 43 36 22 9 >gnl|ti|1779718824 2840 4434 3614 2444 1436423 2813 28 28 2742 2835 2621 26735 26 28 2828 2843 2836 2622 28928 2 27 >gnl|ti|1779718824 3 name:EI1W3PE02ILQXT 28 27 4240 3534 21186 28 4328 3628 2227 1033 2724 4226 3528 2028628 28 40 43 33 36 14 22 28 9 36 27 28 2828 2843 28 26 28 22827 35 21 26 26 37928 29 28 28 28 28 27 2813 2828 3742 28 27 27735 2826 36 28 37 36 22 2840 4434 3614 2444 1436 423 28 28 27 28 26 26 27 4240 3534 21 6 28 43 36 22 10 2724 4228 3528 2027 2828 4337 3629 2236 9 27 28 28 28 27 28 28 18 3 28 28 28 27 28 33 24 26 28 28628 28 40 33 14 28 36 27 27 28 27 28 4326 3626 22 9 28 44 36 24 14 4 28 28 28 27 28 26 26 35 26 28 33 23 28 33 23 28 36 27 33 23 28 35 25 28 28 36 37 29 28 28 28 28 27 28 28 28 37 28 27 27 28 36 2827 3736 27 40 3428 1828328 28 28 28 27 33 24 26 28 28 28 40 33 14 28 36 27 28 28 24 28 37 29 28 19 28 26 37 29 26 39 33 13 37 28 27 28 28 28 24 28 28 27 28 28 37 29 36 27 27 28 2728 28 26 2628 3733 29 28 28 28 28 27 28 28 28 37 28 27 27 28 36 28 37 28 21 24 28 27 41 34 15 28 36 27 26 28 24 35 27 28 23 28 33 23 28 36 27 33 23 28 35 25 28 28 36 27 3640 2734 15 28 2828 2828 2728 2824 2828 2837 2429 2828 2819 2728 2826 2837 3729 2926 3639 2733 2713 2837 2728 28 28 3328 2321 2824 3328 2327 2841 3634 2715 3328 2336 2827 3526 2528 2824 2835 3627 2728 3640 2734 15 28 28 28 24 28 37 29 28 19 28 26 37 29 26 39 33 13 37 28 28 28 21 24 28 27 41 34 15 28 36 27 26 28 24 35 27 28 40 34 15 SNP genotype calls Hapmap genotypes 90 209342 16 F 0 0 2110001?0100210010011002122201210211?1221220212000 90 209342 16 F 0180 F 15 16 21100012010021001001100?100201?10111110111?0212000 2110001?0100210010011002122201210211?1221220212000 15 16 M00 18 F 15 90 20934221120010012001201001120010110101011111011110212000 16 F 021100012010021001001100?100201?10111110111?0212000 0 7M 00 15 M 00 2110001?0100210010011002122201210211?1221220212000 2110001001000200122110001111011100111?121210222000 21120010012001201001120010110101011111011110212000 18 F 15 168 F 0 0 7M 00 21100012010021001001100?100201?10111110111?0212000 011202100120022012211200101101210211122111?0120000 15 M2110001001000200122110001111011100111?121210222000 00 8 F 0 12 0 F 9 10 21120010012001201001120010110101011111011110212000 21100010010002001221100010110111001112121210220000 7 M 0011202100120022012211200101101210211122111?0120000 0 M 00 12 F 99 10 2110001001000200122110001111011100111?121210222000 011?001?012002201221120010?10121021112211110120000 8 F 0 21100010010002001221100010110111001112121210220000 0 11 M 7 8 9 M 0 0 011202100120022012211200101101210211122111?0120000 21100210010002001221100012110111001112121210222000 011?001?012002201221120010?10121021112211110120000 12 F 9 10 11 M 7 8 21100010010002001221100010110111001112121210220000 9 M 021100210010002001221100012110111001112121210222000 0 011?001?012002201221120010?10121021112211110120000 11 M 7 8 21100210010002001221100012110111001112121210222000 … … … … … … … rs12095710 T T 9.988139e-01 rs12127179 C T 9.986735e-01 rs11800791 G G 9.977713e-01 rs11578310 G G 9.980062e-01 rs1287622 G G 8.644588e-01 rs11804808 C C 9.977779e-01 rs17471528 A G 5.236099e-01 rs11804835 C C 9.977759e-01 rs11804836 C C 9.977925e-01 rs1287623 G G 9.646510e-01 rs13374307 G G 9.989084e-01 rs12122008 G G 5.121655e-01 rs17431341 A C 5.290652e-01 rs881635 G G 9.978737e-01 rs9700130 A A 9.989940e-01 rs11121600 A A 6.160199e-01 rs12121542 A A 5.555713e-01 rs11121605 T T 8.387705e-01 rs12563779 G G 9.982776e-01 rs11121607 C G 5.639239e-01 rs11121608 G T 5.452936e-01 rs12029742 G G 9.973527e-01 rs562118 C C 9.738776e-01 rs12133533 A C 9.956655e-01 rs11121648 G G 9.077355e-01 rs9662691 C C 9.988648e-01 rs11805141 C C 9.928786e-01 rs1287635 C C 6.113270e-01 Genotype Calling Accuracy vs. Coverage Watson/454 reads NA18507/Illumina reads 100 90 90 80 80 70 70 % Concordance % Concordance 100 60 50 40 60 50 40 30 Binomial (Homo) HMM-Posterior (Homo) 30 20 Binomial (Het) HMM-Posterior (Het) 20 Binomial (Homo) HMM-Posterior (Homo) Binomial (Het) HMM-Posterior (Het) 10 10 0 0 0 1 2 3 4 Avg. Coverage 5 6 0 1 2 3 Avg. Coverage 4 5 6 Conclusions & Ongoing Work Exploiting LD information yields significant improvements in genotyping calling accuracy and/or cost reduction Accuracy achieved by previously proposed binomial test is achieved by HMM-based posterior decoding algorithm using less than 1/4 of the reads Ongoing work Modeling ambiguities in read mapping Haplotype inferrence Extension to population sequencing data (removing need for reference panels) ACKNOWLEDGEMENTS This work was supported in part by NSF under awards IIS-0546457 and DBI-0543365 to IM and IIS-0803440 to YW. SD and YH performed this research as part of the Summer REU program “Bio-Grid Initiatives for Interdisciplinary Research and Education" funded by NSF under award CCF-0755373.