Viral Proteins

Each metric is computed using sliding windows from 10 to 100 genes wide, starting at every gene along the sequence, and all scores greater than 2 are stored. Local maxima of significance score are then searched and the associated set of genes is defined as a putative viral region. These different predictions (based on the metrics above) are then merged when overlapping (extending the regions to include all predicted windows), leading to a list of putative viral regions associated with a (set of) metric(s). These regions are classified into three categories: (i) category 1 (“most confident” predictions) regions have significant enrichment in viral-like genes or non-Caudovirales genes on the whole region and at least one hallmark viral gene detected; (ii) category 2 (“likely” predictions) regions have either enrichment in viral-like or non-Caudovirales genes, or a viral hallmark gene detected, associated with at least one other metric (depletion in PFAM affiliation, enrichment in uncharacterized genes, enrichment in short genes, depletions in strand switch); and (iii) category 3 (“possible” predictions) regions have neither a viral hallmark gene nor enrichment in viral-like or non-Caudovirales genes, but display at least two of the other metrics with at least one significance score greater than 4. Finally, if a predicted region spans more than 80% of predicted genes on a contig, the entire contig is considered viral. A summary of VirSorter detection types is displayed in Fig. 1B.
Next, higher confidence predictions are used to refine the sequence space search. Specifically, sequences from all open reading frames from category 1 predictions that do not match a viral protein cluster are clustered and added to the reference database (RefSeqABVir or Viromes depending on the initial user choice). This updated database is then used in another round of search by VirSorter. This iteration where category 1 sequences are used to refine the searches is continued until no new genes are added to the database. Once no new genes are added, the final VirSorter output is provided to the user and includes nucleotide sequences of all predicted viral sequences in fasta files, an automatic annotation of each prediction in genbank file format, and a summary table displaying for each prediction the associated category and significance scores of all metrics. By providing the predictions and the underlying significance scoring, users can evaluate each prediction and apply custom thresholds on significance scores through a simple text-parsing script, even for large-scale datasets.
VirSorter is available as an application (App) in the iPlant discovery environment (https://de.iplantcollaborative.org/de/) under Apps/Experimental/iVirus (see Fig. S1 for a step-by-step guide of VirSorter app on iPlant). This application allows users to search any set of contigs for viral sequences using either the RefSeqABVir or the Viromes database. The reference values of VirSorter metrics will be evaluated on the complete set of input sequences, hence mixed datasets should be sorted (when possible) by type of bacteria or archaea in order to get the most accurate result possible. In addition to these reference databases, the VirSorter App on iPlant allows users to input their own reference viral genome sequence already assembled or to-be assembled using iPlant Apps prior to analysis with VirSorter. Assembled sequences are processed as follows: (i) genes are predicted with MetaGeneAnnotator (Noguchi, Taniguchi & Itoh, 2008 (link)), (ii) predicted proteins are clustered with sequences from the user-selected database (either RefSeqABVir or Viromes), and (iii) unclustered proteins are added to the “unclustered” pool. VirSorter scripts are also available through the github repository https://github.com/simroux/VirSorter.git.

Two reference databases of viral protein sequences were built for VirSorter and are available in the iPlant Discovery Environment (Data/Community_Data/iVirus/VirSorter/Database). The first includes 114,297 proteins from viruses infecting bacteria or archaea in RefSeqVirus genomes (as of January 2014), hereafter named “RefSeqABVir.” Protein clusters (PCs) were defined using MCL clustering (Enright, Van Dongen & Ouzounis, 2002 (link)) of these proteins (inflation 2.0) based on their reciprocal blastp comparisons (threshold of 50 on bit score and 10⁻⁰³ on E-value). The 9,735 PCs with at least 3 sequences were used to define a profile database searchable with HMMER3 tools (Eddy, 2011 (link)). The remaining 34,668 unclustered sequences were formatted for a blastp search. All PCs that did not contain any sequences from Caudovirales and unclustered sequences from viruses other than Caudovirales were marked as “Non-Caudovirales.”
The RefSeqABVir database was then augmented by virome sequences sampled from freshwater, seawater, and human gut, lung and saliva, resulting in an extended version of the reference database (hereafter named “Viromes”) which includes both virome and RefSeqABVir sequences. This combined reference dataset should help to detect new viruses for which no cultivated reference sequence is available. When only raw reads were available, viromes were assembled using Newbler (threshold of 98% identity on 35bp). The resulting contigs were then checked for the presence of cellular genome sequences, and only the 68 viromes for which no 16S rRNA genes were retained (see Table S1 for a complete list of these viromes). Contigs assembled from these 68 viromes were then manually inspected (through annotations generated by Metavir; Roux et al., 2014a (link)) and revealed no identifiable cellular genome sequences (i.e., no sequence contained more than 2 genes that matched a cellular genome and were not found in any known virus). A total of 146,521 complete predicted proteins from this quality-controlled dataset were then clustered with the 114,297 proteins from RefSeqABVir, leading to 15,673 clusters with 3 sequences or more, and 88,052 unclustered sequences. PCs from the combined Viromes database were used to create a profile database searchable with HMMER3, and the 34,338 unclustered sequences from RefseqABVir were formatted for BLAST search (unclustered sequences from viromes were not added to the database to prevent the inclusion of contaminating sequences).
Within these databases, viral “hallmark” genes were defined though a text-searching script looking for “major capsid protein,” “portal,” “terminase large subunit,” “spike,” “tail,” “virion formation” or “coat” annotations. After a manual curation step removing genes with more general annotation such as “protease” or “chaperone,” 826 PCs or single genes were identified as “viral hallmark genes.” This latter point meant removing domains also matching “protease” or “chaperone” domains and was conducted to minimize false positives for our viral hallmark genes category by extra-cautiously avoiding PCs that might include domains that could derive from either both viruses or microbes.

Predicted proteins from reference viral genomes from NCBI and VOG database viral proteins were combined to generate v-scores, which resulted in a total of 633,194 proteins. Redundancy was removed from the viral protein dataset using CD-HIT (v4.6) [56 (link)] with an identity cutoff of 95%, which resulted in a total of 240,728 viral proteins. This was the final dataset used to generate v-scores. All KEGG HMM profiles were used to annotate the viral proteins. A v-score for each KEGG HMM profile was determined by the number of significant (e-value < 1e−5) hits by hmmsearch, divided by 100, and a maximum value was set at 10 after division. The same v-score generation was done for Pfam and VOG databases. Any HMM profile with no significant hits to the virus dataset was given a v-score of zero. For KEGG and Pfam databases, any annotation that was given a v-score above zero and contained the keyword “phage” was given a minimum v-score of 1. To highlight viral hallmark genes, any annotation within all three databases with the keyword portal, terminase, spike, capsid, sheath, tail, coat, virion, lysin, holin, base plate, lysozyme, head, or structural was given a minimum v-score of 1. Non-prokaryotic virus annotations (e.g., reovirus core-spike protein) were not considered. Each HMM is assigned a v-score and represents a metric of virus association (i.e., do not take into account virus specificity or association with non-viruses) and are manually tuned to put greater weight on viral hallmark genes (Additional File 4: Table S4). Overall, annotations that are likely non-viral will have a low v-score whereas annotations that are commonly associated with viruses will have a high v-score. Raw HMM table outputs for v-score generation can be found in Additional Files 5, 6, and 7 for KEGG, Pfam, and VOG, respectively (Additional File 5: Table S5, Additional File 6: Table S6, and Additional File 7: Table S7).