加载中...
宏病毒组分析流程1-VirSorter2
发表于:2021-12-02 | 分类: 生物信息
字数统计: 4k | 阅读时长: 19分钟 | 阅读量:

安装软件

  • 安装主程序及依赖

    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 中发现的最大病毒片段数的数字;
      • 带有标志基因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 %
    seqnamedsDNAphageNCLDVRNAssDNAlavidaviridaemax_scoremax_score_grouplengthhallmarkviralcellular
    NODE_5_length_17317_cov_8.385876||full0.9930.8470.0050.0600.4670.993dsDNAphage17315264.7005.900
    NODE_6_length_16611_cov_115.615064||full0.9200.2070.0350.0870.0530.920dsDNAphage1661003.2000.000
    NODE_8_length_14848_cov_778.417157||full1.0000.2200.1050.3800.6271.000dsDNAphage1484815100.0000.000
    NODE_16_length_12563_cov_14.331948||full0.9730.2000.1650.2730.2270.973dsDNAphage1208308.0000.000
    NODE_17_length_11885_cov_350.043956||full0.6530.5130.0500.0800.0470.653dsDNAphage1188509.1000.000
    NODE_21_length_11527_cov_216.405073||full0.6200.4070.0000.0130.0600.620dsDNAphage11526010.5005.300
    NODE_23_length_11316_cov_8.144303||full0.3670.5400.0100.0000.4000.540NCLDV11313123.1007.700

    不同病毒类群的分类器并非相互排斥,它们的目标病毒序列空间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
    seqnametrim_orf_index_starttrim_orf_index_endtrim_bp_starttrim_bp_endtrim_prtrim_pr_maxprox_orf_index_startprox_orf_index_endprox_bp_startprox_bp_endprox_prprox_pr_maxpartialfull_orf_index_startfull_orf_index_endfull_bp_startfull_bp_endpr_fullarcbaceukvirmixunalignedhallmark_cntgroupshapeseqname_new
    NODE_999_length_4026_cov_7.610929112140250.5470.54711214025nannan0112140250.5470.00.00.00.00.0100.00dsDNAphagelinearNODE_999_length_4026_cov_7.610929||full
    NODE_9999_length_1276_cov_11.59869013312740.9550.9551331274nannan013312740.9550.00.00.00.00.0100.00RNAlinearNODE_9999_length_1276_cov_11.598690||full
    NODE_99999_length_314_cov_4.0000001212870.570.57121287nannan01212870.570.00.00.00.050.050.00RNAlinearNODE_99999_length_314_cov_4.000000||full
    NODE_99992_length_314_cov_4.3899611212750.7470.747121275nannan01212750.7470.00.00.00.00.0100.00ssDNAlinearNODE_99992_length_314_cov_4.389961||full
    NODE_9997_length_1276_cov_44.11384112112400.980.981211240nannan012112400.980.00.00.00.050.050.00RNAlinearNODE_9997_length_1276_cov_44.113841||full

    原病毒provirus提取过程中,为了获得更好的敏感性,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 的结果进行质量控制,并修剪在原病毒proviruses末端留下的潜在宿主区域。可以根据 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 注释鉴定的病毒序列,以用于后续人工整理manual curation。可通过 "--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

要查看病毒基因、宿主基因、评分和序列的特征标记hallmark of sequences,您可以合并 "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 注释筛选

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

手动检查

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

  • 结构基因、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。
  • 细胞样基因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-ch)/ 还原酶(ain)脱水酶(dehydratase)
  • 注释为 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” 获取下载链接。

生信之巅微信公众号生信之巅小程序码

敬告:使用文中脚本请引用本文网址,请尊重本人的劳动成果,谢谢!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
上一篇:
利用GTDB-TK对细菌和古菌基因组进行物种分类
下一篇:
Linux中使用tar将大文件压缩为多个小的压缩包
本文目录
本文目录