The use of informativity in the development of robust viromics-based examinations

Materials and Methods

Development of the informativity metric

Establishing a taxonomic signal threshold

To ascertain the presence/absence of a specific taxon within a metagenome, we suggest a threshold to differentiate between informative and uninformative hits. The taxonomic signal threshold T is determined through a two-step process prior to evaluation of the metagenomic data. In the first step, each annotated coding region for a given taxon of interest is compared to all annotated sequences within the genome(s) of a known relative. Thus, each coding region’s sequence x (x ∈ X, where X is the set of sequences for all coding regions annotated within the genome of the taxon of interest) is compared to each coding region’s sequence g (g ∈ G, where G is the set of sequences for all coding regions annotated within the genome of a known relative). The use of a known relative genome(s) establishes if and how conserved the coding region is between known, related strains/species. Where sequence homology is detected, the sequence identity and query coverage of the match is recorded: S₁ and Q₁, respectively.

In the second step, each coding region’s sequence is compared again, this time to the sequences for all annotated coding regions for the group assayed by the metagenome (e.g., all phages, viruses, bacteria, archaea, etc.), however, those belonging to the taxonomic group containing the taxon of interest and the known relative considered in step one are omitted. Many hits may be recorded for a particular gene x. Thus the best hit, the highest scoring hit both with respect to the sequence identity and the query coverage of the match, is selected; S₂ and Q₂ denote this best match’s sequence identity and query coverage, respectively. A taxonomic signal threshold T is defined as T = {S₁−S₂, Q₁−Q₂} where the subscripts 1 and 2 represent the sequence identity and query coverage of the match detected from steps one and two, respectively. Figure 1 illustrates the two-step process and the T values produced.

It is important to note that the taxonomic group used for comparison is user defined. For instance, in order to ascertain if a gene can be used to distinguish between the presence/absence of a particular species, one may consider the taxonomic group to be inclusive only of strains of the species. Therefore, in this case, the most distant relative belonging to the taxonomic group in step one would be the closest related species. If a more distant relative, say the most distantly related species of the same genus, were to be investigated, then the taxonomic signal threshold T would serve as a means to distinguish between the presence/absence of a subset of the species (inclusive of the taxon of interest) within the genus. This flexibility enables the researcher to define and control the granularity of his/her analyses. If a particular taxa of interest lacks available genomes capturing the phylogenetic diversity of the species (or genus or subfamily, etc.), a more distant relative can be selected. In addition to the intended purpose of establishing the taxonomic signal threshold, the two-step process can provide insight into putative horizontally acquired elements and gene loss events, e.g., instances in which the gene did not include a homolog in the most distant relative but did exhibit sequence similarity to a gene within the genome of another taxonomic group.

Implementation

The method for assessing the informativity of viromic hits was implemented using a series of BLAST databases and BLAST searches. A collection of all coding regions (nucleotide sequences) for the taxon of interest (X) and all genes (amino acid sequences) annotated within the genome of the selected relative (G) are supplied by the user. A local BLAST database is created for G, and the genes belonging to X are queried against the local database via blastx. The sequence identity and query coverage of the match detected for the best hit for each gene is then parsed from the BLAST results quantifying each gene’s S₁ and Q₁ values. Next, a BLAST database is created for the annotated coding regions (amino acid sequences) provided for step 2 of this method (set Z), again supplied by the user. Each of the genes for the taxon of interest X is queried against this second local database via blastx; the results are again parsed for each gene’s S₂ and Q₂ values so that the taxonomic signal threshold T can be calculated.

A metagenomic dataset can next be evaluated, comparing each read or contig against a collection of annotated gene sequences. To accommodate the variation between characterized sequences in databases and environmental samples, contigs are translated—generating all six open reading frames—and a protein database representative of the metagenomic dataset is produced. Each BLAST hit is next assessed with respect to its scores {S_H, Q_H} relative to that of the gene’s threshold T. For each gene in the genome of interest X, the values for S₁, Q₁, S₂, Q₂, S_H, and Q_H are written to file. The user can then evaluate the likelihood of a particular taxon’s or taxonomic group’s presence within the metagenomic sample based upon the I values for informative genes. Note that for the analyses presented here we have weighted S and Q values equally in the calculation of T; the two values are, however, reported separately such that users can select their own weighting of the contributions of sequence identity and query coverage.

The described process has been automated via functionality developed in C++ (available for both Windows and Unix OS). Users must supply or specify the FASTA format files for the taxon of interest (X), the known relative (G), and the group assayed (less the taxonomic group of interest) (Z). If metagenomic comparisons are to be conducted, as this is optional in the current implementation, the user must also supply the metagenomic dataset. The code has been designed for both ease of use, speed, and flexibility, such that analyses can be tailored to the environmental niche and/or hypothesis under investigation. Most importantly, this is a light-weight solution which can be integrated into the standard method of viral metagenomic analyses. Source code, documentation, and sample data are publicly available at https://github.com/putonti/informativity.

Datasets examined

Identifying informative genes

The new metric described here, Informativity or I, provides a quantifiable means of identifying if a particular taxonomical group is present/absent within a sequenced community. Developed specifically for the detection of viral sequences in complex community metagenomic data sets, I captures the likelihood of a sequence belonging to taxa. Described in greater detail within the Methods, Fig. 2 provides an overview of how informative genes are identified. Users must supply the query sequence(s) (likely a contig or set of contigs from a sequenced community), at least two representative sequences for a taxon of interest, and lastly a set of sequences representative of ‘non-relatives’ (sequences belonging to other taxa of, e.g., viruses). The taxon of interest can be, e.g., a species, a genus, or a subfamily.

Informative genes for Caudovirales taxa

All protein coding sequences were collected for species belonging to 70 tailed-virus (Caudovirales) taxa identified by NCBI Taxonomy (see Methods). Using the ICTV type species as a representative of the taxa, each gene sequence (x) of the type species’ genome (X) was compared to all other gene sequences for species of the same taxa. For each gene, the sequence identity S₁ and query coverage of the match Q₁ for the most dissimilar homologous gene sequence within the taxa is calculated. This captures the sequence variation for the gene within the species of the taxon. Thus, the S₁ and Q₁ scores for one gene x _i may be from homology detected in one species of the taxa, while the scores for another gene x _j may be to a homolog in another species’ genome. If the gene is unique to the type species’ genome, then S₁ = Q₁ = 0. The sequence identity, S₂, and query coverage, Q₂, scores were next calculated for each gene in the type species’ genome; each gene was compared to: (1) genes belonging to species classified within other genera within the order Caudovirales and (2) genes belonging to species of other taxonomic orders. In contrast to the S₁ and Q₁ scores, the S₂ and Q₂ scores are for the best hit or the most similar homolog found. Using these two values, the taxonomic signal threshold T can be calculated (see ‘Methods’). This threshold value signifies how reliable the particular gene is as an indicator of the presence/absence of the species. Genes which are found in multiple species and taxa would thus have a low threshold value T and perform poorly as an indicator of the taxon.

Figure 3 illustrates the thresholds for Myoviridae and Podoviridae type species; Siphoviridae is included in Fig. S1. (Type species names and accession numbers as well as scores are listed in Table S1). In these maps, each gene’s taxonomic signal threshold is shown; dark gray boxes indicate uninformative genes; these uninformative genes either exhibit greater homology to species belonging to other phage taxa or lack homology to other representative genomes of the taxon of interest (i.e., are present only within the type species’ genome). Also listed for each taxon is the number of genome sequences included in the comparisons. Those taxonomical groups with more phylogenetic diversity represented within available genome sequences tend to have less informative genes. This is quite prominent when evaluating the 10 Podoviridae taxon: the well sampled subfamily of Autographivirinae species have significantly less informative genes than the undersampled Podoviridae genera of, e.g., F116virus and Bpp1virus. It is important to note, however, that taking into consideration numerous genome sequences does not necessarily mean that the phylogenetic diversity of the taxon was examined. In contrast to classifying unknown sequences by a single marker, the informativity metric provides a multiple gene marker strategy. Thus, taxonomical ‘calls’ for a sequence can be made with greater confidence by reporting the aggregate of informative markers found, not just the presence/absence of a single gene.

Targeting specific phages in environmental samples

The Pseudomonas phage PB1 was selected for examination. Each gene annotated for the PB1 genome (Accession Number: NC_011810) (Ceyssens et al., 2009) was compared first to the set of genes for the most distant relative of PB1 within its genus Pbunavirus (previously called Pbunalikevirus), Burkholderia phage BcepF1 (Accession Number: NC_009015). For each gene the S₁ and Q₁ values were computed. Next, all 93 annotated PB1 genes were compared to all genes from phage species—other than those classified as Pbunaviruses (see ‘Methods’). Homologous sequences were identified, the S₂ and Q₂ values. The similarity observed (the S₂ and Q₂ values) for each of the PB1 genes is shown in the heatmap of Fig. 4. Several PB1 gene sequences (as indicated by the color scale) exhibited sequence homology to genes within phage genomes of other taxa. Dark gray blocks in the heatmap signify that no homologs were detected. The upper chart in Fig. 4 details the percent sequence identity (bars) and percent query coverage (circles) values observed for the best hits to GenBank records. PB1 genes with homologies to other phage taxa include conserved genes (e.g., gp47 = tail fiber component and gp50 = DNA ligase), amongst other conserved “hypothetical proteins”.

The methodology developed here was then applied to seven freshwater DNA viromes (Table 1); a list of the SRA datasets from each study is provided in Table S3. Each of the 56 samples examined was first assembled (see ‘Methods’ for details). The PB1 coding regions were then compared to the 56 collections of contigs. The heatmap shown in Figure 5A graphically represents these results; each row represents a single sample (Methods). Again, each gene’s best hit within each virome’s sample was qualified (colored) with respect to its conservation amongst the Pbunaviruses, the gene’s S₁ and Q₁ value. Nevertheless, not all genes provide an equal signal as to the presence or absence of PB1 within the sample: some serve as better markers. As shown in Fig. 4, there are several “non-Pbunavirus” species which contain homologs to PB1 genes. Thus, the informativity I of each BLAST hit within the seven viromes was calculated. In doing so, individual genes that provide a strong signal for the Pseudomonas phage PB1 can readily be identified. Figure 5B represents the results of this computation, in which each hit to a PB1 gene is now assessed in light of the taxonomic signal threshold T.

In an effort to assess the strength of the metric presented here, we evaluated the raw BLAST results of the datasets and a BLAST score-based analysis. The BLAST results of Viromes II, IV, V, and VII are publicly available through the web service MetaVir (Roux et al., 2014). Nine of the samples from Virome I are also available through MetaVir. It is important to note that in contrast to the uniform method in which the viral metagenomes were preprocessed here (see ‘Methods’), the sequences submitted to MetaVir, or comparable online resources, may be assembled or raw sequences. Furthermore, MetaVir conducts BLAST comparisons against the RefSeq viral database (O’Leary et al., 2016), whereas here we have included all partial and complete phage sequences from GenBank. Nevertheless, hits to the Pbunavirus (Table S2) genomes were identified in all five MetaVir datasets; the Lake Michigan and Lake Bourget samples (nine samples from Virome I and both samples from Virome II, respectively) produced the most hits in MetaVir to the Pbunavirus genomes (hundreds to thousands), many which were the best hits identified. Hits from MetaVir metagenomic samples, including Viromes I, II, IV, V, VII and additional sampling sites not included in our proof-of-concept work, to the Pseudomonas phage PB1 genome are shown in Fig. S2.

Virome I, the Lake Michigan viral metagenomes generated by our group (Watkins et al., 2015; Sible et al., 2015), includes many informative genes (Fig. 5B) indicative of the presence of a Pbunavirus similar to PB1. Thus, with confidence, one can predict its presence within this sample. Viromes II, V, and VII contain far fewer hits to informative genes (one, two, and one PB1 genes respectively). Furthermore, their informativity scores are low, {S_H, Q_H}≈ T. This would suggest that PB1 (or a close relative) is not present within the sample: rather a homolog of the gene is present, within an uncharacterized species. As viral sequence databases expand through the isolation and characterization of additional viruses, the threshold T is likely to change thus providing greater confidence in the evaluation of BLAST hits for OTU calling.