RFfiller: a robust and fast statistical algorithm for gap filling in draft genomes

Gap detection

The algorithm iterates across the scaffold in this stage to look for gaps or an unknown sequence (i.e., the letter N). When a gap is discovered, the surrounding sequences are retrieved based on an overlap threshold. There are three scenarios involved in the extraction, as illustrated in Fig. 2. The first scenario determines if there are enough neighboring sequences to extract from both ends of a gap, which results in an equal number of overlapping sequences surrounding the gap. The second scenario determines whether there are sufficient sequences on the gap’s left side, which results in the extraction of as many sequences as possible from the left side. The third scenario is identical to the second but for the right side of the gap. The second and third scenarios also evaluate whether both ends of the gap contain an adequate sequence for extraction. The remaining sequence from both sides is retrieved in this case.

When a gap is discovered and its neighboring sequence is extracted, it is stored in a file with the ID of the sequence in which the gap was discovered.

Long-read alignment on the gap region

The gap files that are stored are the focus of this stage. The gap files goes through each saved file one by one, searching for a matching sequence by using a long-read alignment on the gap region.

Long-read alignment can be accomplished using a variety of bioinformatics methods. The Basic Local Alignment Search Tool (BLAST) (McGinnis & Madden, 2004) is a well-known aligner. BLAST searches for sequences that share local similarities. Comparing nucleotide or protein sequences to databases and calculating the statistical significance of matches are the primary functions of this program. With the use of BLAST, it is possible to infer functional and evolutionary links between sequences, as well as to identify members of gene families.

BLAST uses a heuristic method to find short matches between two sequences, so it does not consider the entire sequence space when doing its work. The long-read and the gap regions are input sequences. To find sequence information quickly, we create an index of the long-read. Then, we create a BLAST database to search against the long-read sequences. The long read is compared with the gap region, and any matches found are recorded. BLAST tries to attempts to begin local alignments from the initial matches after the first match.

To reduce the number of similar reads, we limit the length of the gap region; if all of the gap regions were allowed to align with the long read, the likelihood of having thousands of similar matches for a specific location would be high. We also set a limit on the percentage of aligned sequences that are identical. Percentage identity is a number that expresses how similar the query and target sequences are (how many characters in each sequence are indistinguishable). The more significant the match, the higher the percentage identity. Based on this premise, we set the percentage identity to 97% to fill the gap with a proper sequence, and any sequence having a match of 98 percent or greater is considered a high-quality match. This later helps the Markov chain to infer accurate reads. The more accurate the reads generated by the blast, the more accurately the Markov chain infers.

The length and location of the match are also restricted. Even if a short-length match were found far away from the gap, it would not be considered a match, regardless of its percentage identity. Due to its alignment length and location, such a sequence does not provide valuable information about the gap; if the aligned sequence were similar to the gap and gap region, it would have longer matches that span the entire gap region or at least one end. Any aligned sequence that fails to meet this constraint is removed from the alignment. The alignment statistics generated by BLAST are stored in a separate directory for use after the alignment process is completed.

Read extension on aligned sequences

At this stage, we go through each gap’s stored filtered alignment files one by one. Our confidence in the sequence alignment allows us to filter the alignment so that it only contains the information we require (i.e., subject and query id, subject and query start, and subject and query end).

We hypothesize that if there is an aligned sequence with sufficient information that meets the constraint set at stage two, then by extending that sequence within the gap’s range, we have a good chance of filling the gap. We then perform the left and right extension procedures for all match sequences based on this hypothesis, as illustrated in Fig. 3.

Before the extension procedure, we perform a sequence directionality check on the alignment statistics file. Sequence directionality (Lodish et al., 2008) is the end-to-end chemical orientation of a single strand of nucleic acid. In a single strand of DNA, the chemistry convention of naming carbon atoms in the nucleotide pentose-sugar-ring means that there will be a 5′-end, which frequently contains a phosphate group that is attached to the 5′  carbon of the ribose ring, and a 3′-end, which is typically unmodified from the ribose -OH substituent. In a DNA double helix, the strands run in opposite directions to permit base pairing between them. A gene is read in the 3′  to 5′  direction (forward sequence), and so its REVERSE complementary strand is read in the 5′  to 3′  direction (reverse sequence). Based on sequence directionality, we split the filtered alignment into the forward and reverse sequence directions. We do this by checking if the query start is larger than the query end. If the query start is greater, the directionality of the sequence is reversed.

After we determine the sequence directionality, we pass it on to the left and right extension procedures.

Left extension

This procedure is only applicable to sequences that are aligned with the gap’s right end (right flank). It extends the aligned sequence from the aligned sequence’s start index to the gap’s start index.

This procedure first determines whether the given sequence is moving in the correct direction. If the sequence is moving in the forward direction, subject start—(gap length + (query start - gap end))—is used to create a new start index. RFfiller then checks if the newly generated start index is contained within the subject sequence’s length. We must check if the extension index is within the sequence because we are extracting the gap-filling sequence from the subject sequence (long-read alignment). If the new start index is not within the range of the subject length, it means that the extended read resides at the beginning of the subject length. If this is the case, we set the new start index to zero. After confirming that the index can be used, the new extension end index is created by simply adding the gap length to the newly created start index.

If the sequence directionality is reversed, we check if the newly generated end index is less than the sequence length. Because the query start is greater than the query end, the extension of the gap length usually extends the limit of the subject sequence. In this case, we set the new end index to the subject length, and the new start index becomes the difference between the new start index and the gap length. Now that we have the final left extended sequence index, we write the subject id, new start index, new end index, and its directionality to a file.

Construction of a transition table

At this stage, we access the stored BED files. The BED files contain candidate sequence indexes to fill all of the scaffold’s gaps. We use BEDTools, a specialized tool for parsing BED files, to convert the index to a sequence. We pass the BED files and the long-read sequence as input parameters. The BEDTools look for indexes in the long-read sequence and then check the sequence’s directionality. If the sequence directionality is forward, the tool searches the index on the long-read sequence directly and outputs the found sequence. However, if the directionality is reversed, it applies the reverse complement to the discovered sequence by searching the index.

We now have enough information about the sequences to fill the gap after converting the indexes into reads. RFfiller generates a transition table for each of the candidate sequences for further analysis. The transition table is a table that shows how the transition function works. It takes two inputs (a state and a symbol) and outputs a state (the “next state”). A transition table is represented by columns corresponding to input symbols.

RFfiller builds a transition table by counting the frequency of each nucleotide in the candidate sequence. To calculate the frequency, we must first convert the sequences into numerical representations: A = 0, C = 1, G = 2, and T = 3. We then count the frequency of each nucleotide pair (how many times a nucleotide follows another nucleotide). We perform this operation for all the candidate sequences, as illustrated in Fig. 4.

After constructing the frequency table, we use probability to convert it to a transition table. In a frequency table, the total sum of a row is used as a divider for the columns within the row. To prevent a division-by-zero error, the sum should always be greater than zero. In some cases, however, we allow the sum to be zero; when a candidate sequence is only one nucleotide long. In this case, we set the rest of the rows to zero to maintain the integrity of the transition table dimension. The final transition matrix is indexed with A, C, G, and T to complete this step.

Markov chain

This stage takes the transition matrix, gap length, and last character before the gap stars (last character) as input parameters. RFfiller generates the final sequence based on the transition matrix, as illustrated in Fig. 5. It first converts the last character to uppercase. This is because some scaffolds are of mixed-case nucleotides (lower- and uppercase), and it disrupts the transition matrix because it is indexed on the uppercase of the nucleotides. It then tries to guess, given the last character, only one nucleotide. This is one of the most crucial steps in this stage. Providing the last character base assures the sequence to be generated is following a pattern from the gap region. RFfiller also guarantees that the Markov chain does not start at a random location. This means that there will always be prior knowledge of the gap before sequences that fill it are generated.

The Markov chain generates the same number of nucleotides as the number of gaps present. We repeat this stage 5–10 times for a single gap because it is a statistical procedure. Now that we have multiple sequences, the next step is to calculate the probability of each newly generated sequence and to select the one with the highest probability. The transition matrix is used to calculate the sequence probability. It checks the probability of the first nucleotide to its preceding nucleotide and records the transition cost, given one of the final sequences. There are n-1 transition costs for an n-length sequence. The sequence rank is determined by the average transition cost of a sequence.

Gap filling

The algorithm creates an output file similar to the scaffold sequence now that we have gathered the sequences to fill the gap. If there are gaps in the scaffold, the algorithm searches for the sequence ID within the newly generated sequences before creating the output file. If the sequence ID exists, the newly generated sequence fills the gap. After that, it is saved to the output file. RFfiller does not affect the rest of the sequence, even if a newly created sequence is not found. It only replaces unknown sequences (labeled by the letter ‘N’).

RFfiller algorithm

Long reads are used to generate a sequence based on a statistical outcome, followed by sequence assembly using a Markov chain-based gap-filling technique. By providing sufficient information about the gap region, RFfiller can forecast the next nucleotide given the current state using the Markov chain concept.

The RFfiller algorithm fills gaps using both a scaffold and a long-read sequence from the same organism. When the scaffold contains unknown sequences or gaps, as in this case, the algorithm is efficient. The algorithm attempts to identify gaps by examining the neighborhood sequence in which they occur. Due to the presence of duplicates in the gap region, it is necessary to inspect the gap region itself. In the gap region, a pattern is determined. This pattern serves as a starting point for evaluating the newly discovered sequence.

The long read may generate read fragments that are identical to those found in the gap region, indicating that the long read may contain the same sequence as the gap region. When the algorithm discovers similar reads, it replaces the unknown sequence with read fragments from those candidates using a Markov chain. The output of the Markov chains is used to close the gap. As illustrated in Fig. 6, RFfiller consists of five stages: gap detection, long-read alignment, read extension, transition table construction, Markov chain construction, and gap filling. Each stage is discussed in detail in the Methods section.

Experiment

To assess accuracy, assemblies were compared with finished sequences derived from independent sequencing experiments. Assemblies were created using popular assembly programs with default parameters as a way to compare the naively generated assemblies. Some gap fillers used short reads while others used long reads. We intended to compare how well a gap filler fills a gap given the same input. In total, there were 36 assemblies from six assembly programs and six datasets. Using the output of the assemblers, we performed 144 gap-filling experiments using four gap-filling algorithms. All experiments were performed on the same server (two Intel Xeon E5-2695 v4 processors with a memory limit of 128 GB). All the parameters used to conduct this experiment can be found in the supplementary material. To keep the article concise, we have chosen to present only the key metrics from each gap filler result in the subsequent evaluations.

Evaluation

This section intends to demonstrate which gap-filling algorithm generates the fewest mismatches and misassemblies, as well as the greatest number of NGA50 and total aligned length, on the filled gap sequence. Additionally, gap-filling algorithms are evaluated in terms of different flanks and gap length, including an equal number of flanks, a single-sided flank, a shorter flank, and a gap length of 1 bp. Different flanks and gap lengths in the same dataset may fall into different categories. According to our analysis of the sequences, those sequences having a large number of 1 bp gaps fell into the 1 bp gap-length category. The same was true for the remaining categories.

Equal number of flanks

To validate this case, the Arabidopsis thaliana, Bacillus cereus, and Staphylococcus aureus datasets were used. RFfiller generated the fewest misassemblies and mismatches. Additionally, it generated the greatest number of total aligned lengths compared with the other gap fillers, as illustrated in Fig. 7. These datasets contain an equal number of flanking sequences. RFfiller took advantage of this to amass sufficient overlap information to construct a well-defined transition table. The precise unknown sequences were then generated using the well-defined transition table. On the three datasets used to validate this concept, RFfiller outperformed the other gap-filling algorithms for an equal number of flanks. Table 2 summarizes the output of the gap fillers.

Table 2:

Difference in assembly statistics of all gap fillers and de novo assemblers across all datasets.

Assemblers	Gap filling algorithms	Arabidopsis thaliana					Bacillus cereus					Escherichia coli
		# gaps	Misassembly	Mismatch	TAL	NGA50	# gaps	Misassembly	Mismatch	TAL	NGA50	# gaps	Misassembly	Mismatch	TAL	NGA50
A5-MiSeq	RFfiller	0	10	1009	−17	–	0	0	20	−160	–	0	0	4	−7	0
	GapCloser	63610	−24	−24172	229264	–	799	0	−3	−108	–	558	0	−127	0	0
	Sealer	220	2	332	12974	–	8	0	−6	−141	–	1	0	1	22	0
	TGS-GapCloser	0	0	0	0	–	0	0	0	0	–	0	0	0	0	0
ABySS	RFfiller	0	0	0	0	–	303	0	19	−516	–	288	2	26	−578	328
	GapCloser	0	0	0	0	–	855	0	−16	−409	–	952	−2	−10	−372	328
	Sealer	0	0	0	0	–	789	0	−6	−194	–	718	−1	−19	−556	328
	TGS-GapCloser	0	0	0	0	–	100	0	0	0	–	202	0	−48	−1924	0
FLYE	RFfiller	0	0	0	0	0	0	5	0	0	–	0	0	0	0	0
	GapCloser	0	0	0	0	0	0	0	0	0	–	0	0	0	0	0
	Sealer	0	0	0	0	0	0	5	0	0	–	0	0	0	0	0
	TGS-GapCloser	0	0	0	0	0	0	0	0	0	–	0	0	0	0	0
SGA	RFfiller	0	0	0	0	–	50	0	8	−1060	–	0	0	0	0	1464
	GapCloser	0	0	0	0	–	864	0	−2	61	–	0	0	0	0	1464
	Sealer	0	0	0	0	–	0	0	0	0	–	0	0	0	0	1464
	TGS-GapCloser	0	0	0	0	–	175	0	−55	−1032	–	0	0	0	0	1464
SOAPdenovo2	RFfiller	763	2	−64	−188	–	786	0	0	−163	–	1343	1	14	−11545	−51
	GapCloser	2938	3	−16	−1637	–	5030	0	−16	−713	–	4103	−2	−234	−6465	40
	Sealer	2880	4	−24	−4577	–	4857	0	32	−408	–	3807	−1	−415	−11068	−29
	TGS-GapCloser	897	0	−67	−628	–	87	0	0	7	–	393	0	−58	−1355	−29
SPAdes	RFfiller	200	0	199	−805	–	100	0	6	−45	–	200	0	0	−300	4722
	GapCloser	80	0	−39	−258	–	695	0	2	297	–	447	0	−7	−134	4722
	Sealer	0	0	0	0	–	600	0	5	153	–	300	0	−6	6	4722
	TGS-GapCloser	0	0	0	0	–	0	0	0	0	–	200	0	−9	212	4722
		Homo sapiens					Saccharomyces cerevisiae					Staphylococcus aureus
A5-MiSeq	RFfiller	118	21	2007	38	–	−5647	2	−8	−1000	0	6097	0	−9990	0	0
	GapCloser	697827	−562	−347997	1186684	–	24	3	−3885	−90	0	6137	0	−10012	0	0
	Sealer	704	12	589	69560	–	−5965	1	−80	−567	0	6097	0	−9998	265	0
	TGS-GapCloser	187	0	−6	102	–	−5729	1	−42	−859	0	6097	0	−10000	0	0
ABySS	RFfiller	0	0	0	0	–	−5472	4	2096	−16125	−100	4917	−25	−300	−1200	−16603
	GapCloser	0	0	0	0	–	199	−3	−1781	−15215	−1	5068	−25	−322	−130	−16603
	Sealer	0	0	0	0	–	−5790	0	2024	−15692	−81	5075	−25	−308	0	−16603
	TGS-GapCloser	0	0	0	0	–	−5554	−3	2062	−15984	3147	5022	−25	−310	−401	−16603
FLYE	RFfiller	0	0	0	0	–	−300	0	0	0	–	300	0	0	0	0
	GapCloser	0	0	0	0	–	−200	0	13	262	–	300	0	0	0	0
	Sealer	0	0	0	0	–	−300	0	0	0	–	300	0	0	0	0
	TGS-GapCloser	0	0	0	0	–	−300	0	0	0	–	300	0	0	0	0
SGA	RFfiller	0	0	0	0	–	−2100	8	−83	−20018	0	10926	4	104	−3911	−603
	GapCloser	0	0	0	0	–	4737	−13	−511	−11565	0	14173	−3	−63	−3771	−25
	Sealer	0	0	0	0	–	−3803	2	−229	−3706	3297	9268	0	0	0	0
	TGS-GapCloser	0	0	0	0	–	−66	−9	−411	−2887	6978	12401	1	−917	5649	−537
SOAPdenovo2	RFfiller	171396	−2	−17853	123776	–	−22310	−7	−477	−41611	1213	39137	−2	−202	−894	0
	GapCloser	164706	14	−1209	43796	–	2992	−23	−1495	−89720	2098	41044	−10	−228	−10522	178
	Sealer	0	0	0	0	–	−7817	−6	−1320	−75616	−820	40348	−10	−178	−7659	178
	TGS-GapCloser	2521	1	−62	−4041	–	−32770	−10	−1455	−32457	694	38116	0	−422	−7659	−71
SPAdes	RFfiller	1534	6	34	−209	–	−4639	3	−136	−10592	0	7439	0	9	−1069	0
	GapCloser	1153	−12	−2	−1393	–	−658	0	−115	−3607	0	7733	0	−46	−638	0
	Sealer	0	0	0	0	–	−51	2	−66	1739	0	7538	0	−45	−965	0
	TGS-GapCloser	190	6	−753	−49977	–	−439	2	−168	2094	0	7339	0	−17	−379	0

DOI: 10.7717/peerj.14186/table-2

Notes:

Values in italicized bold indicate the best result in each category.

The full results of gap fillers and de novo assembly statistics on each dataset can be found in the Supplementary Material.

Performance

We compared RFfiller with the three gap fillers in terms of time and memory consumption. In comparison with the other gap fillers, RFfiller only utilizes reads that are in the neighboring region of the gap sequences, thus RFfiller outperformed other gap fillers in terms of memory usage. For comparing in terms of time, gap filling is dependent on the gap filler’s objective. For example, if a gap filler’s objective is to fill all gaps within a scaffold at the expense of accuracy, the gap filler will run faster. However, if its objective is accuracy rather than speed, it will make extensive use of memory and time. RFfiller’s objective is to find accurate overlapping sequences to tackle mismatches and misassemblies. It achieves its object at the expense of time. Nonetheless, by extracting sufficient overlapping information from gap regions, a highly accurate read can be achieved within a short period of time (Fig. 8). For this reason, we compared gap fillers based on their ability to fill gaps with fewer mismatches and misassemblies relative to the original sequence. RFfiller was the most resource-efficient algorithm of all the gap fillers. The average result of the comparison is illustrated in Fig. 8.