Genome assembly

Note

Genome assembly consists in reconstructing the genome from short sequencing reads. It’s inverse process consists into mapping the reads to a known reference genome, and it’s called read mapping.

There are two major classes of assembly algorithms, OLC and DBG.

Model of a genome assembly

Note

Let’s consider an idealised genome, that represents a long random sequence of four bases and that doesn’t contain repeats or other complex structures. Let’s also consider a simple and error free sequencing strategy. We sample equal length fragments with starting points randomly distributed across the genome.

The “shotgun” process can be compared to a process that samples bases from all genomes positions at random. The chance that any particular base is sampled is very low in a single sampling process. However we perform the sampling process a very large number of times.

The Poisson distribution expresses the probability of a given number of events occurring in a fixed interval of time if these events are independent and identically distributed.

Where refers to the number of reads that overlap a certain genomic position, and refers to the mean sequencing depth.

Let’s consider as the genome size, as the read length, as the number of reads and as the total number of sequenced bases. We know that:

We define a -mer as a subsequence with nucleotides, in general there are -mer subsequence in a sequence of length with . In general the total number of -mers win out whole genome sequencing (WGS) is:

While the coverage depth in term of -mers is:

Therefore the ration between the coverage depth for bases () and that for -mers is then:

It’s possible to estimate the genome size and base coverage as:

We can the estimate the probability that a given base will not be covered:

And therefore, the probability of seeing at leas one read at a given position:

Contigs

Contigs are combinations of overlapping reads that represent contiguous sequence. The result of an assembly is a a set of contigs with gaps.

It’s possible to piece together contigs, with paired-end sequencing we can assemble some contigs into “scaffolds”.

Let’s consider an interval that is as long as a read ( nucleotides), the probability that at least one read starts in is:

Now lets consider a nucleotide at genomic position , this nucleotide is in a gap between contigs if no read starts in the interval:

This interval has length and therefore the probability that no read starts in it is , we can then estimate the number of nucleotides in gaps across the entire assembly as:

Each contig has a unique rightmost read , and the probability that a given read is the rightmost read is the same as the probability that no other read starts within that read. Also, the number of contigs must be equal to the number of rightmost reads, thus the expected number of contigs is:

The expected number of reads per contig is then:

We can also get the expected size of a contig:

Let refer to the minimum portion of that is required to detect an overlap, a group of reads combine to a contig if they are connected by overlaps of length greater or equal to . Given this we can estimate the expected number of contigs, given that we demand an overlap of at least nucleotides:

N50

The N50 is the measure used to estimate the quality of a genome assembly, to find it we arrange the contigs from largest to smallest, we find the position where the contigs cover of the total genome size, and the length of the contig in this position is defined as the N50.

The longer the N50 is, the better the assembly.

We need to find an algorithm that will allow us to take a collection of short NGS sequence reads, and to output a longer string representing the Genome that was sequenced.

De Bruijn graph

Note

Imagine we have a function called that takes a DNA sequence and returns a set of all -mers contained in it, however since we don’t know the original order of the -mers in the genome we show the lexicographically. Let us now put each of the -mers into the node of a graph and connect the graph by edges. To find the sequence based only on a collection of -mers subsequences we search for overlaps between them:

In our graph the solution to our problem was a path that visited every node exactly once, that is an Hamiltonian path.

We can instead label the nodes with the -mer subsequences, and label the edges with these -mers.

We can now merge identically labeled nodes in this graph, whilst retaining the edges. This is the de Bruijn graph of the string, and to find the sequence we need a path that visits every edge exactly once, that is the Eulerian path problem, which fortunately has more efficient algorithms to solve compared to Hamiltonian paths.

Hamiltonian paths

An Hamiltonian path is a path in an undirected or directed graph that visits each vertex exactly once (no need to use all edges).

An Hamiltonian cycle is a Hamiltonian path that is a cycle.

Determining whether Hamiltonian paths and cycles exist in a given graph is NP-complete, as .

Eulerian cycles

An Eulerian cycle is a path that traverses every edge exactly once and returns at the end of the traversal to the start node.

A traversable graph is one that can be drawn without taking a pen from the paper and without retracing the same edge. In such a case the graph is said to have an Eulerian path.

Given that the degree of a vertex is the number of edges starting/ending at the vertex, when the degree of all the vertices is even, the graph is traversable (we have an Eulerian cycle), and we can draw it starting at any vertex (in directed graphs all vertices have ).

If there are exactly two vertices of odd degree and all other vertices are of even degree, there is an Eulerian path (starting at one of the odd vertices), but no cycle (in directed graphs one vertex has , one vertex has , while all others have ).

If there are more than two odd vertices the graph cannot be traversed without repeating an edge.

The minimum read overlap must be bases, such that the corresponding nodes of size can be merged.

To find an Eulerian path in a De Bruijn we use the Hierholzer’s algorithm, with :

1. Make sure the Graph satisfies the degree requirements for an Eulerian path (or cycle) to exist
2. Initialise two stacks, a temporary tpath and a epath for the final solution
3. Choose a suitable starting vertex v (the one with one more outgoing edge for an Eulerian path, or any for a cycle)
4. Push v to tpath
5. Let u=tpath.TOP
6. Check all outgoing edges from u, if they’ve all been visited pop u from tpath and push it to epath, if not select a random outgoing edge , push x to tpath and delete the edge from the set of available edges
7. Repeat from step 5 until tpath is empty

The Eulerian path will be in epath.