Optimal reference sequence selection for genome assembly using minimum description length principle

Reference assisted assembly requires the use of a reference sequence, as a model, to assist in the assembly of the novel genome. The standard method for identifying the best reference sequence for the assembly of a novel genome aims at counting the number of reads that align to the reference sequence, and then choosing the reference sequence which has the highest number of reads aligning to it. This article explores the use of minimum description length (MDL) principle and its two variants, the two-part MDL and Sophisticated MDL, in identifying the optimal reference sequence for genome assembly. The article compares the MDL based proposed scheme with the standard method coming to the conclusion that “counting the number of reads of the novel genome present in the reference sequence” is not a sufficient condition. Therefore, the proposed MDL scheme includes within itself the standard method of “counting the number of reads that align to the reference sequence” and also moves forward towards looking at the model, the reference sequence, as well, in identifying the optimal reference sequence. The proposed MDL based scheme not only becomes the sufficient criterion for identifying the optimal reference sequence for genome assembly but also improves the reference sequence so that it becomes more suitable for the assembly of the novel genome.


Introduction
Rissanen's minimum description length (MDL) is an inference tool that learns regular features in the data by data compression. MDL uses "code-length" as a measure to identify the best model amongst a set of models. The model which compresses the data the most and presents the smallest code-length is considered the best model. MDL principle stems from Occam's razor principle which states that "entities should not be multiplied beyond necessity", http://www.cs.helsinki.fi/group/cosco/ Teaching/Information/2009/lectures/lecture5a.pdf, stated otherwise, the simplest explanation is the best one, [1][2][3][4][5]. Therefore, MDL principle tries to find the simplest explanation (model) to the phenomenon (data).

Methods
The relevance of MDL to Genome assembly can be realized by understanding that Genome assembly is an inference problem where the task at hand is to infer the novel genome from read data obtained from sequencing. Genome assembly is broadly divided into comparative assembly and de-novo assembly. In comparative assembly, all reads are aligned with a closely related reference sequence. The alignment process may allow one or more mismatches between each individual read and the reference sequence depending on the user. The alignment http://bsb.eurasipjournals.com/content/2012/1/18 of all the reads creates a "Layout", beyond which the reference sequence is not used any more. The layout helps in producing a consensus sequence, where each base in the sequence is identified by simple majority amongst the bases at that position or via some probabilistic approach. Therefore, this "Alignment-Layout-Consensus" paradigm is used by genome assemblers to infer the novel genome, [27][28][29][30][31][32][33][34][35].
Comparative assembly, therefore, is an inference problem which requires to identify a model that best describes the data. It begins the process by identifying a model, the "reference sequences", most closely related to the set of reads. It then uses the set of reads to build on this model producing a model which overfits the data, the "novel genome", [27,28,34,[36][37][38][39][40][41]. The task of MDL is to identify the model that best describes the data and within comparative assembly framework the same meaning applies to finding the reference sequences that best describes the set of reads.
MDL presents three variants Two-Part MDL, Sophisticated MDL and MiniMax Regret [1]. The application of these will be briefly discussed in what follows.

Two-part MDL
Also called old-style MDL, the two-part MDL chooses the hypothesis which minimizes the sum of two components: A) The code-length of the hypothesis. B) Code-length of the data given the hypothesis.
The two-part MDL selects the hypothesis which minimizes the sum of the code-length of the hypothesis and code-length of the data given the hypothesis, [1,[42][43][44][45][46][47]. The two-part MDL fits perfectly to the comparative assembly problem. The potential hypothesis which is closely related to the data, in comparative assembly, happens to be the reference sequence whereas the data itself happens to be the read data obtained from the sequencing schemes.

Sophisticated MDL
The two components of the two-part MDL can be further divided into three components: where p θ (X ) denotes the distribution of the Data X according to the model θ. The three part code-length assessment process again can be converted into a two-part code-length assessment by combining steps B and C into a single step B.
A) Encoding the model class: l(M i ), where M i belongs to any Model class. B) Code-length of the Data given the hypothesis class , where X stands for any data set.
Item (B) above, i.e., the 'length of the encoded data given the hypothesis' is also called the "stochastic complexity" of the model. Furthermore, if the data is fixed, or if item (B) is constant, then the job reduces to minimizing l(M i ), otherwise, reducing part (A), [1,[48][49][50][51][52][53].

MiniMax regret
MiniMax Regret relies on the minimization of the worst case regret, [49,50,[53][54][55][56][57][58][59]: where M can be any model, M represents the best model in the class of all models and X denotes the data. The Regret, R M i ,X , is defined as Here the loss function, loss(M i , X ), could be defined as the code-length of the data X , given the model class M i . The application of Sophisticated MDL in the framework of comparative assembly will be discussed in what follows.

Sophisticated MDL and genome assembly
In reference assisted assembly, also known as comparative assembly, a reference sequence is used to assemble a novel genome from a set of reads. Therefore, the best model is the reference sequence most closely related to the novel genome and the data at hand are the set of reads.
However, it should be pointed out that the aim is not to find a general model, rather, the aim is to find a "model that best overfits the data" since there is just one or maybe two instances of the data, based on how many runs of the experiment took place. One "run" is a technical term specifying that the genome was sequenced once and the data was obtained. The term "model that best overfits the data" can be explained using the following example.
Assume one has three Reads {X, Y, and Z} each having n number of bases. Say reference sequences (L) and (M), where (L) = XXYYZZ and (M) = XYZ contains all three reads placed side by side. Since both models contain all the three reads, the stochastic complexity of both (L) and (M) is the same and both overfit the data perfectly. However, since (M) is shorter than (L), therefore (M) is the model of choice on account of being the model that "best" overfits the data. http://bsb.eurasipjournals.com/content/2012/1/18 The table shows that choosing the reference sequence which has the highest number of reads present is not a sufficient condition. Just by looking at the "Data given the model" ≡ "Number of reads found" one ends up choosing Human Chromosome 21. However, looking at the fact that Chromosome 21 is about 9× larger than S85 one realizes that actually S85 is the model of choice. Furthermore, S85 is a bacterial genome whereas Chromosome 21 comes from a eukaryote genome. PAb1 is also a bacteria, therefore, S85 is most definitely the model of choice.
To formalize the MDL process, the first step would be to identify the following considerations: Data given the hypothesis is typically defined as "Number of reads that align to the Ref. sequence". In the case presented below "data given the hypothesis" is defined in an inverted fashion as the "Number of reads that do not align to the reference sequence". These two are interchangeable as the "Total number of reads" is the sum total of the "number of reads that aligned to the Ref. " and the "number of reads that do not align to the Ref. ". Table 1 shows that choosing the reference sequence having the highest number of reads present is not a sufficient condition for selecting the optimal reference sequence. The simulation carried out compared two reference sequences Fibrobacter succinogenes S85 (NC 013410.1), [60,61], and Human Chromosome 21 (AC 000044.1), [62][63][64], with the reads of Pseudomonas aeruginosa PAb1 (SRX000424), [48,65,66]. It shows that in order to choose the optimal reference sequence one has to take into account both the "Code-length of the model" and "Number of reads found" to be the sufficient conditions for choosing the optimal reference sequence.
Therefore, a simple yet novel scheme is proposed for the solution to the problem, see Figure 1 and Table 2. The proposed scheme follows the three assessment process of Sophisticated MDL. The MDL based proposed scheme stores the model class (Ref. sequence), the parameters of the model (where each base of the sequence is flagged with {−1, 0, 1}) and the data given the hypothesis (reads of the novel genome that do not align to the Ref. sequence) is one file. The file is than encoded using either Huffman Coding [67][68][69][70] or Shannon-Fano coding [68][69][70][71] to determine the code-length. For a simplistic three bits per character coding the code-length is measured according to Equation (3).

Figure 1 MDL proposed scheme:
The output of the system shows that the three components of the encoding scheme are separated from one another by ">". The scheme follows the format "Model > Model given the Data > Data given the hypothesis". In the genome assembly framework the scheme mentioned above translates into "Reference Sequence > Reference Sequence according to the set of reads > Set of reads according to the Reference sequence". "Model given the Data" is identified using {−1, 0, 1}. "1"(s) represent the base locations where the reads are found. "0"(s) represents the locations which are not covered by any read. "−1"(s) represents the locations of the genome that are inverted.
Here the loss function, loss(M i , X ), happens to be code-length of the data X , given the model class M i . Whereas, "Data given the hypothesis", is the code-length of the "Reads that do not align to the reference sequence". The code-length in the last column is measured according to Equation (3). The experiment shows that given the MDL proposed scheme Ref. 1 is the optimal choice for a reference sequence. http://bsb.eurasipjournals.com/content/2012/1/18 The proposed scheme not only allows to determine the best model, amongst the pool of models to choose from, but also improves the model to be better suited according to the novel genome to be assembled. This is done by identifying all insertions and inversions, larger than one read length. It then removes those insertions and rectifies those inversions to get a better model, better suited to assemble the novel genome compared to what was started from, see

MDL algorithm
The pseudo code for analysis using sophisticated MDL and the scheme proposed in Section 2.4 is shown in Algorithm 1.
Given the reference sequence S R and K set of reads, {r 1 , r 2 , . . . , r K } ∈ R, obtained from the FASTQ [72,73] file, the first step in the inference process is to filter all low quality reads. Lines 3-10 filters all the reads that contain the base N in them and also the reads which are of low quality leaving behind a set of O reads to be used for further analysis. This pre-processing step is common to all assemblers. Once all the low quality reads are filtered out, the remaining set of O reads are sorted and then collapsed so that only unique reads remain.
Lines 13-27 describe the implementation of the proposed scheme as defined in Section 2.4. Assume that S R is l bases long, and the length of each read is p. Therefore, φ S R picks up p bases at a time from S R and checks whether or not φ S R is present in the set of collapsed reads R . In the event φ S R ∈ R then the corresponding location on S R , i.e., j → j + p are flagged with "1(s)". If φ S R ∈ R , then invert The sequence generated has both yellow and blue regions rectified. Notice that using a simple ad-hoc scheme of counting the number of reads in the reference sequence one would have made use of (b) for assembly of novel genome. However, using MDL one can now use (c) for the assembly of the novel genome.

Algorithm 1 MDL Analysis of a Ref. sequence given a set of reads of the unassembled genome.
1: Input reference sequence S R ; 2: Input read data set {r 1 , r 2 , . . . , r K } ∈ R; 3: for i : 1 → K do 4: if r i contains base N then 5: remove r i from the set of reads; 6: end if 7: if r i has low quality bases then 8: remove r i from the set of reads; 9: end if 10: end for 11: Sort remaining set of reads {r 1 , r 2 , . . . , r O } ∈ R 12: Collapse duplicated reads. 13: for j : 1 → l do 14: read 15: if φ S R = r k ∈ R then 16: flag 1(s) in locations j → j + p 17: flag read r k to be present. 18: if ψ S R = r q ∈ R then 21: flag −1(s) in locations j → j + p 22: flag read r q to be present 23: else 24: flag 0(s) in locations j → j + p 25: end if 26: end if 27: end for 28: for j : 1 → l do 29: modified sequence ← S R 30: identify all inversions by looking at −1 flags 31: start = start of an inversion 32: end = end of an inversion 33: invert genome start → end 34: end for 35: for j : 1 → l do 36: identify all insertions by looking at 0 flags 37: start = start of an insertion 38: end = end of an insertion 39: if τ 1 < end − start < τ 2 then 40: remove segment of genome start → end φ S R → ψ S R and check whether or not ψ S R ∈ R . If yes, then mark the corresponding location on S R , i.e., j → j + p with "−1(s)" and flag φ S R to be present in R . Otherwise, mark the corresponding locations on S R as "0(s)".
Lines 28-34 generates a modified sequence which has all the inversions rectified in the original sequence S R . Lines 35-44 identifies all insertions larger than τ 1 and smaller than τ 2 and removes them, see Figure 3. Here τ 1 and τ 2 are user-defined. Care should be taken to avoid removing very large insertions as this may affect the overall performance in deciding the best sequence for genome assembly. Lines 45-47 removes all the reads that are present in the original S R and the modified sequence identified by flags 1 and −1. In the end the code-lengths are identified by any popular encoding scheme like Huffman [67][68][69][70] or Shannon-Fano coding [68][69][70][71]. If ξ is the smallest code-length amongst all models then use as a reference for the assembly of the unassembled genome rather than using S R .

Results
Simulations were carried out on both synthetic data as well as real data. At first, the MDL process was analyzed on synthetic data on four different sets of mutations by varying the number and length of {Single nucleotide polymorphisms (SNPs), Inversions, Insertions, and Deletions}. The experiments using synthetic data were carried out by generating a sequence S N . The set of reads were derived from S N and sorted using quick sort algorithm [74,75]. Each experiment modified S N to produce two reference sequences S R1 and S R2 by randomly putting in the four set of mutations. The choice of the best reference sequence was determined by the code-length generated by the MDL process. See Tables 3, 4, 5, and 6 for results.
Once the robustness of MDL scheme on each of the four types of mutations was confirmed two-set of experiments were carried out on real data using Influenza viruses A, B, and C which belong to the Orthomyxoviridae group. Influenza virus A has five different strains, i.e., {H1N1, H5N1, H2N2, H3N2, H9N2}, while Influenza viruses B and C each have just one. The genomes of Influenza viruses is divided into a number of segments. Influenza virus A and B each have eight segments while virus C has seven segments, [76][77][78]. Amongst the first segments of each of the viruses only one was randomly selected and then modified to be our novel genome, S N . Reads were then derived from S N and compared  The location and length of these insertions was chosen randomly. 136 196 shows that out of 196 insertions in S R1 only 136 were removed. The remaining insertions were not recovered due to the choice of τ 1 and τ 2 . S R2 has higher number of insertions as opposed to S R1 . The code-length suggests that S R1 is the model of choice as it has a smaller code-length. http://bsb.eurasipjournals.com/content/2012/1/18 The location and length of these deletions was chosen randomly. S R2 has higher number of deletions as opposed to S R1 . The code-length suggests that S R1 is the model of choice as it has a smaller code-length. The experiment show that although no insertions were put in the actual sequence yet still two and three insertions were found for S R1 and S R2 , respectively. This may be due to a large section of reads that could not align to the reference sequence on the edges of these deletions.
with all the seven reference sequences. See Table 7 for results.
The second-set of experiments analyzed the performance of the MDL proposed scheme on reference sequences of various lengths. The test was designed to check whether the proposed scheme chooses smaller reference sequence with more number of unaligned reads or does it choose the optimal reference sequence for assembly. The reads were derived from Influenza A virus (A Puerto Rico 834 (H1N1)) segment 1. All the reference sequences used in this test were also derived from the same H1N1 virus, however, with different lengths, see Tables 8 and 9.

Discussion
The MDL proposed scheme was tested using two-set of experiments. In the first set the robustness of the proposed scheme was tested using reference sequences, both real and simulated, having four types of mutations {Inversions, Insertions, Deletions, SNPs} compared to the novel genome. This was done with the help of a program called change sequence. The program 'change sequence' requires the user to input ϒ m , the probability of mutation, in addition to the original sequence from which the reference sequences are being derived. It start by traversing along the length of the genome, and each time it arrives at a new base, a uniformly distributed random generator generates a number between 0 and 100. If the number generated is less than or equal to ϒ m a mutation is introduced. Once the decision to introduce a mutation is made, the choice of which mutation still needs to be made. This is done by rolling a biased four sided dice. Where each face of the dice represents a particular mutation, i.e., {inversion, deletion, insertion, and SNPs}. The percentage bias for each face of the dice is provided by the user as four additional inputs, ϒ inv , for the percentage bias for inversions, ϒ indel , representing percentage bias for insertions and deletions and ϒ SNP for SNPs. If Table 6  Both S R1 and S R2 have the same code-length. This is because the MDL scheme not only detected all the inversions for S R2 but also recovered all of them. So effectively S R2 ≡ S R1 after the MDL process as explained in Figure 2. , and is therefore, the model of choice. The experiment also shows that given the optimal reference sequence, in this case Influenza virus B, the MDL process rectifies all inversions (4/4). However, given non-optimal reference sequences, the proposed MDL process is not able to rectify the inversions (0/4). So the proposed algorithm chooses the optimal reference sequence, and given the optimal reference sequence if not all, at least most of the inversions are also corrected. http://bsb.eurasipjournals.com/content/2012/1/18 the dice chooses inversion, insertion or deletion as a possible mutation it still needs to choose the length of the mutation. This requires one last input from the user, ϒ len , identifying the upper threshold limit of the length of the mutation. A uniformly distributed random generator generates a number between 1 and ϒ len , and the number generated corresponds to the length of the mutation. The proposed MDL scheme is shown to work successfully, as it chooses the optimal reference sequence to be the one which has smaller number of SNPs, see Table 3, smaller number of insertions, see Table 4, and smaller number of deletions compared to the novel genome, see Table 5. The proposed MDL scheme is also seen to detect and rectify most, if not all, of the inversions present in the reference sequence, see Table 6. Since the code-length of S R1 is the same as S R2 , and all the inversions of S R2 are rectified, the corrected S R2 sequence and S R1 sequence are equally good for reference assisted assembly.
The experiment carried out using Influenza viruses is shown in Table 7. One sequence was randomly chosen amongst the seven sequences and modified at random locations, using the same 'change sequence' program, to form the novel sequence S N . The novel sequence contained {SNPs = 7, inversions = 4, deletions = 1, insertions = 3} as compared to the original sequence. The MDL process used the reads derived from S N to compare seven sequences and determined Influenza virus B to be optimal reference sequence as it had the smallest code-length. The MDL process rectified all inversions while only one insertion was found. This meant that the remaining two Seq. 125% has a quarter of the actual genome concatenated with the complete H1N1 genome making the total length 125% of H1N1. All other reference sequences were derived in a similar way, see Table 8. The unique set of reads and the reference sequences were tested using the MDL proposed scheme, where the code-length was calculated using Equation (3). The results show that the MDL scheme does not choose smaller reference sequences with more unaligned reads rather it chooses the correct reference sequence, Ref. Seq. 7. It was Ref. Seq. 7 from which all the reads were derived from. Since the MDL scheme chooses Ref. Seq. 7 as the optimal sequence, this experiment further proves the correctness of the reference sequence chosen.
Lastly, the above experiment was repeated using a single set of reads derived from the same H1N1 virus segment 1, but this time containing mutations. The set of reads, 390 in total, were derived using the ART read simulator for NGS with read length 30, standard deviation 10, and mean fragment length of 100, [PUT ART Reference], see Table 9. The results show that the MDL proposed scheme chooses the correct reference sequence, Ref. Seq. 100%, even when all the contending reference sequences are closely related to one another in terms of their genome and length.
All simulations were carried out on Intel Core i5 CPU M430 @ 2.27 GHz, 4 GB RAM. Execution time of MDL proposed scheme have been provided in Table 8.

Conclusions
The article explored the application of Two-Part MDL qualitatively and the application of Sophisticated MDL both qualitatively and quantitatively for selection of the optimal reference sequence for comparatively assembly. The article compared the MDL scheme with the standard method of "counting the number of reads that align to the reference sequence" and found that the standard method is not sufficient for finding the optimal sequence. Therefore, the proposed MDL scheme encompassed within itself the standard method of 'counting the number of reads' by defining it in an inverted fashion as 'counting the number of reads that did not align to the reference sequence' and identified it as the 'data given the hypothesis' . Furthermore, the proposed scheme included the model, i.e., the reference sequence, and identified the parameters (θ M i ) for the model (M i ) by flagging each base of the reference sequence with {−1, 0, 1}. The parameters of the model helped in identifying inversions and thereafter rectifying them. It also identified locations of insertions. Insertions larger than a user defined threshold τ 1 and smaller than τ 2 were removed. Therefore, the proposed MDL scheme not only chooses the optimal reference sequence but also fine-tunes the chosen sequence for a better assembly of the novel genome.
Experiments conducted to test the robustness and correctness of the MDL proposed scheme, both on real and simulated data proved to be successful.