Detection, annotation, and taxonomic classification of viral contigs in metagenomic and metatranscriptomic assemblies

public public 1yr ago Version: Version 1 0 bookmarks
View Workflow
Report
virify
virify
virify
virify
virify
virify
View Workflow

VIRify

Sankey plot

VIRify is a recently developed pipeline for the detection, annotation, and taxonomic classification of viral contigs in metagenomic and metatranscriptomic assemblies. The pipeline is part of the repertoire of analysis services offered by MGnify . VIRify’s taxonomic classification relies on the detection of taxon-specific profile hidden Markov models (HMMs), built upon a set of 22,014 orthologous protein domains and referred to as ViPhOGs.

VIRify was implemented in CWL.

What do I need?

The current implementation uses CWL version 1.2 dev+2. It was tested using Toil version 4.10 as the workflow engine and conda to manage the software dependencies.

Docker - Singularity support

Soon…

Setup environment 

conda env create -f cwl/requirements/conda_env.yml
conda activate viral_pipeline

Basic execution 

cd cwl/
virify.sh -h

A note about metatranscriptomes

Although VIRify has been benchmarked and validated with metagenomic data in mind, it is also possible to use this tool to detect RNA viruses in metatranscriptome assemblies (e.g. SARS-CoV-2). However, some additional considerations for this purpose are outlined below:
1. Quality control: As for metagenomic data, a thorough quality control of the FASTQ sequence reads to remove low-quality bases, adapters and host contamination (if appropriate) is required prior to assembly. This is especially important for metatranscriptomes as small errors can further decrease the quality and contiguity of the assembly obtained. We have used TrimGalore for this purpose.

2. Assembly: There are many assemblers available that are appropriate for either metagenomic or single-species transcriptomic data. However, to our knowledge, there is no assembler currently available specifically for metatranscriptomic data. From our preliminary investigations, we have found that transcriptome-specific assemblers (e.g. rnaSPAdes ) generate more contiguous and complete metatranscriptome assemblies compared to metagenomic alternatives (e.g. MEGAHIT and metaSPAdes ).

3. Post-processing: Metatranscriptomes generate highly fragmented assemblies. Therefore, filtering contigs based on a set minimum length has a substantial impact in the number of contigs processed in VIRify. It has also been observed that the number of false-positive detections of VirFinder (one of the tools included in VIRify) is lower among larger contigs. The choice of a length threshold will depend on the complexity of the sample and the sequencing technology used, but in our experience any contigs <2 kb should be analysed with caution.

4. Classification: The classification module of VIRify depends on the presence of a minimum number and proportion of phylogenetically-informative genes within each contig in order to confidently assign a taxonomic lineage. Therefore, short contigs typically obtained from metatranscriptome assemblies remain generally unclassified. For targeted classification of RNA viruses (for instance, to search for Coronavirus-related sequences), alternative DNA- or protein-based classification methods can be used. Two of the possible options are: (i) using MashMap to screen the VIRify contigs against a database of RNA viruses (e.g. Coronaviridae) or (ii) using hmmsearch to screen the proteins obtained in the VIRify contigs against marker genes of the taxon of interest.

Contact us

MGnify helpdesk

Code Snippets

13
14
15
"""
viral_contigs_annotation.py -o . -p ${faa} -t ${tab}
"""
NextFlow From line 13 of modules/annotation.nf 
15
16
17
"""
contig_taxonomic_assign.py -i ${tab} -d ${db} --factor ${factor} --taxthres ${params.taxthres}
"""
NextFlow From line 15 of modules/assign.nf 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
"""
if [[ "${set_name}" == "all" ]]; then
  # concat all taxonomy results per assembly regardless of HC, LC or PP
  cat *_annotation_taxonomy.tsv | grep -v contig_ID > tmp.tsv
  NAME=all
else
  cp ${tbl} tmp.tsv
  NAME=${set_name}
fi

# genus
grep -v contig_ID tmp.tsv | awk -v SAMPLE="${name}" 'BEGIN{FS="\\t"};{if(\$2!="" && \$2 !~ /^0/){print SAMPLE"\\tgenus\\t"\$2}}' | sort | uniq -c | awk '{printf \$2"\\t"\$3"\\t"\$4"\\t"\$1"\\n"}' > \$NAME"_summary.tsv"

# subfamily
grep -v contig_ID tmp.tsv | awk -v SAMPLE="${name}" 'BEGIN{FS="\\t"};{if(\$3!="" && \$3 !~ /^0/){print SAMPLE"\\tsubfamily\\t"\$3}}' | sort | uniq -c | awk '{printf \$2"\\t"\$3"\\t"\$4"\\t"\$1"\\n"}' >> \$NAME"_summary.tsv"

# family
grep -v contig_ID tmp.tsv | awk -v SAMPLE="${name}" 'BEGIN{FS="\\t"};{if(\$4!="" && \$4 !~ /^0/){print SAMPLE"\\tfamily\\t"\$4}}' | sort | uniq -c | awk '{printf \$2"\\t"\$3"\\t"\$4"\\t"\$1"\\n"}' >> \$NAME"_summary.tsv"

# order
grep -v contig_ID tmp.tsv | awk -v SAMPLE="${name}" 'BEGIN{FS="\\t"};{if(\$5!="" && \$5 !~ /^0/){print SAMPLE"\\torder\\t"\$5}}' | sort | uniq -c | awk '{printf \$2"\\t"\$3"\\t"\$4"\\t"\$1"\\n"}' >> \$NAME"_summary.tsv"

if [ -s \$NAME"_summary.tsv" ]; then
  balloon.R "\${NAME}_summary.tsv" "\${NAME}_balloon.svg" 10 8 
fi
"""
NextFlow From line 12 of modules/balloon.nf 
18
19
20
"""
  imgvr_merge.py -f ${blast_filtered} -d IMG_VR_2018-07-01_4/IMGVR_all_Sequence_information.tsv -o \$(basename ${blast_filtered} .blast).meta
"""
NextFlow From line 18 of modules/blast_filter.nf 
22
23
24
"""
  imgvr_merge.py -f ${blast_filtered} -d ${db}/IMG_VR_2018-07-01_4/IMGVR_all_Sequence_information.tsv -o \$(basename ${blast_filtered} .blast).meta
"""
NextFlow From line 22 of modules/blast_filter.nf 
18
19
20
21
22
23
24
25
"""
  sed -i "s/ /|/" ${fasta}
  HEADER_BLAST="qseqid\\tsseqid\\tpident\\tlength\\tmismatch\\tgapopen\\tqstart\\tqend\\tqlen\\tsstart\\tsend\\tevalue\\tbitscore\\tslen"
  printf "\$HEADER_BLAST\\n" > ${confidence_set_name}.blast

  blastn -task blastn -num_threads ${task.cpus} -query ${fasta} -db IMG_VR_2018-07-01_4/IMGVR_all_nucleotides.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend qlen sstart send evalue bitscore slen" >> ${confidence_set_name}.blast
  awk '{if(\$4>0.8*\$9){print \$0}}' ${confidence_set_name}.blast >> ${confidence_set_name}.filtered.blast
"""
NextFlow BLAST From line 18 of modules/blast.nf 
27
28
29
30
31
32
33
34
"""
  sed -i "s/ /|/" ${fasta}
  HEADER_BLAST="qseqid\\tsseqid\\tpident\\tlength\\tmismatch\\tgapopen\\tqstart\\tqend\\tqlen\\tsstart\\tsend\\tevalue\\tbitscore\\tslen"
  printf "\$HEADER_BLAST\\n" > ${confidence_set_name}.blast

  blastn -task blastn -num_threads ${task.cpus} -query ${fasta} -db ${db}/IMG_VR_2018-07-01_4/IMGVR_all_nucleotides.fna -evalue 1e-10 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend qlen sstart send evalue bitscore slen" >> ${confidence_set_name}.blast
  awk '{if(\$4>0.8*\$9){print \$0}}' ${confidence_set_name}.blast >> ${confidence_set_name}.filtered.blast
"""
NextFlow BLAST From line 27 of modules/blast.nf 
14
15
16
17
        """
	checkv end_to_end ${fasta} -d ${database} -t ${task.cpus} ${confidence_set_name}
        cp ${confidence_set_name}/quality_summary.tsv ${confidence_set_name}_quality_summary.tsv 
        """
NextFlow From line 14 of modules/checkV.nf 
19
20
21
22
"""
checkv end_to_end ${fasta} -d ${database} -t ${task.cpus} ${confidence_set_name} 
cp ${confidence_set_name}/quality_summary.tsv ${confidence_set_name}_quality_summary.tsv 
"""
NextFlow From line 19 of modules/checkV.nf 
25
26
27
28
29
"""
mkdir negative_result_${confidence_set_name}.tsv
echo "contig_id	contig_length	genome_copies	gene_count	viral_genes	host_genes	checkv_quality	miuvig_quality	completeness	completeness_method	contamination	provirus	termini	warnings" > ${confidence_set_name}_quality_summary.tsv
echo "pos_phage_0	146647	1.0	243	141	1	High-quality	High-quality	97.03	AAI-based	0.0	No" >> ${confidence_set_name}_quality_summary.tsv   
"""
NextFlow From line 25 of modules/checkV.nf 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
"""
# combine
if [[ ${set_name} == "all" ]]; then
  cat *.tsv | grep -v Abs_Evalue_exp | sort | uniq > ${id}.tsv
  cat *.fasta > ${id}.fasta
else
  cp ${annotation_table} ${id}.tsv
  cp ${assembly} ${id}.fasta
fi

# now we remove contigs shorter 1500 kb and very long ones because ChromoMap has an error when plotting to distinct lenghts
# we also split into multiple files when we have many contigs, chunk size default 30
filter_for_chromomap.rb ${id}.fasta ${id}.tsv ${id}.contigs.txt ${id}.anno.txt 1500
"""
NextFlow From line 14 of modules/chromomap.nf 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
"""
#!/usr/bin/env Rscript

library(chromoMap)
library(ggplot2)
library(plotly)

contigs <- list()
annos <- list()
contigs <- dir(pattern = "*.contigs.txt")
annos <- dir(pattern = "*.anno.txt")

for (k in 1:length(contigs)){
  c = contigs[k]
  a = annos[k]

  # check if a file is empty
  if (file.info(c)\$size == 0 || file.info(a)\$size == 0) {
    next
  }

  # check how many categories we have
  categories <- c("limegreen", "orange","grey")
  df <- read.table(a, sep = "\\t")
  set <- unique(df\$V5)
  if ( length(set) == 2 ) {
    if ( set[1] == 'High confidence' && set[2] == 'Low confidence') {
      categories <- c("limegreen", "orange")
    }
    if ( set[1] == 'High confidence' && set[2] == 'No hit') {
      categories <- c("limegreen", "grey")
    }
    if ( set[1] == 'Low confidence' && set[2] == 'No hit') {
      categories <- c("orange", "grey")
    }
  }
  if ( length(set) == 1 ) {
    if ( set[1] == 'High confidence') {
      categories <- c("limegreen")
    }
    if ( set[1] == 'Low confidence') {
      categories <- c("orange")
    }
    if ( set[1] == 'No hit') {
      categories <- c("grey")
    }
  }

  p <-  chromoMap(c, a,
    data_based_color_map = T,
    data_type = "categorical",
    data_colors = list(categories),
    legend = T, lg_y = 400, lg_x = 100, segment_annotation = T,
    left_margin = 100, canvas_width = 1000, chr_length = 8, ch_gap = 6)
  htmlwidgets::saveWidget(as_widget(p), paste("${id}.chromomap-", k, ".html", sep=''))
}    
"""
NextFlow ggplot2 plotly chromoMap From line 44 of modules/chromomap.nf 
 9
10
11
"""
fastp -i ${reads[0]} -I ${reads[1]} --thread ${task.cpus} -o ${name}.R1.fastp.fastq.gz -O ${name}.R2.fastp.fastq.gz
"""
NextFlow fastp From line 9 of modules/fastp.nf 
 8
 9
10
11
"""
mkdir fastqc
fastqc -t ${task.cpus} -o fastqc *.fastq.gz
"""
NextFlow FastQC From line 8 of modules/fastqc.nf 
12
13
14
"""
hmmscan_format_table.py -t ${hmmer_tbl} -o ${set_name}_modified
"""
NextFlow From line 12 of modules/hmm_postprocessing.nf 
13
14
15
16
17
18
19
20
21
22
23
24
25
"""
if [[ ${params.db} == "viphogs" ]]; then
  if [[ ${params.version} == "v1" ]]; then
    hmmscan --cpu ${task.cpus} --noali -E "0.001" --domtblout ${set_name}_${params.db}_hmmscan.tbl ${db}/${db}.hmm ${faa}
  else
    hmmscan --cpu ${task.cpus} --noali --cut_ga --domtblout ${set_name}_${params.db}_hmmscan_cutga.tbl ${db}/${db}.hmm ${faa}
    #filter evalue for models that dont have any GA cutoff
    awk '{if(\$1 ~ /^#/){print \$0}else{if(\$7<0.001){print \$0}}}' ${set_name}_${params.db}_hmmscan_cutga.tbl > ${set_name}_${params.db}_hmmscan.tbl
  fi
else
  hmmscan --cpu ${task.cpus} --noali -E "0.001" --domtblout ${set_name}_${params.db}_hmmscan.tbl ${db}/${db}.hmm ${faa}
fi
"""
NextFlow hmmscan (EBI) From line 13 of modules/hmmscan.nf 
13
14
15
16
17
18
19
20
21
22
"""
if [[ "${set_name}" == "all" ]]; then
  grep contig_ID *.tsv | awk 'BEGIN{FS=":"};{print \$2}' | uniq > ${name}.tmp
  grep -v "contig_ID" *.tsv | awk 'BEGIN{FS=":"};{print \$2}' | uniq >> ${name}.tmp
  cp ${name}.tmp ${name}.tsv
  generate_counts_table.py -f ${name}.tsv -o ${name}.krona.tsv
else
  generate_counts_table.py -f ${tbl} -o ${set_name}.krona.tsv
fi
"""
NextFlow From line 13 of modules/krona.nf 
34
35
36
37
38
39
40
"""
if [[ ${set_name} == "all" ]]; then
  ktImportText -o ${name}.krona.html ${krona_file}
else
  ktImportText -o ${set_name}.krona.html ${krona_file}
fi
"""
NextFlow Krona From line 34 of modules/krona.nf 
15
16
17
18
"""
sed -i "s/ /|/" ${fasta}
mashmap -q ${fasta} -r ${reference} -t ${task.cpus} -o ${confidence_set_name}_mashmap_hits.tsv --noSplit -s ${params.mashmap_len}
"""
NextFlow MashMap From line 15 of modules/mashmap.nf 
 9
10
11
"""
multiqc -i ${name} .
"""
NextFlow MultiQC From line 9 of modules/multiqc.nf 
18
19
20
21
"""
touch virsorter_metadata.tsv
parse_viral_pred.py -a ${fasta} -f ${virfinder} -p ${pprmeta} -s ${virsorter}/Predicted_viral_sequences/*.fasta &> ${name}_virus_predictions.log
"""
NextFlow From line 18 of modules/parse.nf 
14
15
16
17
18
19
20
21
22
23
24
25
26
    """
  	# get only contig IDs that have at least one annotation hit
    cat ${tab} | sed 's/|/VIRIFY/g' > virify.tmp 
	  IDS=\$(awk 'BEGIN{FS="\\t"};{if(\$6!="No hit"){print \$1}}' virify.tmp | sort | uniq | grep -v Contig)
	  head -1 ${tab} > ${set_name}_prot_ann_table_filtered.tsv
	  for ID in \$IDS; do
		  awk -v id="\$ID" '{if(id==\$1){print \$0}}' virify.tmp >> ${set_name}_prot_ann_table_filtered.tsv
	  done
    sed -i 's/VIRIFY/|/g' ${set_name}_prot_ann_table_filtered.tsv
    mkdir -p ${set_name}_mapping_results
    cp ${set_name}_prot_ann_table_filtered.tsv ${set_name}_mapping_results/
    make_viral_contig_map.R -o ${set_name}_mapping_results -t ${set_name}_prot_ann_table_filtered.tsv
    """
NextFlow From line 14 of modules/plot_contig_map.nf 
16
17
18
19
"""
[ -d "pprmeta" ] && cp pprmeta/* .
./PPR_Meta ${fasta} ${name}_pprmeta.csv
"""
NextFlow From line 16 of modules/pprmeta.nf 
38
39
40
41
42
"""
git clone https://github.com/mult1fractal/PPR-Meta.git
mv PPR-Meta/* .  
rm -r PPR-Meta
"""
NextFlow From line 38 of modules/pprmeta.nf 
15
16
17
"""
prodigal -p "meta" -a ${confidence_set_name}_prodigal.faa -i ${fasta}
"""
NextFlow prodigal From line 15 of modules/prodigal.nf 
14
15
16
17
"""
[ -d "models" ] && cp models/* .
ratio_evalue_table.py -i ${modified_table} -t ${model_metadata} -o ${set_name}_modified_informative.tsv -e ${params.evalue}
"""
NextFlow From line 14 of modules/ratio_evalue.nf 
13
14
15
16
17
18
19
20
"""
krona_table_2_sankey_table.rb ${krona_table} ${set_name}.sankey.tsv

# select the top ${params.sankey} hits with highest count because otherwise sankey gets messy
sort -k1,1nr ${set_name}.sankey.tsv | head -${params.sankey} > ${set_name}.sankey.filtered.tsv

tsv2json.rb ${set_name}.sankey.filtered.tsv ${set_name}.sankey.filtered-${params.sankey}.json
"""
NextFlow From line 13 of modules/sankey.nf 
37
38
39
40
41
42
43
44
45
46
47
48
"""
#!/usr/bin/env Rscript

library(sankeyD3)
library(magrittr)

Taxonomy <- jsonlite::fromJSON("${json}")

# print to HTML file
sankey = sankeyNetwork(Links = Taxonomy\$links, Nodes = Taxonomy\$nodes, Source = "source", Target = "target", Value = "value", NodeID = "name", units = "count", fontSize = 22, nodeWidth = 30, nodeShadow = TRUE, nodePadding = 30, nodeStrokeWidth = 1, nodeCornerRadius = 10, dragY = TRUE, dragX = TRUE, numberFormat = ",.3g", align = "left", orderByPath = TRUE)
saveNetwork(sankey, file = '${id}.sankey.html')
"""
NextFlow magrittr From line 37 of modules/sankey.nf 
 9
10
11
12
'''
spades.py --meta --only-assembler -1 !{reads[0]} -2 !{reads[1]} -t !{task.cpus} -o assembly
mv assembly/contigs.fasta !{name}.fasta
'''
NextFlow SPAdes From line 9 of modules/spades.nf 
17
18
19
20
"""
run_virfinder.Rscript ${model} ${fasta} .
awk '{print \$1"\\t"\$2"\\t"\$3"\\t"\$4}' ${name}*.tsv > ${name}.txt
"""
NextFlow From line 17 of modules/virfinder.nf 
36
37
38
"""
wget ftp://ftp.ebi.ac.uk/pub/databases/metagenomics/viral-pipeline/virfinder/VF.modEPV_k8.rda
"""
NextFlow From line 36 of modules/virfinder.nf 
17
18
19
"""
wrapper_phage_contigs_sorter_iPlant.pl -f ${fasta} --db 2 --wdir virsorter --ncpu ${task.cpus} --data-dir ${database} --virome
"""
NextFlow virsorter From line 17 of modules/virsorter.nf 
21
22
23
"""
wrapper_phage_contigs_sorter_iPlant.pl -f ${fasta} --db 2 --wdir virsorter --ncpu ${task.cpus} --data-dir ${database}
"""
NextFlow virsorter From line 21 of modules/virsorter.nf 
17
18
19
20
21
22
23
24
25
26
"""
write_viral_gff.py \
-v ${viphos_annotations.join(' ')} \
-c ${quality_summaries.join(' ')} \
-t ${taxonomies.join(' ')} \
-s ${name} \
-a ${fasta}

gt gff3validator ${name}_virify.gff
"""
NextFlow From line 17 of modules/write_gff.nf
ShowHide 22 more snippets with no or duplicated tags.

Login to post a comment if you would like to share your experience with this workflow.

Do you know this workflow well? If so, you can request seller status , and start supporting this workflow.

Free

Created: 1yr ago
Updated: 1yr ago
Maitainers: public
URL: https://github.com/EBI-Metagenomics/emg-viral-pipeline
Name: virify
Version: Version 1
Badge:
workflow icon

Insert copied code into your website to add a link to this workflow.

Other Versions:
Version Version 1
Downloaded: 0
Copyright: Public Domain
License: Apache License 2.0
  • Future updates

Related Workflows

psychip_snakemake
Show Details View Workflow
ENCODE pipeline for histone marks developed for the psychENCODE project
public

public

psychip pipeline is an improved version of the ENCODE pipeline for histone marks developed for the psychENCODE project. The o...
raw sequence reads Alignment Sequence alignment report macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
ncov_2
Show Details View Workflow
Near-real time tracking of SARS-CoV-2 in Connecticut
public

public

Repository containing scripts to perform near-real time tracking of SARS-CoV-2 in Connecticut using genomic data. This pipeli...
JSON nextclade Augur Biopython FOCUS Pandas Snakemake bs4 epiweeks geopy matplotlib numpy pycountry pycountry-convert uszipcode
Free
28
cellranger-snakemake-gke
Show Details View Workflow
snakemake workflow to run cellranger on a given bucket using gke.
public

public

A Snakemake workflow for running cellranger on a given bucket using Google Kubernetes Engine. The usage of this workflow ...
macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
atlas
Show Details View Workflow
ATLAS - Three commands to start analyzing your metagenome data
public

public

Metagenome-atlas is a easy-to-use metagenomic pipeline based on snakemake. It handles all steps from QC, Assembly, Binning, t...
raw sequence reads Genome assembly Annotation track checkm2 gunc prodigal snakemake-wrapper-utils MEGAHIT Atlas BBMap Biopython BioRuby Bwa-mem2 cd-hit CheckM DAS Diamond eggNOG-mapper v2 MetaBAT 2 Minimap2 MMseqs MultiQC Pandas Picard pyfastx SAMtools SemiBin Snakemake SPAdes SqueezeMeta TADpole VAMB CONCOCT ete3 gtdbtk h5py networkx numpy plotly psutil utils metagenomics
Free
175
rna-seq-star-deseq2
Show Details View Workflow
RNA-seq workflow using STAR and DESeq2
public

public

This workflow performs a differential gene expression analysis with STAR and Deseq2. The usage of this workflow is described ...
raw sequence reads Gene expression profiling ID list snakemake-wrapper-utils biomaRt PCAtools BiocParallel BWA Cutadapt DESeq2 GFFutils MultiQC Pandas SAMtools Snakemake STAR cli stringr tidyverse RNA-Seq
Free
41
dna-seq-gatk-variant-calling
Show Details View Workflow
This Snakemake pipeline implements the GATK best-practices workflow
public

public

This Snakemake pipeline implements the GATK best-practices workflow for calling small germline variants. The usage of thi...
VCF raw sequence reads Variant calling genetic variants gatk rust-bio-tools snakemake-wrapper-utils tabix BCFtools BWA FastQC MultiQC Pandas Picard SAMtools Snakemake Trimmomatic Variant Effect Predictor (VEP) common matplotlib numpy seaborn DNA
Free
57