Identifying Viral Contigs II

Overview

Teaching: 0 min
Exercises: 180 min

Objectives

Identify phage contigs in an assembled virome

Interpret geNomad’s output

Assess assembly completeness

Select medium-complete and low-contaminated contigs

Select viral contigs for further analysis

In this section, we will identify viral sequences among our assembled contigs. We will use geNomad, which identifies both viruses and plasmids by combining a gene-marker classifier with a neural network operating on the nucleotide sequence. geNomad also assigns a taxonomy to the sequences it calls viral, so we will test how well it performs ourselves by running it on our virome contigs.

Identify viral contigs

We will start to identify viral contigs with geNomad.

Challenge

Why is it important to use tools like geNomad on viromes?

Note that geNomad needs quite some memory to run, so allocate at least 20 GB for the sbatch job. Allocate 10 threads. geNomad runs on ordinary CPU nodes, so in the parameter --partition of the sbatch script you can use short,standard; no GPU allocation is needed. geNomad should run in a few minutes on our dataset.

Exercise

Run geNomad using the assembly as input:

Read and interpret the output following geNomad’s quickstart page
source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh
conda activate genNomad_v1.11.2
genomad end-to-end -t 20 <assembly> <output_directory> /veodata/03/databases/geNomad/v1.11.2/genomad_db
geNomad writes its results into several subfolders of the output directory. The summary we will work with is <outdir>/assembly_summary/assembly_virus_summary.tsv, which lists one row per sequence that geNomad classified as viral, together with its length, topology, virus score, number of hallmark genes and taxonomy. The corresponding sequences are in
Inspect the output files
after running geNomad, and answer the questions below.

```bash
head <output>

sbatch script for geNomad

#!/bin/bash
#SBATCH --tasks=1
#SBATCH --cpus-per-task=10
#SBATCH --partition=short,standard
#SBATCH --mem=20G
#SBATCH --time=2:00:00
#SBATCH --job-name=genomad
#SBATCH --output=10_genomad/genomad.slurm.%j.out
#SBATCH --error=10_genomad/genomad.slurm.%j.err

# this will run very quickly

source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh && conda activate genNomad_v1.11.2

# geNomad needs three positional arguments: the input assembly, the output
# directory, and the path to the geNomad database.
assembly='../1.2_assembly/10_assembly_flye/assembly.fasta'
outdir='10_genomad'
database='/veodata/03/databases/geNomad/v1.11.2/genomad_db'

genomad end-to-end -t 20 $assembly $outdir $database

Questions

Do the results corroborate your expectations?

Was a contig classified as not viral? Take one of these contigs from the dataset and BLAST it using blastn. What are the top hits? Are they expected?

How many viral contigs are there in the virome?

bash command for getting the number of phage contigs

geNomad’s virus summary contains only the sequences it called viral, so the number of viral contigs is simply the number of data rows in that file (the first line is a header):
tail -n +2 10_genomad/assembly_summary/assembly_virus_summary.tsv | wc -l
To see how the whole assembly was split between chromosome, plasmid and virus, use the aggregated classification instead:
cut -f 2- 10_genomad/assembly_summary/assembly_aggregated_classification.tsv | head

Estimating genome completeness

There are tools to assess the completeness of bacterial and viral genome sequences. For viruses we use CheckV. The tool identifies genes on the query sequence and compares them to a database of viral and bacterial marker genes. Each taxonomic group of bacteria or phages, e.g. a species or a family, has certain marker genes on its genome. So based on the number and types of marker genes, CheckV can figure out to which taxon a query contig belongs to and estimates how much of the genome it represents (estimated completeness). Note that, we will improve the taxonomic annotation in the “Viral Taxonomy and Phylogeny” section next week. CheckV also checks for unexpected marker genes on the sequence (estimated contamination), and whether part of a phage contig likely represents bacterial sequences, in which case the fragment could be part of a host genome with an integrated prophage.

Allocate 20 threads and at least 20 GB memory for the sbatch job. It should take ~1 minute to run.

# create a folder for the assessment (or let sbatch create it when you assign the output and error log files)
$ mkdir 20_results_assessment_checkv

# activate the conda environment containing the checkv installation
$ source /vast/groups/VEO/tools/anaconda3/etc/profile.d/conda.sh && conda activate checkv_v1.0.1

# run checkV on both assemblies
$ checkv end_to_end ...

sbatch script for running checkV

#!/bin/bash
#SBATCH --tasks=1
#SBATCH --cpus-per-task=20
#SBATCH --partition=short,standard
#SBATCH --mem=20G
#SBATCH --time=02:30:00
#SBATCH --job-name=checkv
#SBATCH --output=20_checkv/checkv.slurm.%j.out
#SBATCH --error=20_checkv/checkv.slurm.%j.err

# run CheckV to assess the completeness of single-contig virus genomes.
# First, activate the conda environment which holds the CheckV installation on draco:
source /vast/groups/VEO/tools/anaconda3/etc/profile.d/conda.sh && conda activate checkv_v1.0.1

# CheckV parameters (https://bitbucket.org/berkeleylab/checkv/src/master/#markdown-header-running-checkv)
# checkv end-to-end runs the CheckV pipeline from end to end :). It expects an input fasta file 
# with the assembly and an output path.
# -t: threads

# assigning variables for readability

database='/veodata/03/databases/checkv/v1.5'
outdir='20_checkv'
assembly='../1.2_assembly/10_assembly_flye/assembly.fasta'

checkv end_to_end -t 20 -d $database $assembly $outdir

Go through the CheckV results

CheckV produces many output files, and also saves files for the intermediate steps of the tools that are used to find viral marker genes (diamond and hmmsearch). A summary of all results can be found in the file quality_summary.tsv. Open the file and familiarize yourself with the information presented in the table. On CheckV’s website, you can find information about the output in the sections “How it works” and “Output files”.

How are completeness and length related? (qualitative answer, name examples)

What are possible reasons for contigs with low completeness to appear in the assembly?

Filter contigs

Exercise - select high-quality phage contigs

Use geNomad’s virus predictions and CheckV’s estimates of completeness and contamination to separate phage from non-phage contigs and reduce our assembly to high-quality contigs of high completeness. If you are an experienced programmer and have enough time, write (a) script(s) for that. If not, use the solutions below directly. For the solutions below, we used the file assembly_virus_summary.tsv from geNomad’s summary folder, which already contains only the sequences classified as viral. The results of CheckV are located in the file quality_summary.tsv located in CheckV’s output folder.

As a rule of thumb, you could keep all contigs with completeness >50% and contamination <5%. These values could be changed depending on the data and on the project. Note that filtering for low completeness can remove some conserved regions. Allocate 2 threads and 1GB memory for this job. It should take only a few seconds to run.

Note that geNomad renames sequences in which it detected a provirus: such a row has a seq_name like contig_1|provirus_1200_15600 instead of contig_1. Strip everything from the | onwards to get back the original contig name, otherwise those contigs will never match the CheckV table.

To read the tabular files, you can use python’s pandas package. Then, you can load a .tsv file and select content from it like this:
import pandas as pd

# the .tsv format separates cells by a tab ('\t')
genomad_df = pd.read_csv(genomad_results_path, sep='\t')

# every row in the virus summary is a viral prediction; split off the provirus suffix
genomad_selection = {row['seq_name'].split('|')[0] for index, row in genomad_df.iterrows()}

# do the same for the checkv results and use set operations to get the contigs selected by both tools
joint_selection = genomad_selection.intersection(checkv_selection)

# modify the code from yesterday (rename and filter contigs) to go through the assembly and save contigs in the joint selection
...

python script for selecting viral contigs

import os, sys
import pandas as pd
from Bio import SeqIO

def main():

    assembly_path = os.path.abspath(sys.argv[1])
    assert assembly_path.endswith(".fasta")

    genomad_results_path = os.path.abspath(sys.argv[2])
    assert genomad_results_path.endswith(".tsv")
    
    checkv_results_path = os.path.abspath(sys.argv[3])
    assert checkv_results_path.endswith(".tsv")
    
    out_fasta = os.path.abspath(sys.argv[4])
    assert out_fasta.endswith(".fasta")

    # read the tsv files as pandas dataframes
    genomad_df = pd.read_csv(genomad_results_path, sep='\t')
    checkv_df = pd.read_csv(checkv_results_path, sep='\t')
	
    # collect the sets of contigs which stick to our selection cutoffs
    genomad_selection = {row['seq_name'].split('|')[0] for index, row in genomad_df.iterrows()}
    checkv_selection = {row['contig_id'] for index, row in checkv_df.iterrows() if row['completeness'] > 50 and row['contamination'] < 5}
	
    # use set operation union to get the contigs in the geNomad set AND in the checkv set
    joint_selection = genomad_selection.intersection(checkv_selection)
	
    # print some numbers
    print(f"Predicted viral contigs: {len(genomad_df.index)}, selected by geNomad: {len(genomad_selection)}, selected by checkv: {len(checkv_selection)}, joint selection: {len(joint_selection)}")
	
    # define list of records to keep and fill it by comparing the contig id of each record to the joint set of selected contigs
    out_records = []
    with open(assembly_path) as handle:
        for record in SeqIO.parse(handle, "fasta"):
            if record.id in joint_selection: out_records.append(record)
    
    # write the selected records into a new file
    with open(out_fasta, "w") as fout:
        SeqIO.write(out_records, fout, "fasta")

if __name__ == "__main__":
    main()

sbatch script for submitting the python script

#!/bin/bash
#SBATCH --tasks=1
#SBATCH --cpus-per-task=2
#SBATCH --partition=short
#SBATCH --mem=1G
#SBATCH --time=00:30:00
#SBATCH --job-name=filter_contigs
#SBATCH --output=30_filter_contigs/filter_contigs.slurm.%j.out
#SBATCH --error=30_filter_contigs/filter_contigs.slurm.%j.err

# activate the python virtual environment with the packages we need
source ../py3env/bin/activate

# in this sbatch script, its not necessarry to create the directory,
# we already told sbatch to create it for the log files.
# mkdir -p 30_results_filter_contigs

# In this solution, our script takes the assembly and the files
# 'assembly_virus_summary.tsv' from geNomad and 
# 'quality_summary.tsv' from CheckV as an input. Set them as variables
# for readability
assembly='../1.2_assembly/10_assembly_flye/assembly.fasta'
genomadresults='10_genomad/assembly_summary/assembly_virus_summary.tsv'
checkvresults='20_checkv/quality_summary.tsv'
outdir='30_filter_contigs'

# run our script for filtering contigs based on the output of geNomad
# and CheckV as well as the output path. 
python ../python_scripts/1.3_filter_contigs.py $assembly $genomadresults $checkvresults $outdir/assembly.fasta

# deactivate the environment
deactivate

Compare the results

CheckV and geNomad follow slightly different approaches, and we expect their outputs not to match perfectly. How different are their predictions?

At the selected cutoffs, how many contigs get chosen by geNomad, how many by CheckV?

How many geNomad hits were annotated by CheckV as having high-quality?

To find how many geNomad hits were annotated by CheckV as high quality:

tail -n +2 10_genomad/assembly_summary/assembly_virus_summary.tsv | cut -f1 | cut -d'|' -f1 | sort -u > phage_contigs_genomad.list

grep -w -Ff phage_contigs_genomad.list  20_results_assessment_checkv/cross_assembly/quality_summary.tsv | grep "High-" | wc -l

Key Points

Filtering contigs by completeness and contamination is crucial to obtain an informative dataset

Tools like geNomad classify your contigs, enabling you to understand your samples

No wet-lab or dry-lab technique is perfect. Filtering non-viral contigs from your data improves its quality, helping you obtain better results

previous episode

Viromics2026

next episode

Identifying Viral Contigs II

Overview

Identify viral contigs

Challenge

Exercise

sbatch script for geNomad

Questions

bash command for getting the number of phage contigs

Estimating genome completeness

sbatch script for running checkV

Go through the CheckV results

Filter contigs

Exercise - select high-quality phage contigs

python script for selecting viral contigs

sbatch script for submitting the python script

Compare the results

To find how many geNomad hits were annotated by CheckV as high quality:

Key Points

previous episode

next episode