Virseqimprover: an integrated pipeline for viral contig error correction, extension, and annotation

Overview

Portions of this text were previously published as part of a preprint (https://doi.org/10.21203/rs.3.rs-3318217/v1). Virseqimprover was formulated for recovering complete viral genomes or viral sequences containing terminal repeats starting from single viral contigs in fasta file format and the metagenomic short reads in fastq file format. However, the steps described here could theoretically be applied to a number of difficult-to-assemble elements, including plasmids or transposable elements. The workflow can be divided into three main steps, the error correction step, the extension step, and the annotation step. The error correction step checks for both the circularity of the contig and the uniformity of the depth of coverage. The error correction and extension steps are done iteratively until the contig cannot be extended anymore. Then the final extended contig is annotated for its protein-coding regions. The output contains the extended contig along with the protein annotation. Figure 1 outlines the three steps and details are described in the following.

Error correction

During the error correction step, the circularity of the contig is checked. Circularity is used as an indicator that the contig recovers the complete genome or it reaches to a terminal repeat, and therefore, if Virseqimprover finds that a contig is circular, it goes to the annotation step, trims the redundant part of the contig, and outputs the contig as the final contig. On the other hand, since repetitive sequences can lead to redundant local assembling process, checking circularity can avoid potential repetitive extension. Figure 2 shows how Virseqimprover checks the circularity of the contig. Assume L_r is the read length and L_s is the length of the contig, Virseqimprover divides the sequence into two parts, G_a and G_b, where G_a starts from (L_s − 2 × L_r) bp to (L_s − L_r) bp and G_b starts from the beginning or from 1 bp to the beginning of G_a or to (L_s − 2 × L_r) bp. G_a is aligned against G_b using BLAST (Ye, McGinnis & Madden, 2006). If any part of G_a aligns with G_b with 95% identity and 95% alignment length, Virseqimprover will try to extend the alignment on both sides of the sequences to get the similar region with the maximum length. Then one of the similar regions is trimmed since having the same region twice is redundant.

After checking the circularity of the contig, Virseqimprover checks the uniformity of the read coverage and uses it as an indicator for chimerism in the contig. First, per base depth of coverage of the contig is calculated using Samtools (Li et al., 2009). For every base position, if its coverage is within 15th to 85th percentile of all the base depth of coverages, it is considered to be within the normal range and the position is marked as normal, otherwise marked as suspicious. Consecutive bases marked as suspicious form suspicious regions and those longer than 1,000 bps are considered to be true suspicious regions. All the regions other than true suspicious regions are flagged as true non-suspicious regions. Virseqimprover chooses the longest true non-suspicious region to extend during the extension step.

Virseqimprover is also applicable to linear viruses via its second extension termination criterion. During iterative assembly, edges are checked for circularity and, if found, trimmed up to 2,000 base pairs and subjected to additional rounds of iterative assembly (maximum 100 iterations). If a region cannot be extended beyond this, and the longest edge is not indicative of circularity, the extension will stop and produce a linear contig, putatively representing an extended linear viral sequence.

Extension

During the extension step, Virseqimprover first extracts the start and end edges of the contigs using BEDTools (Quinlan & Hall, 2010). For each edge region, read length * 1.5 is used as the default edge length. Then, for each contig, all the reads are mapped to the edges of the contig using Salmon (Patro et al., 2017). SPAdes (Bankevich et al., 2012) is used for the local assembly process. The extraction-mapping-assembly step is run iteratively for each contig until it stops growing. When the contig cannot be extended, Virseqimprover trims some bps from both ends of the contig and tries to extend the trimmed contig again. The length of the trimming part ranges from 300 bps to 2,000 bps, depending on how much trimming enables the contig extendable. The logic for trimming the ends is that our empirical investigation shows that assemblers often misassemble in one or both ends of the sequence, causing the assembler to stop prematurely which in turn leads to sequence segmentation. It is observed that trimming some bases from both ends often helps the assembler to continue the assembly in the right direction. After trimming and extending the contig, if the contig gets extended, the new extended contig goes back to the error correction step; if it cannot be extended after trimming, the extension step ends and the contig is moved to the annotation step.

Annotation

During the annotation step, the contig is checked for circularity again. If it is circular, the contig is trimmed to remove the redundant sequence. Then Pilon (Walker et al., 2014) is applied to the contig to improve the assembly by correcting single nucleotide polymorphisms (SNPs), insertions and deletions. The inputs of Pilon include a genome/contig in FASTA format and reads mapped to the contig in BAM format. From the alignment information, Pilon creates a pileup structure and then corrects the base based on the frequency of each nucleotide in a position. During the base correction step, Pilon also considers if the reads are properly paired or not and the mapping quality of the base. If the alignment of read pairs indicates a discrepancy in the assembly, Pilon tries to fix the assembly by doing a local reassembly in those places. The improved contig after Pilon is sent to Prodigal (Hyatt et al., 2010) for ORF prediction and eggNOG-Mapper (Huerta-Cepas et al., 2017) for protein function annotation using the virus database. Results can be visualized using Proksee (CGView) (Grant & Stothard, 2008) and Clinker (Gilchrist & Chooi, 2021). Other than the method, other databases like Prokaryotic Virus Orthologous Groups (pVOGs) (Grazziotin, Koonin & Kristensen, 2017), and annotation tools like DRAM (Shaffer et al., 2023) and Pharokka (Bouras et al., 2023) can also be used to annotate the final contigs.

Contigs generated by FVE-novel

In this experiment, the raw assemblies was prepared by Tithi et al. (2023). To be more specific, the GOV database containing 24,411 contigs as the reference “genomes” were taken and FVE-novel (Tithi et al., 2023) was applied to an ocean metagenomic sample (NCBI; Sayers et al., 2022) accession number SRX2912986 (Aylward et al., 2017) to generate viral contigs (Tithi et al., 2018). The sample contained 18,471,506 paired-end reads with an average read length 151 bp. Through FVE-novel, we chose five contigs (hereafter labeled as S0, S1, S2, S3, and S4) and applied Virseqimprover to all of them to see whether the contigs could be either further extended and/or corrected for any error. Among the five contigs, S0, S1, and S2 were highly similar to each other whereas S3 and S4 were not.

To validate the results, we reassembled the contigs using the “Map to Reference” algorithm implemented in Geneious 11.0.4 (Kearse et al., 2012) together with multiple rounds of manual inspection and processing. Through this semi-automated process, we hope to examine whether there are multiple viral strains or species and if there are, whether a complete assembly of the dominant strain can be generated. Specifically, the metagenome reads were aligned to contigs S0, S3 and S4 using the “Low sensitivity/Fastest” setting allowing for 10% mismatches. Then the consensus sequence from the alignment was segmented into contigs with the highest depth of coverage >40x. These contigs were binned into lists of contigs with similar depth of coverage for further assembly. Next, the contigs in each bin were iteratively grown using Geneious by mapping reads to the ends with high stringency. To be more specific, all of the phage metagenome paired-end reads were aligned to these high depth of coverage contigs using “Map to Reference” with stringent “Custom Sensitivity” settings allowing no more than 1% “Mismatches per Read” and 1% “Gaps per Read” and requiring that both of the paired-end reads map to the new consensus sequence. This process was iteratively continued until the extended contigs merged together, maintained approximately uniform depth of coverage, and could no longer be extended. Using this laborious and semi-automated approach, we recovered a 153 kb contig from contig S0, a 177 kb contig from contig S3, and a 151 kb contig from contig S4, respectively.

Contigs generated by metaSPAdes and MEGAHIT

We also ran metaSPAdes (Nurk et al., 2017) and MEGAHIT (Li et al., 2015) on a marine plankton metagenome sample (NCBI Sayers et al., 2022) accession number SRX7079549 (Beaulaurier et al., 2020) collected from Station ALOHA in the North Pacific Subtropical Gyre to generate contigs and then applied Virseqimprover to see whether it can extend the contigs and correct any assembly errors. These two programs (metaSPAdes and MEGAHIT) have been shown to have less misassemblies compared to some other metagenome assemblers (e.g., IDBA-UD Peng et al., 2012) and Faucet (Rozov et al., 2018) as well as have good performances at the strain-level (Wang et al., 2020). However, many contigs generated by the two programs are highly fragmented due to uneven abundances or repeat regions (Wang et al., 2020). In this original study, the authors generated not only metagenomic sequencing data but also nanopore long read data which we can use to examine the performance of Virseqimprover (Beaulaurier et al., 2020). Based on the nanopore sequencing approach, many complete virus genome sequences were recovered. We BLASTed the contigs generated by metaSPAdes and MEGAHIT against the two recovered complete virus genomes AFVG_25M466 and AFVG_25M409, and identified three contigs that are highly similar to the genome sequences. The three contigs, referred to as S5, S6, and S7, were sent to Virseqimprover for further extension.

Evaluation of virseqimprover

A simulation experiment was conducted to evaluate the capability of correction and extension of Virseqimprover. We injected the genome of a papillomavirus (genome accession GCF_002826985.1) into the middle of contig S7 as a misassembly, generated 2,000 number of reads targeting 8X depth, spiked the reads into the original metagenomic sample that we used to generate S7 (SRX7079549), and ran MEGAHIT to assemble the spike-in sample. We identified the viral contig from the assemblies and ran Virseqimprover to evaluate whether it could correct remove the misassembly and extend it further.

Moreover, with regards to the state-of-the-art contig extension tool, ContigExtender is a popular tool that is used to extend assembled contigs following de novo assembly. ContigExtender uses a recursive overlap layout candidates strategy that explores multiple extending paths to achieve longer contigs (Deng & Delwart, 2021). Here we ran ContigExtender and compared the extended contigs with those extended by Virseqimprover using CheckV (Nayfach et al., 2021) for contig completeness and quality. CheckV determines the completeness and quality of assembled contigs by comparing them to IMG/VR (Camargo et al., 2023) database, which is a large database of complete virus genomes and has been used widely to evaluate the quality of the assembly.

Contigs S0, S1, S2

Among S0, S1, and S2, S0 is the longest with length 193,112 bp. Figure 3A shows that for S0, depth of coverage from around 24,500 bp to 25,500 bp is lower than the average depth of coverage and varies a lot after 150 kb. Virseqimprover checked for the uniformity of the depth of coverage and extracted the longest region with uniform depth of coverage, which was a region with length 127,423 bp from 25,567 bp to 152,989 bp. This longest non-suspicious region went through the iterative extension and error correction steps. When the iterative extension step was done, the output contig (length 163,662 bp) was checked for circularity. When the circular region was trimmed and Pilon was applied to the output contig, an improved contig (length 152,707 bp, denoted as S0′) was generated and became the final output of Virseqimprover as shown in Fig. 3A. Some missing part on the left was due to the trimming of the circular region.

Similarly, examination of the depth of coverage along S1 revealed a non-uniform depth of coverage region or suspicious region from 133,376 bp to 134,543 bp. Thus, Virseqimprover then extracted the longest non-suspicious region which was a region of length 133,375 bp from 1 bp to 133,375 bp and then applied the iterative extension and error correction steps to generate an extended contig. When all these steps were finished, the circularity of this contig was checked and Pilon was applied to the contig. The final output was a contig with length 152,362 bp. Figure 3A shows the depth of coverage along both original (S1) and final sequences (S1′).

For contig S2 with length 80,620 bp, the per base depth of coverage was checked and Virseqimprover identified no non-uniform depth of coverage region or suspicious region. Hence this contig directly went through the iterative extension and error correction steps. Finally, after using Pilon to improve the assembly, a contig of 151,828 bp was generated. Figure 3A shows that after applying Virseqimprover, we extended the original contig on both ends and nearly doubled the total length of the original contig to get a greatly extended contig (S2′) which has a uniform depth of coverage along the sequence.

The improved versions of contigs, S0′, S1′, and S2′, were compared to the 153 kb strain of a novel uncultured virus. After predicting the genes and visualizing the gene cluster comparison, as shown in Fig. 3B, we can see that S1′ and S2′ are very similar to the 153 kb reference strain, whereas S0′ is a bit different from all other contigs. This shows that Virseqimprover has correctly recovered the whole virus sequences. Additionally, since some regions in S0′ do not match with any of the regions of all the other contigs, it could be an indication that S0′ is a different strain of the same virus species.

Contig S3

In the same manner, S3, with length 132,604 bp was checked for the uniformity of the depth of coverage and found that it had a region with low depth of coverage at around 49 kb to 66 kb as shown in Fig. 4A. BLAST search shows that this region aligns best with Cyanophage P-RSM3 and Prochlorococcus phage P-SSM4 whereas the other parts of the contig aligns best with Cyanophage P-RSM1 and Synechococcus phage metaG-MbCM1, indicating that this low depth of coverage region is probably a misassembly. Virseqimprover flagged all the suspicious regions and extracted the longest region with uniform depth of coverage from 66,390 bp to 132,603 bp. Then this 66 kb region was extended through iterative extension and error correction and a contig with length 160,821 bp was generated. After checking for the circularity of this contig and applying Pilon to this contig, an improved contig (denoted as S3′) with length 160,744 bp was produced which was the final output of Virseqimprover. Figure 4A shows that the non-uniform regions in the original contig are filtered out so that the final corrected and extended contig has a uniform depth of coverage along its length. Using BLAST search, we identified four phages that show the highest sequence similarity to S3′. Based on the gene cluster comparisons, as shown in Fig. 4B, we can see that the improved contig S3′ does have some similar proteins with these phages. However, DNA sequence alignment of S3′ with these genomes also reveals some dissimilar regions, with great sequence identity variation along the entire sequence, ranging from 75.38% to 82.72%. Hence the improved S3′ might be from a novel phage species.

Figure 4C shows the similarity of S3′ (length 160,744 bp) to the viral sequence (177,631 bp) recovered by the semi-automated assembly process in Geneious. The protein identity threshold is 30%, which means that two proteins are considered to belong to the same group if their protein identity value is above 30%. Apart from the beginning part in the 177 kb strain that does not have many alignments in S3′, a small region in the middle of the sequence also shows difference between these two sequences (colored as the gray arrows). We thus further compared S3′ with the 177 kb strain to find out the difference in this specific region. Based on the pairwise DNA sequence alignment by EMBOSS Stretcher (Madeira et al., 2022), the comparison of these two contigs reveals that in the 177 kb strain, a 1,282 bp region from 56,831 bp to 58,112 bp does not have many matches with S3′. In this part instead of this 1,282 bp region, S3′ contains a 1,342 bp region from 22,163 bp to 23,504 bp. Analysis of the depth of coverage of these two sequences in the area where they are different reveals that those areas have a relatively lower depth of coverage (about 150x) compared to the average depth of coverage (about 300x). Moreover, based on the BLASTP search, the specific protein sequence corresponding to this area in contig S3′ aligns best with Synechococcus phage metaG-MbCM1, whereas the 177 kb strain aligns better with Synechococcus phage S-SM2. The differences between S3′ and the 177 kb strain suggest that they may represent different strains of the same phage.

Contig S4

S4 has 136,254 bps. After Virseqimprover’s contig extension, error correction, and circularity check, S4 was extended to 151,190 bps. Figure 5 shows that the depth of coverage of the original S4 is rather uniform, and was extended for both sides of the sequence, with both left and right ends of S4′ showing higher depth of coverage than nearby regions, which indicates the presence of repeats. Closer examination of the end sequences reveals that the regions indeed are repeats.

Contig S5, S6

S5 is generated by MEGAHIT and contig S6 is generated by metaSPAdes. These two contigs both have a 99% identity to the marine virus with ID AFVG_25M466, covering 12% and 40% of the viral genome, respectively. After applying Virseqimprover, S5 got extended from 4,179 bp to 14,374 bp, and S6 extended from 13,396 bp to 22,526 bp (Fig. 6). The extended S5′ covers 43% of the marine virus genome, and the extended S6′ covers 68%. Pairwise alignment between the extended sequences and the marine virus genome shows that the extended parts are identical to the corresponding parts of the marine virus genome, suggesting that Virseqimprover can accurately extend the contig sequences generated by other assemblers.

Contig S7

S7 is generated by metaSPAdes. It has a 99% identity to the marine virus with ID AFVG_25M409. After applying Virseqimprover, S7 was extended from 23,114 bp to 32,035 bp, which covers 98% of the AFVG_25M409 viral genome (32,812 bp). Figure 7 shows that Virseqimprover successfully extends the original fragmented contig to nearly complete genome with high accuracy.

Simulation experiment evaluation

In order to examine the capability of Virseqimprover to remove misassembly within sequences, a papillomavirus genome was injected into the sequence of original S7 at a random position. After generating and spiking in the 2,000 synthetic reads and running MEGAHIT to assemble the spike-in sample, a contig with length 31,745 was found to have 99.97% identity with Marine virus AFVG_25M409 and 100.00% identity with the papillomavirus, respectively. Then we ran Virseqimprover on the contig with the insertion and got a sequence with length 32,395 bp. This revealed that the injected genome was entirely removed and the pure AFVG_25M409 viral genome was fully recovered. Figure 8 shows the depth of coverages of the original and recovered sequences, and the low depth of coverage region within the original sequence corresponds to the injected viral genome.

Contig completeness and quality evaluation

Figure 9 shows the completeness of both the original and refined sequences by ContigExtender and Virseqimprover. Based on CheckV, only one of the eight contigs generated by either FVE-novel, metaSPAdes or MEGAHIT is complete. In contrast, four contigs recovered by Virseqimprover are complete, thus achieving 37.5% improvement over these assembly tools in contig completeness. ContigExtender did not generate any complete contigs (S0 is already 100% complete and ContigExtender extended another 256 bp). Remarkably, after applying Virseqimprover, the completeness of S2 improved from 48.19% to 100%, S4 from 60.72% to 100%, and S7 from 68.89% to 100%. S5 also shows a significant improvement of completeness from 13.18% to 44.92%, compared to 15.16% by ContigExtender. This shows that Virseqimprover is effective in extending contigs and improving the completeness of contigs generated by other assemblers or contig extension tools.

Table 1 shows the sequence quality assessed by CheckV as well as length information for all the contigs generated by the original assemblers and extended by either ContigExtender or Virseqimprover. As for the quality, CheckV has four categories of increasing quality: low, medium, high, and complete. Virseqimprover improved the quality of six of the eight contigs, compared to only one by ContigExtender. Moreover, four of the eight contigs extended/refined by Virseqimprover are complete, two high, one medium, one low, compared to two high, three medium, and three low for the original contigs, and two high, four medium, and two low by ContigExtender. Thus, Virseqimprover produces better quality and more complete sequences than the several commonly used assemblers as well as ContigExtender, a tool designed specifically to extend contigs. Interestingly, due to the nonuniform depth of coverage, Virseqimprover trimmed part of S0 that has low depth of coverage, which made the final sequence shorter. Taking together, Virseqimprover significantly enhances the quality of the contigs and is able to extend fragmented sequences into more complete sequences or full genomes. As ContigExtender and Virseqimprover are both contig extension tools, it warrants a detailed comparison between the two.

ContigExtender focuses on consensus read mapping and does not consider uniform depth of coverage and early stopping conditions. Also, ContigExtender runs only on Docker images, so it is not as flexible as Virseqimprover especially when users on HPCs do not have admin permissions. Instead, Virseqimprover not only provides a user-friendly web server but also a Conda package so that all users can run the tool easily on any platform. Another limitation of ContigExtender is that since the input contig is kept entirely for extension, if there is any misassembly within the sequence, ContigExtender will not correct it and extending the misassembled sequence longer only leads to longer misassembled sequence. In contrast, Virseqimprover considers error-correction and extension together when running the iterative local assembly process, and also trims some bases in the sequence during the iterative assembly process to avoid early stopping when the extension gets terminated. As a result, Virseqimprover will conveniently get longer, more complete and more accurate sequences for the input sequences.