A read mapper maps sequencing reads on a reference sequence. Paired-end reads occur in pairs, such that the distance between the reads in a pair (the insert size) satisfies a certain distribution. This additional information can be used to resolve ambiguities in the mapping process, namely when the reads might individually map to several locations in the reference sequence. Treating them as a pair might lead to a unique alignment (in case of doubt, pick the most probable alignment).

Assignment

In this assignment you are expected to write a paired-end read mapper that uses the FM-index constructed for a reference genome. In particular, you are provided with the FM-index of the reference genome of Salmonella enterica subsp. enterica serovar Typhimurium str. 14028S (CP001363¹) and a set of (simulated) paired-end reads². In addition, the reference genome is provided, which can be used to check your results.

Index structure

The FM-index is provided in plain text format and consists of the BWT, a sparse suffix array, a sparse occurrences table and a counts array. In the sparse suffix array only the entries $$SA[k \cdot 32]$$ for $$k \in [0 \dots \left\lfloor{n/32}\right\rfloor]$$ are included. In the sparse occurrences table only the entries $$O[\sigma, k \cdot 128]$$ for $$k \in [0 \dots \left\lfloor{n/128}\right\rfloor]$$ are included. These sparseness factors are introduced to reduce the memory requirements. The other values $$SA[i]$$ and $$O(\sigma, i)$$ can be computed from these sparse tables and the BWT.

Reads

The reads contain no errors and no mutations with respect to the reference genome. All read pairs can be uniquely mapped to the reference genome when handling them as pairs.

Note that DNA is double stranded, and each read pair consists of one read from both strands. The reverse complement of a DNA sequence is formed by reversing the letters, interchanging A and T and interchanging C and G. For example, the reverse complement of ACCTGAG is CTCAGGT. The FM-index is only constructed from the forward direction of the reference genome, so make sure to also handle reads that align to the reverse complement of the reference genome. This can be done by aligning every read in its forward direction, and then repeating the alignment procedure with the reverse complement of the read. The reads in a pair are oriented as follows:

The insert size refers to the distance between the two extreme ends of the alignments of the reads. In the given read data set, insert sizes follow a normal distribution $$N(800, 100)$$. Note that one of the reads in a pair maps to the forward strand in the reference, and the other maps to the reverse complement strand. As shown in the figure, the read that maps to the forward strand, maps before the other read (w.r.t. that strand). It is however not known a priori which of the reads in a pair maps where. In particular, the read in a pair that comes first in the FASTA file is not necessarily the read that maps to the forward strand of the reference.

Output

The output is expected to be in SAM format, as specified in the Sequence Alignment/Map Format Specification³. An example of the output⁴ is provided for the first 5 read pairs. The output must preserve the order of the reads and should be written to standard output.

Specification

Define a class ReadMapper that implements the use of a given FM-index of a reference sequence to map paired-end reads (exact matching) onto the reference sequence. This class must support at least the following methods:

• An initialization method that takes the path prefix of the files that make up an FM-index. These files use extensions as specified in the Files⁵ section. The method also takes two more optional integer arguments that respectively represent the sparseness factor of the occurrence table (default value: 128), and the sparseness factor of the suffix array (default value: 32).

• A method O that takes a symbol $$\sigma$$ and an integer index $$k$$. The method must return $$O[\sigma, k]$$.

• A method LF that takes an integer index $$k$$. The method must return the index $$j$$ such that $$SA[j] = SA[k] - 1$$.

• A method SA that takes an integer index $$k$$. The method must return $$SA[k]$$.

• A method exactMatches that takes a string $$s$$. The method must return a sorted list containing the start positions (0-based) of all exact substring matches of $$s$$ in the reference sequence.

• A method bestPairedMatch that takes two strings $$s_1$$ and $$s_2$$. The method also takes an optional third argument that represents the mean insert size (default value: 800). The method must return the best alignment (exact match) of the given pair of reads $$s_1$$ and $$s_2$$, with respect to the given mean insert size. The method must return this alignment as a tuple containing three elements: i) start of the best alignment of $$s_1$$ in the reference sequence, ii) start of the best alignment of $$s_2$$ in the reference sequence and iii) a Boolean that indicates whether or not the alignment is on the forward strand of the reference sequence.

• A method map that takes two file locations: i) the location of a FASTA file containing a sequence of paired-end reads and ii) the location of a SAM file. The read that form a pair are in consecutive records of the FASTA file. The method also takes an optional third argument that represents the mean insert size (default value: 800). The method must find the best alignment (exact match) of each pair of reads in the given FASTA file, and write the alignment results to a SAM file at the given location (an existing file at that location must be overwritten). The SAM output should follow these rules:

  • the header section should only contain header lines @SQ and @PG

  • the @SQ line should only have SN and LN tags, with the value of the SN tag being the path prefix passed to the initialization method

  • the @PG line should only have an ID tag, whose value is ReadMapper

  • the alignment section should only contain the 11 mandatory fields

  • the field QNAME should contain the identifier of the read taken from the FASTA file

  • the field RNAME should contain the path prefix passed to the initialization method

  • the fields MAPQ, PNEXT, TLEN, SEQ and QUAL should contain an indication that the information is not available or not stored

  • the field RNEXT should contain an equal sign (=), indicating that RNEXT is identical to RNAME

Example

In the following interactive session, we assume all demo.*⁶ and CP001363.*⁷ text files that make up the FM-indices to be located in the current directory, as do the input files demo.reads.fasta⁸ and CP001363.reads.fasta⁹. See the explanation of the files below to understand the content and formatting of these files. In addition, you can also use the FASTA files demo.fasta¹⁰ and CP001363.fasta¹¹ containing the reference genomes, together with the output files demo.sam¹² and CP001363.sam¹³, for debugging purposes.

>>> demo = ReadMapper('demo', 4, 2)
>>> [demo.SA(index) for index in range(11)]
[3, 4, 0, 6, 8, 2, 5, 1, 7, 9, 10]
>>> [demo.LF(index) for index in range(11)]
[5, 0, 10, 6, 8, 7, 1, 2, 3, 4, 9]
>>> [demo.O('A', index) for index in range(11)]
[0, 0, 1, 1, 1, 1, 1, 2, 3, 4, 5]
>>> [demo.O('C', index) for index in range(11)]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
>>> [demo.O('G', index) for index in range(11)]
[0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 2]
>>> [demo.O('T', index) for index in range(11)]
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
>>> [demo.O('$', index) for index in range(11)]
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
>>> demo.exactMatches('GAT')
[5]
>>> demo.exactMatches('AT')
[6, 8]
>>> demo.bestPairedMatch('GCA', 'TA', 8)
(1, 7, True)
>>> demo.bestPairedMatch('TGC', 'TA', 8)
(1, 7, False)
>>> demo.map('demo.reads.fasta14', 'demo.sam15')

>>> CP001363 = ReadMapper('CP001363', 128, 32)
>>> indices = [0, 247380, 685559, 704037, 1502653, 2759112, 3187255, 4870265]
>>> [CP001363.SA(index) for index in indices]
[1952500, 901989, 2630792, 4626188, 3179929, 573015, 6885, 4870265]
>>> [CP001363.LF(index) for index in indices]
[3706925, 2497723, 3847099, 3850517, 2813008, 1870350, 4399036, 1164822]
>>> [CP001363.O('A', index) for index in indices]
[0, 83953, 215950, 220490, 400412, 661685, 746479, 1164822]
>>> [CP001363.O('C', index) for index in indices]
[0, 50134, 167547, 175400, 405799, 705527, 852759, 1271002]
>>> [CP001363.O('G', index) for index in indices]
[0, 61898, 161887, 164554, 377183, 790736, 895905, 1271100]
>>> [CP001363.O('T', index) for index in indices]
[0, 51395, 140174, 143592, 319258, 601163, 692111, 1163340]
>>> [CP001363.O('$', index) for index in indices]
[0, 0, 1, 1, 1, 1, 1, 1]
>>> CP001363.exactMatches('TAAGTAACGATA')
[293380, 4117066, 4213150, 4368115, 4411497]
>>> CP001363.exactMatches('TCAAGCTATA')
[2851405, 3581795, 4362552]
>>> CP001363.bestPairedMatch('TAAGTAACGATA', 'TCAAGCTATA')
(293380, 294198, True)
>>> CP001363.bestPairedMatch('GTGGGTTGCA', 'GCCCGGGGATT')
(291296, 290485, True)
>>> CP001363.bestPairedMatch('ATCGTCAACTCC', 'GGCGGCCGTA')
(2852374, 2851581, True)
>>> CP001363.bestPairedMatch('CTTAACCTTC', 'CAAGACTCAA')
(291320, 290529, True)
>>> CP001363.map('CP001363.reads.fasta16', 'CP001363.sam17')

In the demo example we attempt to align the sequence GAT to the reference sequence AGCAAGATAT. For the sake of the exposition, the sparseness of the suffix array is $$2$$, whereas the sparseness of the occurrence table is $$4$$. The table below contains the arrays $$BWT$$, $$SA$$ and $$O$$, where bold black numbers indicate values that are contained in the sparse index and greyed out numbers represent values that are not indexed.

The table below contains the array $$C$$.

To match $$P = \texttt{GAT}$$, we recursively match its suffixes $$P_i$$. The empty string is a prefix of every suffix, hence $$I(P_3) = [0, 10]$$. From relation (1) we get that $$L(P_2) = L(\texttt{T}P_3) = C(\texttt{T}) + O(\texttt{T}, 0) = 8$$, whereas from relation (2) we obtain that $$R(P_2) = R(\texttt{T}P_3) = C(\texttt{T}) + O(\texttt{T}, 11) - 1 = 9$$. Hence $$I(P_2) = [8, 9]$$. Repeated application of this procedure gives $$I(P_1) = [3, 4]$$ and $$I(P) = I(P_0) = [6, 6]$$. Hence, there is a single match and it starts at position $$SA[6]$$.

To find $$SA[6]$$ we compute $$LF^{j}(6)$$ until we find the minimal $$j$$ such that $$LF^{j}(6) = 0\!\!\!\mod\!\! 2$$: \[ \begin{align} LF(6) &= C(\texttt{A}) + O(\texttt{A}, 6) = 1 \\ LF^{2}(6) &= LF(1) = C(\texttt{A}) + O(\texttt{A}, 1) = 0 = 0\!\!\!\!\!\mod\!\! 2 \end{align} \] Hence, $$SA[6] = SA[0] + 2 = 5$$. We conclude that GAT occurs in AGCAAGATAT starting at position 5.

Files

The indexes and reads are provided in plain text. Here we describe what data each of the input file types contain.

.bwt

The .bwt file contains the BWT string.

.counts

The .counts file contains the counts $$C[\sigma]$$ in lexicographical order, separated from each other by newlines.

.occ

The .occ file contains the sparse occurrence table $$O(\sigma, k)$$, with $$k$$ a multiple of 128. Each line contains the occurrences in lexicographical order, separated by spaces. Elements $$O(\sigma, k)$$ for $$k$$ not a multiple of 128 can be obtained from this table by parsing the relevant subsection of the BWT.

.sa

The .sa file contains the sparse suffix array $$SA[k]$$, for $$k$$ a multiple of 32, separated by newlines. All other elements of the suffix array can be obtained from the sparse suffix array through application of relation (3).

.fasta

The .fasta files contain the reference genome (genome.fasta) and the paired-end reads (pe.fasta). Concerning paired-end reads: read >i/1 is paired with read >i/2.

SAM-output-example.sam

This file contains the output that is expected for the first 5 read pairs.

Formulas

Definition 1

1. The counts table $$C[\sigma]$$ counts for all $$\sigma \in \Sigma$$ how many symbols in $$BWT$$ are lexicographically smaller than $$\sigma$$: \[ \forall \sigma \in \Sigma: C[\sigma] = \#\{k: BWT[k] < \sigma \} \]

2. The occurrence table $$O(\sigma, i)$$ is the number of times $$\sigma$$ occurs in $$S[0..i - 1]$$: \[ O(\sigma, i) = \#\{ k: k < i \text{ and } S[k] = \sigma \} \]

3. Given a pattern $$P$$ we define the interval $$I(P) = [L(P), R(P)]$$ as the interval containing all suffixes of $$S$$ that start with $$P$$: \[ \begin{align} L(P) &= \min \left\{ k: P \text{ is a prefix of } S_{SA[k]} \right\} \tag{1} \\ R(P) &= \max \left\{ k: P \text{ is a prefix of } S_{SA[k]} \right\} \tag{2} \end{align} \]

4. The last-first relation $$LF$$ maps $$i$$ on $$j$$ if and only if $$SA[i] = SA[j] + 1$$.

Lemma 1

The interval $$I(\sigma P)$$ can be computed from $$I(P)$$ and the occurrence tables in the following way: \[ \begin{align} L(\sigma P) &= C(\sigma) + O(\sigma, L(P)) \\ R(\sigma P) &= C(\sigma) + O(\sigma, R(P) + 1) - 1 \end{align} \]

Proof. $$L(\sigma P)$$ refers to a suffix starting with $$\sigma$$, hence it is of the form $$C(\sigma) + k$$ with $$k \in[0, c_\sigma]$$, where $$c_\sigma$$ is the number of $$\sigma$$ in $$S$$. The first $$O(\sigma, L(P))$$ occurrences of $$\sigma$$ in the $$BWT$$ correspond to positions in the suffix array that do not have a prefix $$P$$, else they would be in $$I(P)$$. Hence $$k = O(\sigma, L(P))$$. Analogous for $$R(P)$$, but now we are interested in the number of prefixes starting with $$P$$ or a prefix smaller than $$P$$, hence $$k = O(\sigma, R(P) + 1) - 1$$.

Lemma 2

$$LF(i) = C(BWT[i]) + O(BWT[i], i)$$

Proof. Since the suffixes in a suffix array are sorted, the following follows immediately from the definition of $$LF$$: $$LF$$ maps the $$n$$-th occurrence of $$\sigma$$ in $$BWT$$ on the $$n$$-th occurrence of $$\sigma$$ in the lexicographical sorting of all characters in $$BWT$$. Hence, $$LF$$ can be computed from the counts and occurrences tables.

Lemma 3

Using the $$LF$$ function the following relation between different suffix array entries can be established: \[ SA(i) = SA[LF^{j}(i)] + j \tag{3} \]

Proof. By definition of $$LF$$ this holds for $$j = 1$$. The formula now follows by induction.