# 安装软件
安装主程序及依赖
- VirSorter2 (version >=2.2.3)
- CheckV (version >=0.7.0)
- DRAMv (version >=1.2.0)
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
# 预测
获取测试数据
wget -O test.fa https://bitbucket.org/MAVERICLab/virsorter2/raw/15a21fa9c1ee1d2af52b0148b167292e549d356e/test/test-for-sop.fa
运行 VirSorter2
命令解析
- 首先以 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 片段。
Time: 31m7.310s with a real dataset of 90.52 MB and 207,544 sequences
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
参数解析
-w 指定输出目录
- i 指定输入序列
--min-length 过滤短序列
--min-score 分数阈值
--keep-original-seq 保留环状和接近完整的病毒 contigs
--include-groups 指定包含的病毒类型,用 “,” 分隔。可选:dsDNAphage,NCLDV,RNA,ssDNA,lavidaviridae
-j 线程数
all 直接写上就可以结果解析
final-viral-combined.fa: 病毒序列
- 鉴定为病毒的完整序列(标识为后缀 ||full);
- 鉴定为病毒的部分序列(用后缀 ||{i}_partial 标识); {i} 可以是从 0 到该 Contig 中发现的最大病毒片段数的数字;
- 带有标志基因的短序列(少于两个基因)被鉴定为病毒(用后缀 ||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 不同病毒类群的分类器并非相互排斥,它们的目标病毒序列空间可能存在重叠,这意味着该信息不应被使用或当作可靠的分类。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 在原病毒提取过程中,为了获得更好的敏感性,VirSorter2 有时会高估病毒序列的大小。建议清除这些前病毒预测,以去除预测病毒区域边缘的潜在宿主基因,例如使用 CheckV 等工具。
How to pick a score cutoff?
一般来说,score>0.9 的人为高置信度。得分在 0.5 到 0.9 之间的可能是病毒和非病毒的混合体。很难确定区分病毒和非病毒的最佳分数,因为它取决于宿主序列和未知序列的百分比。因此,建议使用默认截止值 (0.5) 以获得最大灵敏度,然后使用 checkV 应用质量检查步骤以消除误报(预测完整性除外)。请继续下面的流程。
运行 checkV
命令解析
Score 阈值设为 0.5 时,VirSorter2 结果中可能存在一些非病毒序列或区域。因此,使用 CheckV 对 VirSorter2 的结果进行质量控制,并修剪在原病毒末端留下的潜在宿主区域。可以根据 CPU 内核的可用性调整 - t 选项。
Time: 0m28.795s
checkv end_to_end vs2-pass1/final-viral-combined.fa checkv -t 50 -d /new_data/hualin/db/checkv/checkv-db-v1.0
结果解析
- ./checkv/
- completeness.tsv
- complete_genomes.tsv
- contamination.tsv
- proviruses.fna
- quality_summary.tsv
- viruses.fna
- tmp
再次运行 VirSorter2
命令解析
- 再次利用 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
Time: 18m30.896s
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
结果解析
- ./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
命令解析
使用 DRAMv 注释鉴定的病毒序列,以用于后续人工整理。可通过 "--threads" 控制调用的 CPU 核心数。
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
结果解析
- 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 的大多数是病毒。这些经验筛选标准总结如下:
Keep1: viral_gene >0
Keep2: viral_gene =0 AND (host_gene =0 OR score >=0.95 OR hallmark >2)
Manual check: (NOT in Keep1 OR Keep2) AND viral_gene =0 AND host_gene =1 AND length >=10kb
Discard: the rest
要查看病毒基因、宿主基因、评分和序列的特征标记,您可以合并 "vs2-pass1/final-viral-score.tsv" 和 "checkv/contamination.tsv",并在电子表格中过滤。本尊为各位提供了 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
放在同一目录,终端运行如下命令。
perl get_virus.pl |
对于 Manual_check.txt
中的序列,需要用眼睛和脑袋去 dramv-annotate/annotations.tsv
中找注释信息,然后根据下面的方法判断其属于病毒还是细胞,如果是病毒,就将其所在的那行信息复制到 Virus.txt
文件的末尾,并保存。
# 依据 DRAMv 注释筛选
在病毒和宿主中都有一些共同的基因 (如脂多糖相关) 和移动元件,这些基因可能导致上述 “Keep2” 类别中的假阳性。因此,要谨慎对待带有这些基因的 contigs。专家们已经编制了一份与此相关的 “可疑” 基因列表。我们可以使用 “Keep2” 类别中的 contigs 对 DRAMv 表取子集,并在 DRAMv 子集表中筛选 “可疑” 基因 (忽略大小写,例如在 “grep” 中使用 “-i” 选项),然后将带有这些基因的 contigs 放入 “手动检查” 类别中。
# 手动检查
对于存在于“手动检查”类别中的序列,可以观察其在 "dramv-annotate/annotations.tsv" 中的注释信息。本步骤很难标准化,下面是一些经验之谈:
- 结构基因、hallmark 基因、注释缺失或假设性富集 (~10% 的基因具有 non-hypothetical 注释)。
- 缺乏 hallmarks,但 >=50% 已注释的基因为病毒,且其中至少一半以上的 viral bitcore >100,且 contig 的长度 < 50kb。
- Provirus:整合酶 / 重组酶 / 切除酶 / 阻遏子,在一侧富集了病毒基因。
- Provirus: 基因组中存在 “break”:两个基因之间的间隙对应于一个链开关,更高的编码密度,注释缺失,以及在一侧噬菌体基因的富集。
- 仅有~1-3 个基因有注释,但至少一半命中病毒,且命中基因的 bitscore 不超过病毒 bitscore 的 150% ,且 / 或病毒的 bitscore >100 。
- LPS (脂多糖)外观区域对病毒基因的命中率也非常高,bitscore>100。
- 细胞样基因是病毒基因的 3 倍,几乎所有基因都有注释,没有基因只命中病毒,也没有病毒标志基因。
- 缺乏任何病毒 hallmark genes,且长度 >50kb。
- 许多明显的细胞基因字符串,没有其他病毒标志基因。 在基准测试中遇到的例子包括 1) CRISPR Cas, 2) ABC transporters, 3) Sporulation proteins, 4) Two-component systems, 5) Secretion system。这其中一些可能是由病毒编码的,但在没有进一步证据的情况下并不表明是病毒 contig。
- 多个质粒基因或转座酶,但没有明确的只命中病毒的基因。
- 注释信息很少,仅有~1-3 个基因同时命中了病毒和细胞基因,但有 stronger bitscores 支持其为细胞基因。
- 没有强有力的命中病毒的脂多糖样区域。富含通常与脂多糖相关的基因,如外聚酶、糖 基转移酶、酰基转移酶、短 链 脱氢酶/ 还原酶、脱水酶。
- 注释为 Type IV 和 / 或 Type VI 分析系统,并被非病毒基因围绕。
- 注释信息很少,仅有~1-3 个基因全部命中细胞基因 (即使 bitscore <100) ,且没有命中病毒的基因。
最后,用户要注意,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
。
perl get_virus_seqs.pl |
# 脚本获取
关注公众号 “生信之巅”,聊天窗口回复 “1551” 获取下载链接。
敬告:使用文中脚本请引用本文网址,请尊重本人的劳动成果,谢谢!
# 参考
- VirSorter2
- Viral sequence identification SOP with VirSorter2 V.3