1 / 46

Lecture 5. Sequence Assembly

Lecture 5. Sequence Assembly. The Chinese University of Hong Kong CSCI3220 Algorithms for Bioinformatics. Lecture outline. The sequence assembly problem Several general approaches Related graph problems Hamiltonian path Eulerian path Sequence assembly by using de Bruijn graphs. Part 1.

judah-weber
Télécharger la présentation

Lecture 5. Sequence Assembly

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lecture 5. Sequence Assembly The Chinese University of Hong Kong CSCI3220 Algorithms for Bioinformatics

  2. Lecture outline • The sequence assembly problem • Several general approaches • Related graph problems • Hamiltonian path • Eulerian path • Sequence assembly by using de Bruijn graphs CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  3. Part 1 The Sequence Assembly Problem

  4. Reference sequence • In the last lecture, we studied the problem of short read alignment • Assumptions behind the short read alignment problem: • There is a reference genome • The reference is similar to, but not exactly the same as, the DNA sequence from which the short reads were generated • Sometimes, no good references are available CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  5. Lack of a reference sequence • Some situations in which a good reference sequence is not available: • Sequencing a new type of bacteria • Sequencing a genomic region previously poorly annotated • In these situations, we need to reconstruct the sequence by assembling the short reads. • This process is called sequence assembly, or “de novo assembly”. CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  6. Sequence assembly • Example revisited: Suppose we have got the following short reads from multiple copies of an unknown sequence s: • ACA, ATA, ATA, ATT, TAG, TAT, TTC • How to get back sequence s from these short reads? • One possible formulation: Shortest superstring • Find the shortest string s’ such that every observed read is a substring of s’ • There is no guarantee that s’ must be equal to the actual sequence s • Very difficult problem (NP-hard) – No known polynomial time algorithm exists CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  7. Greedy approach • One heuristic method to solving shortest superstring: Greedy merge of the two strings with the maximum overlap • For example, ATA may be followed by TAG, since the last two characters of ATA are the same as the first two characters of TAG • Similar to playing a jigsaw puzzle Image credit: slideteam.net CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  8. Greedy approach • Detailed steps: • First remove any input short read that is a sub-sequence of another. For example, • If both ACAT and ACA are in the input, remove ACA. • If there are multiple copies of ACA in the input, remove all but one of them. • For the remaining list of short reads, find the ordered pair (x, y) with the maximum overlap between x’s suffix and y’s prefix. If y contains l characters and the overlap involves k characters, remove y and replace x with the merged sequence xy[k+1..l]. • Tie-breaking rule for our discussion: If there are multiple of them, merge the first pair according to lexicographic order. • (1, 2) < (1, 3) < (2, 1) < (2, 3) < (3, 1) < (3, 2) • Repeat steps 1 and 2 until there is only one sequence left. CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  9. Greedy approach example • Input short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • (Removing reads that are substrings of others)ACA, ATA, ATT, TAG, TAT, TTC • ACA, ATA, ATT, TAG, TAT, TTC • ACA, ATAG, ATT, TAT, TTC • ACA, ATAG, ATT, TAT, TTC • ACA, ATAG, ATTC, TAT • ACA, ATAG, ATTC, TAT • ACA, ATTC, TATAG • ACA, ATTC, TATAG • ACATTC, TATAG • ACATTCTATAG • Results: • Reconstructed string different from actual one (TATACATTAG) • Also 1 character more (11 vs. 10) CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  10. How good is the greedy algorithm? • How long is the reconstructed sequence? • It has been proved that if the shortest superstring has length m, then the reconstructed string has length  4m. (Note that m is different from n, the length of the actual sequence s.) • Is there an example with a reconstructed string with length  2m? • It is possible to have m < n, i.e., the shortest superstring may not be the actual sequence s • Is there an example with a reconstructed string with length << n? CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  11. Practical issues • There are many practical issues in sequence assembly: • Non-uniqueness due to short read length • Non-uniqueness due to heterogeneity • Incomplete coverage • Sequencing errors • Ambiguity due to repeats ... CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  12. Practical difficulties and solutions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • Problems: • Non-uniqueness due to short read length: TAT may also be followed by ATT(instead of ATA as in TATACATTAG) • Solution: • Use longer reads and/or paired reads, and assemble reads only when they have substantial overlaps • Limited by current technology: <200nt per read CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  13. Paired end sequencing • Sequence both ends of a fragment • The two reads are called a mate pair • Insert sizes form a distribution due to random fragmentation and manual size selection • One read is likely within a certain distance range from the other • If the location of one read is ambiguous, may use the location information of the other to help • More difficult in practice due to imprecise insert size • The idea of paired end sequencing is also useful in short read alignment Fragment length Read length Insert size Sequencing Read1 Read2 TAT ATA ATT ...CAT ...ATT ...GGG (Suppose s=TATACATTAGGG here) CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  14. Practical difficulties and solutions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • Problems: • Non-uniqueness due to heterogeneity: It is possible that the DNA sample contains two versions of the sequence, one with TATACATTAG and one with TATACATTCG • Solution: • Use an assembly method that can consider multiple possible ways of assembling the reads CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  15. Practical difficulties and solutions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • Problems: • Incomplete coverage: We do not know what is after ACA, since no reads start with CA • Solution: • Produce more reads • The ratio between the total length of useful reads and the length of s is called the average read depth. For sequence alignment, it is now common to perform 30x-60x coverage. For sequence assembly, usually we need >100x. • Sometimes we do not know the actual length of s and need to estimate it • Costs more starting materials, reagents, money and computation CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  16. Practical difficulties and solutions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • Problems: • Sequencing errors: ATT seems to be followed by TTC, but actually s does not contain TTC as a subsequence (recall that s=TATACATTAG), due to a sequencing error (TTC should actually be TTA) • Solutions: • Use longer reads/paired reads so that if a read contains an error, it has low probability of being assembled • Use some statistics to estimate error probability CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  17. Practical difficulties and solutions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • Problems: • Ambiguity due to repeats: Should ATA be put before or after TAT? • This problem is due to the occurrence of two TA’s in s=TATACATTAG • Solution: • Use longer reads/paired reads • Still cannot handle the following cases: • Each unit of a repeat is very long • There are many copies of a repeating unit CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  18. Some concepts in sequence assembly • Key terms: • Contig: A partially assembled sequence from some reads • Scaffold: An arrangement of the contigs with specified order and orientations • Usually the final output does not contain a single sequence, but just some scaffolds • Some descriptive statistics of assembly outputs: • Length of longest contig • Average length of contigs • Total length of contigs • N50: Length of the contig such that shorter contigs amount to 50% or more of the total length of all contigs • If the lengths are (in an arbitrary unit) 10, 8, 6, 5, 3, 3, 2, 1, 1, 1, then the N50 value would be 6, since (10+8+6) = 24, which is larger than 50% of the sum (10+8+6+5+3+3+2+1+1+1) = 40 CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  19. General approaches • Some other proposed approaches to sequence assembly: • Greedy extension: Extend the current contig until there are no more or multiple extensions • Short reads: ACA, ATA, ATA, ATT, TAG, TAT, TTC • ATT ATTC • Overlap/layout/consensus: Perform all-against-all read alignments to find overlaps, deduce approximate layout, and refine layout by multiple sequence alignment • The all-against-all part requires tremendous computational power • de Bruijn graph (More details later) CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  20. Part 2 Related Graph problems

  21. Graph formulations • We can use a graph to provide an abstraction of the short read assembly scenario • Formulation 1: • Each node is a short read • There is an edge from node X to node Y if a suffix of X substantially overlaps with a prefix of node Y • E.g., overlapping at least |X|/2 characters • Goal : Find a path that visits every node exactly once • Because you want every short read to appear exactly once in the resulting sequence • This is a Hamiltonian path problem CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  22. Hamiltonian path formulation • Short reads (new example):ACA, ATA, ATT, CAT, TAC, TAG , TAT, TTA • Suppose a node is connected to another one if the length-2 suffix of the former equals the length-2 prefix of the latter: • Unfortunately, even the decision version (whether such a path exists) is very difficult (NP-complete) ACA ATA TTA ATT TAT CAT TAG TAC CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  23. Graph formulations • Formulation 2: • Each edge is a short read • The two nodes that connect an edge are derived from the prefix and suffix of the corresponding read • Different reads can share nodes • Goal: Find a path that visits every edge exactly once • The Eulerian path problem CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  24. Eulerian path formulation • Short reads:ACA, ATA, ATT, CAT, TAC, TAG , TAT, TTA • Suppose each read is decomposed into two nodes, one containing the length-2 prefix of it and the other the length-2 suffix of it. There is an edge pointing from the former to the latter. • Interestingly, it is much easier to find Eulerian paths AC AT ATT ATA ACA TT AG TTA TAC TAT CAT TAG TA CA CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  25. The Eulerian path problem AC AT • Existence of an Eulerian path: • The in-degree of a node is the number of edges going into it • The out-degree of a node is the number of edges going out of it • If a connected directed graph has an Eulerian path, the followings are true (why?): • At most one node has (out-degree – in-degree) = 1 • At most one node has (out-degree – in-degree) = -1 • All other nodes have (out-degree – in-degree) = 0 • Surprisingly, if a connected graph satisfies these conditions, it must have an Eulerian path, i.e., the three form a set of both necessary and sufficient conditions ATT ATA ACA TT AG TTA TAC TAT CAT TAG TA CA CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  26. The Eulerian path problem AC AT • Finding an Eulerian path (Hierholzer’s algorithm): • Start with the node with an extra out-degree (or if not exist, any node) • Follow any unused edges to visit other nodes, until getting stuck • Either back to the starting node and it has no more unused edge, or the node with an extra in-degree (why are there no other possibilities?) • If any visited node has unused edges, repeat the above with this node as the starting node. The path must end at the same node. Join this new path with the old one. ATT ATA ACA TT AG TTA TAC TAT CAT TAG TA CA CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  27. Example AC AT AC AT 2 TT AG TT AG 3 1 TA CA TA CA 4 AC AT AC AT 6 2 4 iii 2 i TT AG TT AG 7 3 3 ii 1 1 5 iv TA CA TA CA 8 4 CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  28. Comparing the two formulations Actual DNA sequence (which should be unknown): TATACATTAG • The Hierholzer’s algorithm runs in linear time • Does it mean we have proved P=NP? • A general Hamiltonian path problem is NP hard, but due to the equivalence to the Eulerian path formulation, Hamiltonian path problems *for short reads* can be efficiently solved Hamiltonian path formulation: Eulerian path formulation: AC AT ATT ATA ACA ACA ATA TT AG TTA ATT TTA TAC TAT CAT TAG TA CA TAT CAT TAG TAC CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  29. Why the big difference? • Why is the Eulerian path problem much easier than the general Hamiltonian path problem? • Straight necessary and sufficient conditions for an Eulerian path to exist • The conditions are simple: They can be checked very efficiently • For the Eulerian path problem, the solution of a sub-problem (involving only a subset of the edges) can always contribute towards the solution of the original problem • Another example of reusing results CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  30. Back to sequence assembly • While the Eulerian path formulation is elegant, we need to deal with the series of issues when applying it to perform sequence assembly • Non-uniqueness • Errors in data • Repeats • Heterogeneity • ... • We now study some practical methods for solving these issues CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  31. Part 3 de Bruijn Graphs

  32. Problems of the graph formulations • So far we have assumed that each node (in the Hamiltonian path formulation) or each edge (in the Eulerian path formulation) is a short read. • It is problematic if “short” reads are not really that short (but 100-200nt) when determining which reads should be connected in the graph: • If the required overlap is too large, a single error/mutation would make two consecutive reads not connected • TATACATTA, ATAGATTAG • If the required overlap is too small, there will be too many connections and it is easy to have non-unique solutions • There can be a tremendous number of nodes/edges, making the graph not fitting into memory CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  33. de Bruijn graph • In the original definition proposed by Dutch mathematician Nicolaas de Bruijn in 1946: • A de Bruijn graph is a graph that contains every k-mer of a certain alphabet as a node, and there is a directed edge from a node to another one if the length-(k-1) suffix of the former is the same as the length-(k-1) prefix of the latter. • Example: k=3, alphabet={0, 1} Image source: Wikipedia CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  34. de Bruijn graph for sequence assembly • To use de Bruijn graph for sequence assembly: • We consider only k-mers that are subsequences of some observed reads • In practice, only a very small fraction of the 4k possible k-mers appear in the reads • Two nodes are connected only if there are reads that contain both at consecutive positions • So, there will not be an edge from ATA to TAG if ATAG is not observed in any read • The number of reads that support each edge is also recorded CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  35. de Bruijn graph for sequence assembly • Example: • Sequencing reads: a:ACGGC, b:CGGCG, c:CGTGA, d:GACGT, e:GCGTG, f:GGCGT, g:GTGAC and h:TGACG • To construct a de Bruijn graph when each edge corresponds to a 3-mer: 1 ACG CGT 2 1 CGG GTG 2 2 GGC TGA 2 2 2 2 GCG GAC CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  36. Patterns of various issues • The various issues we discussed before will form special patterns in a de Bruijn graph: • Non-uniqueness: frayed rope • s=TACCGGACCGC • Observed reads (not necessarily covering all k-mers of s): TACC, ACCG, CCGG, GACC, CCGC • Incomplete coverage: possibly disconnected graph • s=TATACATTAG • Observed reads: TATA, ATAC, ACAT, CATT, ATTA TAC CGG ACC CCG GAC GGA TAT ATA TAC ACA CAT ATT TTA Image credit: cuttingedgedjs.com CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  37. Patterns of various issues • Tandem repeats: cycle • s=TACCGACCGC • Observed reads: ACCG, CCGA, CGAC, GACC, CCGC • Errors at read ends: spur • s=TATACATTAG • Observed reads: TATA, TATT, ATAC, TACA, ACAT CGC ACC CCG CGA GAC TAT ATA TAC ACA CAT ATT Image credit: Wikimedia CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  38. Patterns of various issues • Heterogeneity/errors at read centers: bubble • s=TATACATTAG; s’=TATAGATTAG • Observed reads: TATA, ATAC, TACA, ACAT, CATT, ATAG, TAGA, AGAT, GATT TAC ACA CAT TAT ATA ATT TAG AGA GAT CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  39. Resolving problems • Some methods for resolving the problems (high-level ideas only): • Handling potential errors: • Pre-filtering k-mers supported by few reads • Bimodal distribution of k-mer frequencies: one peak corresponds to legitimate k-mers, the other (much lower) peak due to errors • May also use base quality scores in filtering • Remove paths supported by few reads • Combine near-identical paths CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  40. Resolving problems • Some methods for resolving the problems (high-level ideas only): • Non-uniqueness, heterogeneity: • Duplicate the shared part in frayed ropes and bubbles into separate paths, if supported by read counts • Repeats: • Use read counts to deduce number of copies • Usually not very accurate due to random fluctuations in read counts CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  41. Epilogue Case Study, Summary and Further Readings

  42. Case study: “Synthetic cell” • Topic of this lecture so far: • We have multiple copies of an unknown (biological) DNA sequence • We cut them down into small fragments • We sequence each of them to get short reads (text strings) • We assemble the short reads to get back the original (text) sequence • Is it possible to do the opposite? • We have a long text string s • We cut it down into small strings • We biochemically synthesize each of them • We assemble them to get a DNA molecule with sequence s CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  43. Case study: “Synthetic cell” • In 2008, a team reported such an experiment: • They took the DNA sequence (text string) of a bacterium called Mycoplasma genitalium • Total length: 582,970 base pairs • They synthesized the DNA molecule in a hierarchical manner, with some changes to the sequence • Cassettes of 5-7kb  intermediate assemblies of ~24kb  ~72kb  ~144kb  full sequence CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  44. Case study: “Synthetic cell” • In 2010, the team reported something more: • They synthesized the DNA molecule of a bacterium, Mycoplasma mycoides (1.1Mb) • They then transplanted it into a cell from a closely related species, Mycoplasma capricolum • The cell did not divide • Was it “alive”? • Later they found that there was a frameshift mutation in an important gene • After correcting it, cells receiving the sequence were able to divide • The study stirred up a lot of heated debates Image credit: Gibson et al., Science 329(5987):52-56, (2010) CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  45. Summary • The sequence assembly problem is to assemble the original sequence from short sequencing reads without a reference • Two graph formulations: • Read as nodes: Hamiltonian path problem • Hamiltonian path problems are NP-hard in general • Read as edges: Eulerian path problem • Eulerian path problems can be solved efficiently • Current standard is to use k-mer-based de Bruijn graphs • Complications due to various issues: • Heterogenetity/heterozygosity • Sequence errors • Non-uniqueness of sub-sequence, repeats • Incomplete coverage CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

  46. Further readings • Review papers: • Flicek and Birney, Sense from Sequence Reads: Methods for Alignment and Assembly. Nature Methods 6(11s):S6-S12, (2009) • Miller et al., Assembly Algorithms for Next-Generation Sequencing Data. Genomics 95(6):315-327, (2010) • Compeau et al., How to Apply de Bruijn Graphs to Genome Assembly. Nature Biotechnology 29(11):987-991, (2011) CSCI3220 Algorithms for Bioinformatics | Kevin Yip-cse-cuhk | Fall 2014

More Related