Lesson 3.9
Assembling Genomes from Read-Pairs
Watch the lectureAssembling Read-Pairs
From reads to read-pairs
Previously, we described an idealized form of genome assembly in order to build up your intuition about de Bruijn graphs. In the rest of the chapter, we will discuss a number of practically motivated topics that will help you appreciate the advanced methods used by modern assemblers.
We have already mentioned that assembling reads sampled from a randomly generated text is trivial, since random strings are not expected to have long repeats. Moreover, de Bruijn graphs become less and less tangled when read length increases (see the figure below). As soon as read length exceeds the length of all repeats in a genome (provided the reads have no errors), the de Bruijn graph turns into a path. However, despite many attempts, biologists have not yet figured out how to generate long and accurate reads. The most accurate sequencing technologies available today generate reads that are only about 300 nucleotides long, which is too short to span most repeats, even in short bacterial genomes.
Figure: The graph DeBruijn4(TAATGCCATGGGATGTT) (top right) is less tangled than the graph DeBruijn3(TAATGCCATGGGATGTT) (top left). The graph DeBruijn5(TAATGCCATGGGATGTT) is a path (bottom).
We saw earlier that the string TAATGCCATGGGATGTT cannot be uniquely reconstructed from its 3-mer composition since another string (TAATGGGATGCCATGTT) has the same 3-mer composition.
STOP and Think What additional experimental information would allow you to uniquely reconstruct the string TAATGCCATGGGATGTT?
Increasing read length would help identify the correct assembly, but since increasing read length presents a difficult experimental problem, biologists have devised an ingenious experimental approach to increase read length by generating read-pairs, which are pairs of reads separated by a fixed distance d in the genome (see figure below). You can think about a read-pair as a long “gapped" read of length k + d + k whose first and last k-mers are known but whose middle segment of length d is unknown. Nevertheless, read-pairs contain more information than k-mers alone, and so we should be able to use them to improve our assemblies. If only you could infer the nucleotides in the middle segment of a read-pair, you would immediately increase the read length from k to 2 · k + d.
Figure: Read-pairs sampled from TAATGCCATGGGATGTT and formed by reads of length 3 separated by a gap of length 1. A simple but inefficient way to assemble these read-pairs is to construct the de Bruijn graph of individual reads (3-mers) within the read-pairs.
Transforming read-pairs into long virtual reads
Let Reads be the collection of all 2N k-mer reads taken from N read-pairs. Note that a read-pair formed by k-mer reads Read1 and Read2 corresponds to two edges in the de Bruijn graph DeBruijnk(Reads). Since these reads are separated by distance d in the genome, there must be a path of length k + d in DeBruijnk(Reads) connecting the node at the beginning of the edge corresponding to Read1 with the node at the beginning of the edge corresponding to Read2, as shown in the figure below. If there is only one path of length k + d + 1 connecting these nodes, or if all such paths spell out the same string, then we can transform a read-pair formed by reads Read1 and Read2 into a virtual read of length 2 · k + d that starts as Read1, spells out this path, and ends with Read2.
For example, consider the de Bruijn graph below, which is generated from all reads present in the read-pairs from the previous step. There is a unique string spelled by paths of length k + d + 1 = 5 between edges labeled AAT and CCA within a read-pair represented by the gapped read AAT-CCA. Thus, from two short reads of length k, we have generated a long virtual read of length 2 · k + d, achieving computationally what researchers still cannot achieve experimentally! After preprocessing the de Bruijn graph to produce long virtual reads, we can simply construct the de Bruijn graph from these long reads and use it for genome assembly.
Figure: Each path of length 3 + 1 + 1 = 5 (connecting node AAT to node CCA) spells out AATGCCA. Thus, the gapped read AAT-CCA can be transformed into a long virtual read AATGCCA, shown by the highlighted path. (There are three such paths because there are three possible choices of edges labeled ATG.
Although the idea of transforming read-pairs into long virtual reads is used in many assembly programs, we have made an optimistic assumption: “If there is only one path of length k + d + 1 connecting these nodes, or if all such paths spell out the same string" In practice, this assumption limits the application of the long virtual read approach to assembling read-pairs because highly repetitive genomic regions often contain multiple paths of the same length between two edges spelling different strings (see the figure below). If this is the case, then we cannot reliably transform a read-pair into a long read. Instead, we will describe an alternative approach to analyzing read-pairs.
Figure: (Left) The highlighted path in DeBruijn3(AATCTGACATATGG) spells out the long virtual read AATCTGACA, which is a substring of AATCTGACATATGG. (Right) The highlighted path in the same graph spells out the long virtual read AATATGACA, which does not occur in AATCTGACATATGG.
From composition to paired composition
Given a string Text, a (k,d)-mer is a pair of k-mers in Text separated by distance d. We use the notation (Pattern1|Pattern2) to refer to a (k,d)-mer whose k-mers are Pattern1 and Pattern2. For example, (AAT|TGG) is a (3,4)-mer in TAATGCCATGGGATGTT. The (k,d)-mer composition of Text, denoted PairedCompositionk,d(Text), is the collection of all (k,d)- mers in Text (including repeated (k,d)-mers). For example, here is PairedComposition3,1(TAATGCCATGGGATGTT):
TAA GCC
AAT CCA
ATG CAT
TGC
ATG
GCC TGG
CCA GGG
CAT GGA
ATG GAT
TGG ATG
GGG TGT
GGA GTT
TAATGCCATGGGATGTT
Since the order of (3,1)-mers is unknown, we list them according to the lexicographic order of the 6-mers formed by their concatenated 3-mers:
(AAT|CCA) (ATG|CAT) (ATG|GAT) (CAT|GGA) (CCA|GGG) (GCC|TGG) (GGA|GTT) (GGG|TGT) (TAA|GCC) (TGC|ATG) (TGG|ATG)
Exercise Break Generate the (3,2)-mer composition of TAATGCCATGGGATGTT in lexicographic order. Include repeats, and return your answer as a list on a single line. As a hint to help you with formatting, your answer should begin "(AAT|CAT) (ATG|ATG)..."
Note that while there are repeated 3-mers in the 3-mer composition of this string, there are no repeated (3,1)-mers in its paired composition. Furthermore, although TAATGCCATGGGATGTT and TAATGGGATGCCATGTT have the same 3-mer composition, they have different (3,1)-mer compositions. Thus, if we can generate the (3,1)-mer composition of these strings, then we will be able to distinguish between them. But how can we reconstruct a string from its (k,d)-mer composition? And can we adapt the de Bruijn graph approach for this purpose?
String Reconstruction from Read-Pairs Problem: Reconstruct a string from its paired composition.
- Input: A collection of paired k-mers PairedReads and an integer d.
- Output: A string Text with (k,d)-mer composition equal to PairedReads (if such a string exists).
Paired de Bruijn graphs
Given a (k, d)-mer (a1 ... ak | b1 ... bk), we define its prefix as the (k − 1, d + 1)-mer (a1 ... ak -1 | b1 ... bk - 1), and its suffix as the (k - 1,d + 1)-mer (a2 ... ak | b2 ... bk). For example, Prefix((GAC|TCA)) = (GA|TC) and Suffix((GAC|TCA)) = (AC|CA).
Note that for consecutive (k, d)-mers appearing in Text, the suffix of the first (k, d)-mer is equal to the prefix of the second (k, d)-mer. For example, for the consecutive (k, d)-mers (TAA|GCC) and (AAT | CCA) in TAATGCCATGGGATGTT, Suffix((TAA | GCC)) = Prefix((AAT | CCA))=(AA | CC).
Given a string Text, we construct a graph PathGraphk,d(Text) that represents a path formed by |Text| - ( k + d + k) + 1 edges corresponding to all (k, d)-mers in Text. We label edges in this path by (k, d)-mers and label the starting and ending nodes of an edge by its prefix and suffix, respectively. The figure below illustrates PathGraph3,1(TAATGCCATGGGATGTT).
Figure: PathGraph3,1( TAATGCCATGGGATGTT). Each (3,1)-mer has been displayed as a two-line expression to save space.
The paired de Bruijn graph, denoted DeBruijnk,d(Text), is formed by gluing identically labeled nodes in PathGraphk,d(Text); the figure below shows the construction of the paired de Bruijn graph DeBruijn3,1(TAATGCCATGGGATGTT). Note that the paired de Bruijn graph is less tangled than the de Bruijn graph constructed from individual reads.
Figure: (Top) PathGraph3,1(TAATGCCATGGGATGTT) is formed by eleven edges and twelve nodes. Only two of these nodes have the same label (TG|AT). (Middle) Bringing the two identically labeled nodes closer to each other in preparation for gluing. (Bottom) The paired de Bruijn graph DeBruijn3,1(TAATGCCATGGGATGTT) is obtained from PathGraph3,1(TAATGCCATGGGATGTT) by gluing the nodes sharing the label (TG|AT). This paired de Bruijn graph has a unique Eulerian path, which spells out TAATGCCATGGGATGTT.
STOP and Think It is easy to construct a paired de Bruijn graph from a string Text. But how can we construct the paired de Bruijn graph from the (k,d)-mer composition of Text?
We define PairedCompositionGraphk,d(Text) as the graph consisting of |Text| - (k + d + k) + 1 isolated edges that are labeled by the (k,d)-mers in Text, and whose nodes are labeled by the prefixes and suffixes of these labels. The figure below illustrates the graph PairedCompositionGraph3,1(TAATGCCATGGGATGTT). As you may have already guessed, gluing identically labeled nodes in PairedCompositionGraphk,d(Text) results in exactly the same de Bruijn graph as gluing identically labeled nodes in PathGraphk,d(Text). In practice, we don’t know Text, but we can form PairedCompositionGraphk,d(Text) directly from the (k,d)-mer composition of Text, and the gluing step will result in the paired de Bruijn graph of this composition. The genome can be reconstructed by following an Eulerian path in this de Bruijn graph.
Figure: The graph PairedCompositionGraph3,1(TAATGCCATGGGATGTT) is a collection of isolated edges. Each edge is labeled by a (3,1)-mer in TAATGCCATGGGATGTT; the starting node of an edge is labeled by the prefix of the edge's (3,1)-mer, and the ending node of an edge is labeled by the suffix of this (3,1)-mer. Gluing identically labeled nodes yields the paired de Bruijn graph from the previous step.
A pitfall of paired de Bruijn graphs
We saw earlier that every solution of the String Reconstruction Problem corresponds to an Eulerian path in the de Bruijn graph constructed from a k-mer com- position. Likewise, every solution of the String Reconstruction from Read-Pairs Problem corresponds to an Eulerian path in the Paired de Bruijn graph constructed from a (k,d)-mer composition.
Exercise Break In the paired de Bruijn graph shown below, reconstruct the genome spelled by the following Eulerian path of (2,1)-mers: (AG|AG) →(GC|GC) → (CA|CT) → (AG|TG) → (GC|GC) → (CT|CT) → (TG|TG) → (GC|GC) → (CT|CA).
Figure: A paired de Bruijn graph constructed from the collection of nine (2,1)-mers (AG|AG), (AG|TG), (CA|CT), (CT|CA), (CT|CT), (GC|GC), (GC|GC), (GC|GC), and (TG|TG).
We also saw that every Eulerian path in the de Bruijn graph constructed from a k-mer composition spells out a solution of the String Reconstruction Problem. But is this the case for the paired de Bruijn graph?
STOP and Think The graph shown below has another Eulerian path: (AG|AG) →(GC|GC) → (CT|CT) → (TG|TG) → (GC|GC) → (CA|CT) → (AG|TG) → (GC|GC) → (CT|CA). Can you reconstruct a genome spelled by this path?
The Exercise Break presented the following Eulerian path formed of nine edges in our graph.
AG-AG → GC-GC → CA-CT → AG-TG → GC-GC → CT-CT → TG-TG → GC-GC → CT-CA
We can arrange the (2,1)-mers in this path into the nine rows shown below, revealing the string AGCAGCTGCTGCA spelled by this path:
AG-AG
GC-GC
CA-CT
AG-TG
GC-GC
CT-CT
TG-TG
GC-GC
CT-CA
AGCAGCTGCTGCA
Now, look again at the same graph, and consider the Eulerian path presented by the STOP and Think question:
AG-AG → GC-GC → CT-CT → TG-TG → GC-GC → CA-CT → AG-TG → GC-GC → CT-CA
An attempt to assemble these (2,1)-mers reveals that not every column has the same nucleotide (see the two columns shown in red below). This example illustrates that not not every Eulerian path in the paired de Bruijn graph constructed from (k,d)-mer composition spells out a solution of the String Reconstruction from Read-Pairs Problem.
AG-AG
GC-GC
CT-CT
TG-TG
GC-GC
CA-CT
AG-TG
GC-GC
CT-CA
AGC?GC?GCTGCA
You are now ready to reconstruct a string from read-pairs and become a genome assembly expert.
Code Challenge Solve the String Reconstruction from Read-Pairs Problem. (Solve at Cogniterra or Rosalind.)Solve it with the autograderFree to start in the interactive text
- Input: Integers k and d followed by a collection of paired k-mers PairedReads.
- Output: A string Text with (k, d)-mer composition equal to PairedReads.
To solve the String Reconstruction from Read-Pairs Problem, you will need to reconstruct a string from its path in the paired de Bruijn graph. Check out "Charging Station: Reconstructing a String in the Paired de Bruijn Graph" in the print companion or at Cogniterra to see how this can be done.