安装软件

安装主程序及依赖

:::default
- VirSorter2 (version >=2.2.3)
- CheckV (version >=0.7.0)
- DRAMv (version >=1.2.0)
  :::
1
2
conda create -n virome virsorter=2 checkv dram
conda activate virome

下载数据库

# vs2 db: db-vs2 ~ 10 min
virsorter setup -d /new_data/hualin/db/db-vs2 -j 50

# checkv db: checkv-db-v1.0  < 5 mins
checkv download_database /new_data/hualin/db/checkv

# DRAMv: db-dramv ~5h and ~60GB of memory
DRAM-setup.py prepare_databases --skip_uniref --output_dir db-dramv

预测

获取测试数据

1	wget -O test.fa https://bitbucket.org/MAVERICLab/virsorter2/raw/15a21fa9c1ee1d2af52b0148b167292e549d356e/test/test-for-sop.fa

运行VirSorter2

+++info 命令解析

首先以0.5 的score阈值运行VirSorter2 以保证最大的灵敏度。
只关注噬菌体 (dsDNA and ssDNA phage)，可选dsDNAphage,NCLDV,RNA,ssDNA,lavidaviridae。
序列最小长度为5000 bp，后续病毒分类的最低要求如此。
可根据自己的CPU核心数自行调整 “-j”。
“–keep-original-seq” 保留了环状和接近完整的病毒contigs (score >0.8 as a whole sequence)，后续将通过checkV修剪其尾端的潜在宿主基因并处理重复的环状contigs片段。
+++

:::default
Time: 31m7.310s with a real dataset of 90.52 MB and 207,544 sequences
:::

1	virsorter run --keep-original-seq -w vs2-pass1 -i test.fa --include-groups dsDNAphage,ssDNA --min-length 5000 --min-score 0.5 -j 50 all

+++primary 参数解析
-w 指定输出目录
-i 指定输入序列
--min-length 过滤短序列
--min-score 分数阈值
--keep-original-seq 保留环状和接近完整的病毒contigs
--include-groups 指定包含的病毒类型，用“,”分隔。可选：dsDNAphage,NCLDV,RNA,ssDNA,lavidaviridae
-j 线程数
all 直接写上就可以
+++

+++success 结果解析

final-viral-combined.fa: 病毒序列
- 鉴定为病毒的完整序列（标识为后缀 ||full）；
- 鉴定为病毒的部分序列（用后缀 ||{i}_partial 标识）； {i} 可以是从 0 到该Contig中发现的最大病毒片段数的数字；
- 带有标志基因hallmark gene的短序列（少于两个基因）被鉴定为病毒（用后缀 ||lt2gene 标识）；
final-viral-score.tsv: 每个病毒序列的评分跨组和一些关键特征，这可以用于进一步过滤
- sequence name
- score of each viral sequences across groups (多列)
- max score across groups
- max score group
- contig length
- hallmark gene count
- viral gene %
- nonviral gene %

seqname	dsDNAphage	NCLDV	RNA	ssDNA	lavidaviridae	max_score	max_score_group	length	hallmark	viral	cellular
NODE_5_length_17317_cov_8.385876\|\|full	0.993	0.847	0.005	0.060	0.467	0.993	dsDNAphage	17315	2	64.700	5.900
NODE_6_length_16611_cov_115.615064\|\|full	0.920	0.207	0.035	0.087	0.053	0.920	dsDNAphage	16610	0	3.200	0.000
NODE_8_length_14848_cov_778.417157\|\|full	1.000	0.220	0.105	0.380	0.627	1.000	dsDNAphage	14848	15	100.000	0.000
NODE_16_length_12563_cov_14.331948\|\|full	0.973	0.200	0.165	0.273	0.227	0.973	dsDNAphage	12083	0	8.000	0.000
NODE_17_length_11885_cov_350.043956\|\|full	0.653	0.513	0.050	0.080	0.047	0.653	dsDNAphage	11885	0	9.100	0.000
NODE_21_length_11527_cov_216.405073\|\|full	0.620	0.407	0.000	0.013	0.060	0.620	dsDNAphage	11526	0	10.500	5.300
NODE_23_length_11316_cov_8.144303\|\|full	0.367	0.540	0.010	0.000	0.400	0.540	NCLDV	11313	1	23.100	7.700

:::warning
不同病毒类群的分类器并非相互排斥，它们的目标病毒序列空间target viral sequence space可能存在重叠，这意味着该信息不应被使用或当作可靠的分类reliable taxonomic classification。VirSorter2的用途仅限于病毒鉴定。
:::

final-viral-boundary.tsv: 带有边界信息的表 (与其他两个文件相比，可能有额外的记录，应该忽略)。

only some of the columns in this file might be useful:
- seqname: original sequence name
- trim_orf_index_start, trim_orf_index_end: start and end ORF index on orignal sequence of identified viral sequence
- trim_bp_start, trim_bp_end: start and end position on orignal sequence of identified viral sequence
- trim_pr: score of final trimmed viral sequence
- partial: full sequence as viral or partial sequence as viral; this is defined when a full sequence has score > score cutoff, it is full (0), or else any viral sequence extracted within it is partial (1)
- pr_full: score of the original sequence
- hallmark_cnt: hallmark gene count
- group: the classifier of viral group that gives high score; this should NOT be used as reliable classification

seqname	trim_orf_index_start	trim_orf_index_end	trim_bp_start	trim_bp_end	trim_pr	trim_pr_max	prox_orf_index_start	prox_orf_index_end	prox_bp_start	prox_bp_end	prox_pr	prox_pr_max	partial	full_orf_index_start	full_orf_index_end	full_bp_start	full_bp_end	pr_full	arc	bac	euk	vir	mix	unaligned	hallmark_cnt	group	shape	seqname_new
NODE_999_length_4026_cov_7.610929	1	12	1	4025	0.547	0.547	1	12	1	4025	nan	nan	0	1	12	1	4025	0.547	0.0	0.0	0.0	0.0	0.0	100.0	0	dsDNAphage	linear	NODE_999_length_4026_cov_7.610929\|\|full
NODE_9999_length_1276_cov_11.598690	1	3	3	1274	0.955	0.955	1	3	3	1274	nan	nan	0	1	3	3	1274	0.955	0.0	0.0	0.0	0.0	0.0	100.0	0	RNA	linear	NODE_9999_length_1276_cov_11.598690\|\|full
NODE_99999_length_314_cov_4.000000	1	2	1	287	0.57	0.57	1	2	1	287	nan	nan	0	1	2	1	287	0.57	0.0	0.0	0.0	0.0	50.0	50.0	0	RNA	linear	NODE_99999_length_314_cov_4.000000\|\|full
NODE_99992_length_314_cov_4.389961	1	2	1	275	0.747	0.747	1	2	1	275	nan	nan	0	1	2	1	275	0.747	0.0	0.0	0.0	0.0	0.0	100.0	0	ssDNA	linear	NODE_99992_length_314_cov_4.389961\|\|full
NODE_9997_length_1276_cov_44.113841	1	2	1	1240	0.98	0.98	1	2	1	1240	nan	nan	0	1	2	1	1240	0.98	0.0	0.0	0.0	0.0	50.0	50.0	0	RNA	linear	NODE_9997_length_1276_cov_44.113841\|\|full

:::warning
在原病毒provirus 提取过程中，为了获得更好的敏感性，VirSorter2有时会高估病毒序列的大小。建议清除这些前病毒预测，以去除预测病毒区域边缘的潜在宿主基因，例如使用CheckV等工具。
:::

:::primary
How to pick a score cutoff?
:::

一般来说，score>0.9的人为高置信度。得分在 0.5 到 0.9 之间的可能是病毒和非病毒的混合体。很难确定区分病毒和非病毒的最佳分数，因为它取决于宿主序列和未知序列的百分比。因此，[建议使用默认截止值 (0.5) 以获得最大灵敏度，然后使用 checkV 应用质量检查步骤以消除误报（预测完整性除外）]{.rainbow}。==请继续下面的流程。==
+++

运行checkV

+++info 命令解析
Score阈值设为0.5时，VirSorter2结果中可能存在一些非病毒序列或区域。因此，使用CheckV对VirSorter2的结果进行质量控制，并修剪在原病毒proviruses末端留下的潜在宿主区域。可以根据CPU内核的可用性调整-t选项。
+++

:::default
Time: 0m28.795s
:::
1
checkv end_to_end vs2-pass1/final-viral-combined.fa checkv -t 50 -d /new_data/hualin/db/checkv/checkv-db-v1.0
+++success 结果解析
- ./checkv/
- completeness.tsv
- complete_genomes.tsv
- contamination.tsv
- proviruses.fna
- quality_summary.tsv
- viruses.fna
- tmp
  +++

再次运行VirSorter2

+++info 命令解析

再次利用checkV-trimmed序列运行VirSorter2以得到”affi-contigs.tab”文件，该文件将作为DRAMv的输入以鉴定AMGs。
注意”--seqname-suffix-off”选项保留了原始的输入序列名称，因为我们确信在本步骤中，不可能从同一条contig中获得>1个原病毒。
“--viral-gene-enrich-off”选项关闭了病毒基因要多于宿主基因的要求，以确保VirSorter2在这一步不做任何筛查。
以上两个选项需要VirSorter2版本>=2.2.1。
可选所有病毒：dsDNAphage,NCLDV,RNA,ssDNA,lavidaviridae
+++

:::default
Time: 18m30.896s
:::

1
2
3

cat checkv/proviruses.fna checkv/viruses.fna > checkv/combined.fna

virsorter run --seqname-suffix-off --viral-gene-enrich-off --provirus-off --prep-for-dramv -i checkv/combined.fna -w vs2-pass2 --include-groups dsDNAphage,ssDNA --min-length 5000 --min-score 0.5 -j 50 all

+++success 结果解析

./vs2-pass2/
inal-viral-combined.fa
final-viral-score.tsv
for-dramv/final-viral-combined-for-dramv.fa
for-dramv/viral-affi-contigs-for-dramv.tab
+++

运行DRAMv

+++info 命令解析
使用DRAMv注释鉴定的病毒序列，以用于后续人工整理manual curation。可通过”--threads”控制调用的CPU核心数。
+++

:::default
Time: 8.81h
:::

# step 1 annotate 耗时步骤，建议投后台运行
DRAM-v.py annotate -i vs2-pass2/for-dramv/final-viral-combined-for-dramv.fa -v vs2-pass2/for-dramv/viral-affi-contigs-for-dramv.tab -o dramv-annotate --skip_trnascan --threads 50 --min_contig_size 1000

#step 2 summarize anntotations
DRAM-v.py distill -i dramv-annotate/annotations.tsv -o dramv-distill

+++success 结果解析

dramv-annotate/annotations.tsv
dramv-annotate/genbank 各条contig的gbk文件
dramv-annotate/genes.faa
dramv-annotate/genes.fna
dramv-annotate/genes.gff
dramv-annotate/rrnas.tsv
dramv-annotate/scaffolds.fna
dramv-annotate/vMAGs 各条contig的fasta文件
+++

过滤

依据病毒和宿主基因数、分数、hallmark基因数以及contig长度进行筛选

来自checkV的病毒和宿主基因计数可用于假阳性筛查。由于checkV在预测病毒基因方面非常保守，那些由checkV预测的病毒基因的序列应该是病毒的，而那些没有被checkV预测到病毒基因的序列更可能是非病毒的。

基于我们对土壤宏基因组的基准测试，(1) 那些没有预测到病毒和宿主基因的序列是病毒；(2) 没有病毒基因但有2个或2个以上宿主基因的大多数为非病毒基因；(3) 那些没有病毒基因和具有1个宿主基因的很难区分其为病毒还是非病毒（可能是可移动的基因元件，类似于VirSorter1中的第3类），除非手动检查，否则应该丢弃。

只选择那些大于10kb的用于手动检查，因为太短的无法分辨。还有那些VirSorter2 得分≥0.95或hallmark基因计数>2的大多数是病毒。==这些经验筛选标准总结如下：==

:::success
Keep1: viral_gene >0
:::
:::success
Keep2: viral_gene =0 AND (host_gene =0 OR score >=0.95 OR hallmark >2)
:::
:::warning
Manual check: (NOT in Keep1 OR Keep2) AND viral_gene =0 AND host_gene =1 AND length >=10kb
:::
:::danger
Discard: the rest
:::

[要查看病毒基因、宿主基因、评分和序列的特征标记hallmark of sequences，您可以合并”vs2-pass1/final-viral-score.tsv”和”checkv/contamination.tsv”，并在电子表格中过滤。]{.aqua}本尊为各位提供了Perl脚本cat_tsv.pl以实现机动合并！直接在终端运行perl cat_tsv.pl即可得到合并后的文件forCheck.txt。

#!/usr/bin/perl
use strict;
use warnings;

my %hash;
my $head_checkv;
my $head_pass1;
my $count = 0;
my $num = 0;
open IN, "checkv/contamination.tsv" || die;
while (<IN>) {
	chomp;
	$_ =~s/[\r\n]+//g;
	$count++;
	if ($count == 1) {
		$head_checkv = $_;
	}else {
		my @lines = split /\t/;
		$hash{$lines[0]} = $_;
	}
}
close IN;

open IN, "vs2-pass1/final-viral-score.tsv" || die;
open OUT, ">forCheck.txt" || die;
while (<IN>) {
	chomp;
	$_ =~s/[\r\n]+//g;
	$num++;
	if ($num == 1) {
		$head_pass1 = $_;
		print OUT "$head_pass1\t$head_checkv\n";
	}else {
		my @line = split /\t/;
		if (exists $hash{$line[0]}) {
			print OUT "$_\t$hash{$line[0]}\n";
		}else {
			print OUT "$_\n";
		}
	}
}
close IN;
close OUT;

接下来按照上面的四条规则对forCheck.txt进行拆分，得到病毒Virus.txt、手动核对Manual_check.txt及抛弃的Discard.txt。可以自己看，也可以用get_virus.pl来完成。

#!/usr/bin/perl
use strict;
use warnings;

my $count = 0;
open IN, "forCheck.txt" || die;
open VIRUS, ">Virus.txt" || die;
open MANUAL, ">Manual_check.txt" || die;
open DISCARD, ">Discard.txt" || die;
while (<IN>) {
	chomp;
	$_ =~s/[\r\n]+//g;
	$count++;
	if ($count == 1) {
		print VIRUS $_ . "\n";
		print MANUAL $_ . "\n";
		print DISCARD $_ . "\n";
	}else {
		my @lines = split /\t/;
		if ($lines[15] > 0) { # virus keep1
			print VIRUS $_ . "\n";
		}elsif (($lines[15] == 0) && ($lines[16] == 0)) { # virus keep2
			print VIRUS $_ . "\n";
		}elsif (($lines[15] == 0) && ($lines[6] >= 0.95)) { # virus keep2
			print VIRUS $_ . "\n";
		}elsif (($lines[15] == 0) && ($lines[9] > 2)) { # virus keep2
			print VIRUS $_ . "\n";
		}elsif (($lines[15] == 0) && ($lines[16] == 1) && ($lines[8] >= 10000)) { # manual check
			print MANUAL $_ . "\n";
		}else { # discard
			print DISCARD $_ . "\n";
		}
	}
}
close IN;
close VIRUS;
close MANUAL;
close DISCARD;

用法：脚本与forCheck.txt放在同一目录，终端运行如下命令。

1	perl get_virus.pl

对于Manual_check.txt中的序列，需要用眼睛和脑袋去dramv-annotate/annotations.tsv中找注释信息，然后根据下面的方法判断其属于病毒还是细胞，如果是病毒，就将其所在的那行信息复制到Virus.txt文件的末尾，并保存。

依据DRAMv注释筛选

在病毒和宿主中都有一些共同的基因(如{脂多糖^LPS}相关)和移动元件，这些基因可能导致上述“Keep2”类别中的假阳性。因此，++要谨慎对待带有这些基因的contigs++{.dot}。专家们已经编制了一份与此相关的[“可疑”]{.red}基因列表。我们可以使用“Keep2”类别中的contigs对DRAMv表取子集，并在DRAMv子集表中筛选“可疑”基因(忽略大小写，例如在“grep”中使用“-i”选项)，然后将带有这些基因的contigs放入“手动检查”类别中。

手动检查{.tabset}

对于存在于“手动检查”manual check类别中的序列，可以观察其在”dramv-annotate/annotations.tsv”中的注释信息。本步骤很难标准化，下面是一些经验之谈：

判定contig为病毒的标准
判定contig为非病毒的标准

;;;id1 判定contig为病毒的标准

结构基因、hallmark基因、注释缺失或假设性富集depletion in annotations or enrichment for hypotheticals (~10%的基因具有non-hypothetical注释)。
缺乏hallmarks，但>=50%已注释的基因为病毒，且其中至少一半以上的viral bitcore >100，且contig的长度<50kb。
Provirus: 整合酶/重组酶/切除酶/阻遏子Integrase/recombinase/excisionase/repressor，在一侧富集了病毒基因。
Provirus: 基因组中存在“break”：两个基因之间的间隙gap对应于一个链开关strand switch，更高的编码密度，注释缺失，以及在一侧噬菌体基因的富集。
仅有~1-3个基因有注释，但至少一半命中病毒，且命中基因的bitscore不超过病毒bitscore的150% ，Few annotations only ~1-3 genes, but with at least half hitting to viruses, and where the genes hitting cells have a bitscore no more than 150% that of the viral bitscores 且/或病毒的bitscore >100 and/or viral bitscores are >100。
LPS (脂多糖) 外观区域对病毒基因的命中率也非常高looking regions if also has very strong hits to viral genes，bitscore>100。
;;;

;;;id1 判定contig为非病毒的标准

细胞样基因cellular like genes是病毒基因的3倍，几乎所有基因都有注释，没有基因只命中病毒，也没有病毒标志基因hallmark genes。
缺乏任何病毒hallmark genes，且长度 >50kb。
许多明显的细胞基因字符串，没有其他病毒标志基因。在基准测试中遇到的例子包括 1) CRISPR Cas, 2) ABC transporters, 3) Sporulation proteins, 4) Two-component systems, 5) Secretion system。这其中一些可能是由病毒编码的，但在没有进一步证据的情况下并不表明是病毒contig。
多个质粒基因或转座酶，但没有明确的只命中病毒的基因。
注释信息很少，仅有~1-3个基因同时命中了病毒和细胞基因，但有stronger bitscores支持其为细胞基因。
没有强有力的命中病毒的脂多糖样区域LPS looking regions if no strong viral hits。富含通常与脂多糖相关的基因，如{外聚酶^epimerases}、{糖基转移酶^glycosyl transferases}、{酰基转移酶^acyltransferase}、{短链脱氢酶/还原酶^dehydrogenase/reductase short-chain}、{脱水酶^dehydratase}。
注释为 Type IV 和/或 Type VI 分析系统，并被非病毒基因围绕。
注释信息很少，仅有~1-3个基因全部命中细胞基因 (即使 bitscore <100) ，且没有命中病毒的基因。
;;;

:::warning
最后，用户要注意，VirSorter 2和/或checkV预测的任何原病毒边界只是一个近似的估计(寻找“ends”在前噬菌体发现中是一个相当具有挑战性的问题)，也需要仔细地手工检查，特别是在AMG研究中。
:::

最终我们需要拿到病毒contig的序列，用get_virus_seqs.pl完成。

#!/usr/bin/perl
use strict;
use warnings;

my %hash;
open IN, "checkv/combined.fna" || die;
local $/ = ">";
<IN>;
while (<IN>) {
	chomp;
	my ($header, $seq) = split(/\n/, $_, 2);
	my $id;
	if ($header =~/(\S+)\|\|.+/) {
		$id = $1;
	}else {
		$id = $header;
	}
	$hash{$id} = $seq;
}
close IN;

local $/ = "\n";
open IN, "Virus.txt" || die;
open OUT, ">Virus.fas" || die;
open NO, ">NoSeqs.ids" || die;
<IN>;

while (<IN>) {
	chomp;
	my @lines = split /\t/;
	$lines[0] =~/(\S+)\|\|/;
	if (exists $hash{$1}) {
		print OUT ">$lines[0]\n$hash{$1}\n";
	}else {
		print NO "$lines[0]\n";
	}
}
close IN;
close OUT;
close NO;

用法：终端运行如下命令即可得到序列文件Virus.fas。

1	perl get_virus_seqs.pl

脚本获取

关注公众号“生信之巅”，聊天窗口回复“1551”获取下载链接。

敬告：使用文中脚本请引用本文网址，请尊重本人的劳动成果，谢谢！Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

参考

graph TD A [Reads] --> B((MetaSpaDES)) B --> C [Contigs] C --> D ((VirSorter)) C --> E ((VirFinder)) E --> F {Score >= 0.7 & p < 0.05} F -- YES --> G [Virus] F -- NO --> H [Contigs] D --> G H --> I ((CAT)) I --> J {< 40%} J -- YES --> G J -- NO --> K [Not Virus] G --> L ((Nucmer)) L --> M {Identity >= 95% & Coverage >= 80%} M -- YES --> N [Grouped] M -- NO --> O [Not Grouped] G --> P ((Prodigal)) P --> Q [Amino Acid Sequences] Q --> R {Length < 10 kb} R -- YES --> S ((Blastp)) R -- NO --> T ((vConTACT2)) S --> U [Annotated Genes] T --> U

了尘兰若的小坑

分享微生物生物信息学分析方法，欢迎加入QQ群交流945751012，不接受群内广告！

安装软件

预测

过滤

依据病毒和宿主基因数、分数、hallmark基因数以及contig长度进行筛选

依据DRAMv注释筛选

手动检查{.tabset}

脚本获取

参考

你的赏识是我前进的动力