Tutorial Objectives
Table of contents:
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!
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.
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.
Below is a screenshot of how to navigate through the SRA page to where the link is located.
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.
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.
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.
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 .
At this point, we're ready to identify what sequences are potentially viral. For this, we'll use VirSorter2.
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.
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.
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.
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!).
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.
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.
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.
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.
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.
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.
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...
And finding some interesting clusters!
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.
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:
And taking a quick peek into one of the above VCs...
It's unsurprising that some of the most abundant organisms in the oceans would have a virus!
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...
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.
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!
For further reading:
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.