Generated October 19, 2024

Introduction

This pipeline processes a single metagenomic dataset (downloaded from NCBI SRA) from the Global Ocean Viromes dataset, assembles it using metaSPAdes, identifies the viral sequences using VirSorter2, and subsequently classifies them with vConTACT2 and annotates with Prokka and DRAM-v. Along the way there are several intermediary steps to convert files from one KBase object to another - and though some of these can be skipped for this particular dataset - we're explicitly doing them here as to be more easily generalized to other datasets.
iVirus 2.0 on KBase
KBase Processing pipeline using iVirus 2.0 apps

Tutorial Objectives

  • Assemble viral metagenome using raw reads
  • Identify viral sequences from assembly data
  • Classify predicted viral sequences

Table of contents:

  1. Scientific background... and why viruses are important!
  2. Importing read data from the SRA
  3. Assessing read quality with FastQC
  4. Quality control of reads with Trimmomatic
  5. Post-QC read quality assesssment
  6. Assembling the data with metaSPAdes
  7. Identifying dsDNA and ssDNA phage using VirSorter2
  8. Annotating viral genomes with Prokka
  9. Classifying viral genomes using vConTACT2
  10. Post-vConTACT2 analysis
  11. vConTACT2 Quick Reference Guide
  12. Conclusions
  13. Version History
  14. Feedback and Help

Background - why study viruses?!

Viruses are important for more than a few reasons!

The "most abundant biological entities on the planet"1,2, with 1031 virus-like particles3

3% (0-18%) of any microbial genome is really virus4,5

On the topic of composition, 8% of the human genome is viral6,7

Move 1029 genes per day, globally 8,9

Lyse between 20-40% of ocean microbes daily10

They steal metabolic genes and can encode key metabolic components (like photosynthesis!)11,12,13

Can infect other viruses (virophages)14

Viruses - not microbes - encode many of the toxins we think of as bacterial (bordtella, cholera, shiga, etc) 15

Fewer than 1% are culturable16,17

Viruses are often hidden in datasets (both viral and microbial) this guide will help you find them!

Good luck!

Importing read data from the SRA

The reads in this dataset were generated from the Global Ocean Virome. and deposited as ERR594369. This is also known as Tara Oceans Expedition Station 36 surface ("SRF"), a coastal area in the Indian Ocean (more specifically, Northwest Arabian Sea), taken from ~5 m depth.

Image of Tara Oceans Expedition Station 36
Figure of Station 36 location. Figure heavily modified from Roux et al (2016) Nature

We first need to upload these reads into KBase. Files can be uploaded through a simple drag and drop interface (to the "left" of the Narrative, in the data window, or through Globus. A full guide on data upload and download can be found at http://kbase.us/data-upload-download-guide/. For this dataset, we'll be importing our reads directly from SRA into KBase. Alternatively, you could go to the link (more below), download the file(s) from SRA directly to your computer, and then upload them into KBase.

Let's stay simple and let KBase do all the work: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=ERR594369

The initial view (below) gives some details about the run, as well as the experiment and project it's a member.

Importing read data from the SRA
Figure of SRA ERR594369 run landing page

Below is a screenshot of how to navigate through the SRA page to where the link is located.

Importing read data from the SRA
Figure of link location to copy-and-paste into "Import SRA File as Reads" app in KBase

Note! Since the initial run of this Narrative, SRA has updated their storage servers. The link has changed from SRA to Amazon AWS. That said, the direct download link (below) is different - but the data retrieved does not change nor does it affect the results in any way, shape or form. Therefore, it was decided not to re-download the large raw files.

Import an SRA file from a web URL into your Narrative as a Reads data object.
This app completed without errors in 2h 15m 4s.
Objects
Created Object Name Type Description
ERR594369 PairedEndLibrary Imported Reads
Links

Assessing read quality

Following import, we always want to check the quality of the data going into an analysis. Unless you have supreme confidence in your viral isolation, extraction, sequencing preparation, and sequencing facility skills, it's always a good idea to know what quality is [eventually] going into an assembly. To quote a populat CS phrase, "Garbage In, Garbage Out." Essentially, this means that if you put poor quality data into your analysis, you're going to get poor quality results.

A quality control application for high throughput sequence data.
This app completed without errors in 31m 2s.
Links
Files
These are only available in the live Narrative: https://narrative.kbase.us/narrative/126189
  • ERR594369_75811_2_1.rev_fastqc.zip - Zip file generated by fastqc that contains original images seen in the report
  • ERR594369_75811_2_1.fwd_fastqc.zip - Zip file generated by fastqc that contains original images seen in the report

Quality control of reads

The FastQC of these enriched reads is already pretty good. You could get away with proceeding directly to assembly, but we'll trim reads to remove low-quality regions at the read ends and any residual adapters.

We'll trim using a popular read trimming tool, Trimmomatic . Defaults are okay, unless your read data has unique adapter or sequencing conditions.

Trim paired- or single-end Illumina reads with Trimmomatic.
This app completed without errors in 2h 8m 60s.
Objects
Created Object Name Type Description
Reads.t_paired PairedEndLibrary Trimmed Reads
Reads.t_unpaired_fwd SingleEndLibrary Trimmed Unpaired Forward Reads
Reads.t_unpaired_rev SingleEndLibrary Trimmed Unpaired Reverse Reads
Links

Post-QC read assessment

Following read trimming, we'll check to ensure that 1) all adapters are removed (if they even existed in the 1st place) and 2) that the reads are of sufficient quality for assembly.

A quality control application for high throughput sequence data.
This app completed without errors in 1h 53m 36s.
Links
Files
These are only available in the live Narrative: https://narrative.kbase.us/narrative/126189
  • Reads.t_paired_75811_5_1.fwd_fastqc.zip - Zip file generated by fastqc that contains original images seen in the report
  • Reads.t_paired_75811_5_1.rev_fastqc.zip - Zip file generated by fastqc that contains original images seen in the report

Assembly with MetaSPAdes

With clean reads, we'll now assemble these reads into contigs using another popular assembler, MetaSPAdes . Of the assemblers currently available in the KBase ecosystem, I'd argue MetaSPAdes performs slightly better than MEGAHIT and better than IDBA-UD. For an excellent review of viral benchmarks regarding assemblers, see this PeerJ article .

Assemble metagenomic reads using the SPAdes assembler.
This app completed without errors in 17h 14m 47s.
Objects
Created Object Name Type Description
SPAdes.contigs Assembly Assembled contigs
Summary
Assembly saved to: bbolduc:narrative_1603837393062/SPAdes.contigs Assembled into 18226 contigs. Avg Length: 5849.286678371557 bp. Contig Length Distribution (# of contigs -- min to max basepairs): 17860 -- 2000.0 to 31009.3 bp 249 -- 31009.3 to 60018.6 bp 67 -- 60018.6 to 89027.9 bp 20 -- 89027.9 to 118037.2 bp 13 -- 118037.2 to 147046.5 bp 7 -- 147046.5 to 176055.8 bp 3 -- 176055.8 to 205065.1 bp 0 -- 205065.1 to 234074.4 bp 4 -- 234074.4 to 263083.69999999995 bp 3 -- 263083.69999999995 to 292093.0 bp
Links

Identifying viral genomes using VirSorter2

At this point, we're ready to identify what sequences are potentially viral. For this, we'll use VirSorter2.

Overview of VirSorter2 framework
VirSorter2 Overview. Figure modified from Guo et al (2021) Microbiome

VirSorter2 will consume the KBase assembly object "SPAdes.contigs" and identify virus sequences (both integrated and non-integrated). It will generate another assembly object, which can be fed into downstream apps.

Tool performances on dsDNA phages from different data sources.
VirSorter2 performance compared to other tools against dsDNA phage. Figure modified from Guo et al (2021) Microbiome

VirSorter2 performs equivalent to OR superior than nearly all other tools at 3-5 kb and greater. That said, it is important to recognize that database and algorithm improvements often bias the most recently published tools. Always be aware of the databases used to train the model(s) and benchmarks employed by every virus identification tool. Many tools often only examine/benchmark using viral RefSeq, which is heavily biased towards human and human-affecting (e.g. commercial crops) viruses.

VirSorter2 default parameters are fine to use. If you want to use DRAM-v on the output, you will need to enable the "Enable DRAM-v outputs" flag. Also, be sure to include the virus groups. Since our data derives from a marine sample and has not undergone any type of RNA-focused extraction/amplification, we only need to use dsDNA phage and ssDNA. In this example, I've also enabled "Only output high confidence viral sequences." Other options are either situational or adjust the level of sensitivity or specificity in virus identification.

Identifies viral sequences from viral and microbial metagenomes
This app completed without errors in 3h 25m 59s.
Objects
Created Object Name Type Description
VirSorter2 Assembly KBase.Assembly object from VirSorter2
Summary
Results from your VirSorter2 run. Above you'll find a report with the identified,*putative* virus genomes, and below, downloadable links to the results files and links to the KBase assembly object. For users who enabled DRAM-v compatibility, the shock ID is 52a46a09-57c5-4601-a4f7-e1928e61b779
Links
Files
These are only available in the live Narrative: https://narrative.kbase.us/narrative/126189
  • final-viral-boundary.tsv.tar.gz - Table with boundary information
  • final-viral-combined.fa.tar.gz - Viral sequences in FASTA format
  • final-viral-score.tsv.tar.gz - Table with scoring information
Output from VirSorter2
The viewer for the output created by this App is available at the original Narrative here: https://narrative.kbase.us/narrative/126189

Evaluating identified virus quality

metaSPAdes generated 18226 contigs and VirSorter2 - with our parameters - identified 5445 putative virus genomes. As a comparison, running VirSorter2 without enabling high confidence genomes and a 0.7 probability yields 6909 virus genomes, an increase of 26%.

However, we can (and should) get a second opinion to evaluate the quality of these putative virus genomes.

For this, we can use CheckV. It's one of the only tools available to access quality, and definitely the only one on KBase.

Overview of CheckV
CheckV Overview. Figure from Nayfach et al (2020) Nature Biotechnology

To summarize, CheckV uses virus marker genes and knowledge about the basic structure of viruses to access overall quality by identifying ends, completeness, and how strongly genes of the contigs look viral.

Do not expect CheckV to agree 100% with any virus prediction tools. They were benchmarked with different viruses, they use different databases, and viruses are weird (don't expect them to always match your expectations!).

CheckV is a fully automated pipeline for assessing the quality of single-contig viral genomes, including identification of host contamination for integrated proviruses, estimating completeness for genome fragments, and identification of closed genomes.
This app completed without errors in 15m 39s.
Links
Files
These are only available in the live Narrative: https://narrative.kbase.us/narrative/126189
  • checkv_result.zip - File(s) generated by CheckV App

Annotate viral genomes with Prokka

vConTACT2 on KBase no longer requires gene prediction from 3rd-party tools. However, it can be useful to have a genome object with predicted genes in KBase. We'll use Prokka here because it does a decent job at prediction - and while it's not the best when it comes to viral genomes, it's good enough for our purposes.

Prokka's options should be adjusted to viral Kingdom. All other options can be left as-is.

Annotate Assembly and Re-annotate Genomes with Prokka annotation pipeline.
This app completed without errors in 9m 49s.
Objects
Created Object Name Type Description
CheckV Genome Annotated Genome
Summary
Annotated Genome saved to: bbolduc:narrative_1663261015855/CheckV Number of genes predicted: 52041 Number of protein coding genes: 51505 Number of genes with non-hypothetical function: 2676 Number of genes with EC-number: 564 Number of genes with Seed Subsystem Ontology: 0 Average protein length: 241 aa.
Output from Annotate Assembly and Re-annotate Genomes with Prokka - v1.14.5
The viewer for the output created by this App is available at the original Narrative here: https://narrative.kbase.us/narrative/126189

DRAM-v

DRAM-v will be used to functionally profile (i.e. annotate) the putative viruses as well as identify Auxiliary metabolic genes (AMGs). To cover AMGs is well beyond the scope of this pipeline, but to grossly simplify, AMGs are virus-encoded microbial metabolic genes that allow metabolic reprogramming of the infected host.

Overview of DRAM and pipeline
DRAM conceptual overview and pipeline. Figure from Shaffer & Borton et al (2020) Nucleic Acids Research
Annotate your viral assembly with DRAM. Annotations will then be distilled to create an interactive functional summary per virus.
This app produced errors in 5d 12h 38m 23s.
No output found.

Classify annotated viral genomes with vConTACT2

Now that we've handled virus identification and other minor details, we can take our assembly and use it in vConTACT2 .

vConTACT2 works by using a gene-sharing network to associate viral genomes. The more genes that are shared between two genomes, the higher the probability of those two genomes being phylogenetically related.

Image of vContact2's gene-sharing network
vConTACT2 virus classification. Figure heavily modified from Jang and Bolduc et al (2019) Nat. Biotech

So what's happening in the background? vConTACT2 will extract each viral genome and its associated gene predictions, and build the Gene2Genome table that underpins the whole analysis. Thankfully, vConTACT2+KBase generates this file in the background! For non-KBase users, this could be a challenge unless you let vConTACT2 handle everything.

There are a lot of options for vConTACT2. As a developer, there's a balance between giving enough options to allow for granular control of how the tool operates, and not over-burdening the user with options most are unlikely to change. In KBase, all the default options have been selected. There's no need to change anything - except if you want to use the most recent version of NCBI's Viral RefSeq. Often, users prefer to use the "older" version as that's what was used in the publication, so they're looking for consistency. If you'd like to use the most recent, then there might be very minor differences.

Viral cluster automatic cluster taxonomy
This app completed without errors in 1h 21m 56s.
Summary
Basic message to show in the report
Links
Files
These are only available in the live Narrative: https://narrative.kbase.us/narrative/126189
  • c1.ntw.tar.gz - ClusterONE network file suitable for import into Cytoscape or other graph tools
  • genome_by_genome_overview.csv - Final summary file generated directly by vConTACT2
Output from vConTACT2 0.9.19
The viewer for the output created by this App is available at the original Narrative here: https://narrative.kbase.us/narrative/126189

A closer look at the vConTACT2 results

After you've run vConTACT2, you'll get a table with ALL of the genomes in the analysis. The table can be a bit unwieldly as it contains a lot of rows and columns.

The easiest way to manage this (in KBase) is to use the filter function to find YOUR viral genomes. For example, all of our genomes contain NODE - it's a byproduct of the SPAdes assembler. By adding NODE in the appropriate filtering row under the column "Genome" will filter out all the reference genomes.

Image of NODE filtering on column
Filtering vConTACT2 results table using "NODE"

There are 5445 viral genomes remaining. This is great, why? It matches the number of genomes that were annotated by Prokka and found in the VirSorter2 assembly. Now, how many of our viral genomes are clustered? Add Clustered to the "VC Status" column. There are 2564 Clustered (or Clustered/Singletons) genomes. Not bad, not great - but this is actual data - not pretty "mock" data.

Image of NODE and VC Status filtering on column
Filtering vConTACT2 results table using "NODE" and "Clustered"

Now let's find out if our viral genomes are associated with any references. This is not fast through the table in KBase, but it is doable, unless you want to download the csv and do some data wrangling. What I do is sort the table by "VC" and scroll through, keeping track of the "Size" and count of the VC. If the "Size" of the VC is greater than the counts of the VC, make note of that VC.

After sorting by VC...

Image of NODE and VC Status filtering on column, then sorted
Sorting vConTACT2 results table after using "NODE" and "Clustered"

And finding some interesting clusters!

Image of NODE and VC Status filtering on column, then sorted
Finding VCs of interest after filtering, sorting, and comparing

One of the first examples we encounter is VC_277_0. It has a VC size of 6, but only 1 member is seen. Remember - we still have the `NODE` filter on! Remove that NODE filter, revealing the Beetreviruses.

Image of VC 277 with only environmental sequence
Revealing VC 277 with "NODE" filter on, revealing its members
Image of VC 277
Revealing VC 277 without "NODE" filter on, revealing its members

This is excellent, as the Cyanophage and Synechococcus phage are incredibly common in the ocean. A literature search revealed that these sequences are indeed found at the station our SRA reads are derived!

Additional searching reveals at least 6 other VCs that include both reference data and environmental sequences:

  • VC_138_0
  • VC_359_0
  • VC_367_0
  • VC_375_0
  • VC_843_0
  • VC_1079_0

And taking a quick peek into one of the above VCs...

Image of VC 367
VC 367 and its Prochlorococcus and Cyanophage phage members

It's unsurprising that some of the most abundant organisms in the oceans would have a virus!

A quick reference guide to vConTACT2</h2>

What is a "VC"?

A Viral Cluster is a unit of classification. vConTACT2 uses two terms - with subtle differences - in order to classify sequences.

A VC is a first-pass classification of viral genomes. This frequently represents a group of genomes within the same genus. An example is VC_115.

A VC Subcluster is a second-pass classification. It uses a pre-calculated distance calculation to refine the VCs. These are high-confidence, genus-level groupings. It is common for there to be no change between VC and VC Subcluster. An example is VC_115_0. The final value ("0") represents the subcluster within the original VC. If all members of a VC Subcluster have VC_xx_0, then there was no change. However, if there were further refinements, then the VC Subcluster would be: VC_115_0, VC_115_1, VC_115_2...

VC Status

How vConTACT2 describes how it classified a sequence

Clustered: Genomes "successfully" placed into a genus-level group alongside at least one other genome.

Singleton: Genome was not found to be related to any other genome in the dataset. Most likely reason? Very weak or no overlap with any other genes found on any other genome. How to fix? Add more related genomes.

Clustered/Singleton: Genome was initially clustered, but distance-based optimization identified its placement in the cluster as not genus-level. However, no other genomes were found to be within the same "subcluster" as this genome, and resulted in the genome being stranded, without another member. For these genomes, you can look at its "VC" to see distantly related members. So a viral genome that is Clustered/Singleton (VC_221_1 or VC_221_2) is related to other VC_221 members, but vConTACT2 does not have confidence that these genomes are related at the genus level.

Overlap (VC_NN/VC_XX): Genomes identified as sharing significant portions of its gene content with multiple VCs. In other words, vConTACT2 cannot confidently assign it to one OR the other VC. This is incredibly common for viral groups that undergo extensive recombination.

Outlier: Genomes clustered by ClusterONE (a tool used internally by vConTACT2) but were not strongly connected to the other VC members. It is common for these genomes to share a single gene or two to the other members of its closest VC, however the other members in that VC are likely be be sharing 20, 30 or 50+ genes. It is not only unlikely that that particular genome is related at the genus level, but could perhaps be a spurious shared gene and is unlikely to be related at anything lower than family or order.

Conclusions

In summary, we've processed a viral metagenome from public reads available on SRA, identified contigs from the assembled sequence data as putative viruses using VirSorter2, and classified them in approximately genus-level clusters with vConTACT2. We've also assessed the quality of the virus genomes using CheckV and gotten some information about them through Prokka and DRAM-v. This analysis revealed 7 VCs where environmental sequence data was found associated with references, and we can have confidence that those sequences at related to those references at the genus level.

Additionally, we've seen a few larger clusters with no associations to reference sequences that could be further investigated using KBase tools. For example - align those genomes in those VCs, identify shared features, extract those features, and make a discovery about certain proteins found throughout your dataset. Or, go one step further and pull in JGI data and then compare against a variety of JGI datasets for global significance - all using existing KBase apps!

Thanks for following along and using these apps for your research!

For further reading:

  • Guo, J. et al. VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses. Microbiome 9, 37 (2021). doi.org/10.1186/s40168-020-00990-y
  • VirSorter2 manuscript
  • Shaffer, M. et al. DRAM for distilling microbial metabolism to automate the curation of microbiome function. Nucleic Acids Res. 48, 8883–8900 (2020). doi:10.1093/nar/gkaa621
  • DRAM and DRAM-v
  • Bin Jang, H. et al. Taxonomic assignment of uncultivated prokaryotic virus genomes is enabled by gene-sharing networks. Nat. Biotechnol. 37, 632–639 (2019). doi:10.1038/s41587-019-0100-8
  • The significantly improved version of vConTACT that's faster, more accurate, and capable of handling larger datasets
  • Roux, S., Enault, F., Hurwitz, B. L. & Sullivan, M. B. VirSorter: mining viral signal from microbial genomic data. PeerJ 3, e985 (2015). doi:10.7717/peerj.985
  • The original VirSorter paper
  • 1. Roux, S., Enault, F., Hurwitz, B. L. & Sullivan, M. B. VirSorter: mining viral signal from microbial genomic data. PeerJ 3, e985 (2015). doi:10.7717/peerj.3243
  • The initial vConTACT paper describing the theory and background

Known Issues & Update History

Known Issues

  • None known

Update History (YYYY-MM-DD)

2022-09-13

  • Initial draft of tutorial built

Feedback & Helpdesk

Was this Narrative helpful? Please provide feedback and let us know:

If you have a question about one of our apps, need to report a bug or have another system-related query, please join our Help Board and post a ticket. Learn about how to do this here: http://kbase.us/help-board.

Released Apps

  1. Annotate and Distill Viral Assemblies with DRAM-v
    • DRAM source code
    • DRAM documentation
    • DRAM Tutorial
    • DRAM publication
  2. Assemble Reads with metaSPAdes - v3.15.3
    • Nurk S, Meleshko D, Korobeynikov A, Pevzner PA. metaSPAdes: a new versatile metagenomic assembler. Genome Res. 2017; 27:824 834. doi: 10.1101/gr.213959.116
    • Prjibelski A, Antipov D, Meleshko D, Lapidus A, Korobeynikov A. Using SPAdes De Novo Assembler. Curr Protoc Bioinformatics. 2020 Jun;70(1):e102. doi: 10.1002/cpbi.102.
  3. Assess Read Quality with FastQC - v0.12.1
    • FastQC source: Bioinformatics Group at the Babraham Institute, UK.
  4. Import SRA File as Reads From Web - v1.0.10
    • Arkin AP, Cottingham RW, Henry CS, Harris NL, Stevens RL, Maslov S, et al. KBase: The United States Department of Energy Systems Biology Knowledgebase. Nature Biotechnology. 2018;36: 566. doi: 10.1038/nbt.4163
  5. Trim Reads with Trimmomatic - v0.36
    • Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30: 2114 2120. doi:10.1093/bioinformatics/btu170

Apps in Beta

  1. Annotate Assembly and Re-annotate Genomes with Prokka - v1.14.5
    • Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30: 2068 2069. doi:10.1093/bioinformatics/btu153
  2. CheckV
    no citations
  3. vConTACT2 0.9.19
    • Bin Jang, H., Bolduc, B., Zablocki, O., Kuhn, J. H., Roux, S., Adriaenssens, E. M., Sullivan, M. B. Taxonomic assignment of uncultivated prokaryotic virus genomes is enabled by gene-sharing networks. Nature Biotechnology. 2019;37(6): 632 639. https://doi.org/10.1038/s41587-019-0100-8
  4. VirSorter2
    • Guo, J. et al. VirSorter2: a multi-classifier, expert-guided approach to detect diverse DNA and RNA viruses. Microbiome 9, 37 (2021). https://doi.org/10.1186/s40168-020-00990-y