病毒组学数据分析 -03 Virsorter2病毒序列识别-CSDN博客

本文链接：https://blog.csdn.net/kun_98/article/details/123056574

VirSorter2 应用多分类器、专家指导的方法来检测不同的 DNA 和 RNA 病毒基因组。它对其之前版本（virsorter）进行了重大更新：

与更多的病毒组合作，包括 dsDNA 噬菌体、ssDNA 病毒、RNA 病毒、NCLDV（Nucleocytoviricota,核细胞病毒）、lavaviridae（virophages,噬病毒体）;
应用机器学习使用基因组特征（包括结构/功能/分类注释和病毒标志基因）来估计病毒性；
使用来自宏基因组或其他来源的高质量病毒基因组进行训练。

安装

基于conda

conda create -n vs2 -c conda-forge -c bioconda virsorter=2
conda activate vs2
#安装数据库，大概需要10min
virsorter setup -d db -j 10

#若安装失败也可手动安装，https://osf.io/v46sc/download
#解压后将路径传给--db-dir
tar -xzf db.tgz
virsorter config --init-source --db-dir=./db

运行

virsorter run -w test.out -i test.fa --min-length 1500 -j 4 all

#基于不同得分（min-score）重新快速运行程序（只运行classify步骤），新文件添加后缀（label）为rerun
virsorter run -w test.out -i test.fa --include-groups "dsDNAphage,ssDNA" -j 4 --min-score 0.9 --label rerun classify

#增加hmmscan线程数提高运行速度
virsorter config --set HMMSEARCH_THREADS=4

结果解读

final-viral-combined.fa

识别的病毒序列，包括两种类型： - 确定为病毒的完整序列（用后缀||full）; - 部分序列被识别为病毒（用后缀||{i}_partial）;这里的{i}可以是从0到在该重叠群中发现的最大病毒片段数的数字; - 短（少于两个基因）序列，其标志性基因被鉴定为病毒（用后缀||lt2gene）;

final-viral-score.tsv

此表可用于进一步筛选结果。它包括以下列： - 序列名称 - 跨组的每个病毒序列的评分（多列） - 跨组的最大评分 - 最大评分组 - 重叠群长度 - 标志基因计数 - 病毒基因百分比 - 非病毒基因百分比

注意：分类器对病毒序列彼此之间不是排他性的，它们的目标病毒序列可能有重叠，这意味着不应使用或将这些信息视为可靠的分类学分类。对VirSorter2 的定位仅限于病毒识别。

final-viral-boundary.tsv

此文件中只有部分列可能有用： - seqname：原始序列名称 - trim_orf_index_start，trim_orf_index_end：已识别病毒序列的原始序列的开始和结束ORF索引 - trim_bp_start，trim_bp_end：已识别病毒序列原始序列的开始和结束位置 - trim_pr：最终修剪的病毒序列的分数 - partial：全序列为病毒或部分序列为病毒;当完整序列具有分数>分数截止值，它是完整的（0），或者在其中提取的任何病毒序列是部分（1） - pr_full：原始序列的分数 - hallmark_cnt：标志基因计数 - group：给出高分的病毒组的分类器;这不应该被用作可靠的分类

注意：VirSorter2有时倾向于在病毒前提取过程中高估病毒序列的大小，以实现更好的敏感性。建议清理这些前病毒预测，以去除预测病毒区域边缘的潜在宿主基因，例如使用CheckV工具

训练自定义分类器

见Bitbucket

建议

1. 应该如何选择分界线？

A：一般来说，得分>0.9的人是高置信度的。那些得分在 0.5 和 0.9 之间的可能是病毒和非病毒的混合物。很难找到区分病毒和非病毒的最佳分数，因为它取决于宿主序列和未知序列的百分比。因此，建议使用默认截止值 (0.5) 以获得最大灵敏度，然后使用 checkV 应用质量检查步骤来消除误报（除了预测完整性）。

2. 如何判断已鉴定的病毒序列是否为原病毒？

A：只有部分病毒序列（以_partial结尾）可以被确认为原病毒。VirSorter2 中的完全病毒序列（以 ||full 结尾）被定义为作为一个完整序列具有显着病毒信号（得分 >=0.8）的重叠群。因此有些也可能是原病毒：1）它们可能是原病毒内部的片段；2）尽管末端有一些宿主基因，但整个序列具有很强的病毒信号。

3.为什么要修剪完全病毒序列（以 ||full 结尾）？

答：在三种情况下可以修剪完整的病毒序列。

1) VirSorter2 基于被prodigal调用的基因。prodigal忽略了第一个和最后一个基因之外的一些碱基突出，默认情况下也忽略了VirSorter2；
2) 环状序列通常在一个基因的中间分裂并且有重复的片段。VirSorter2 修剪重复片段并通过将部分基因从开始移动到结束来修复分裂基因。
3）完全病毒序列仅意味着整个序列具有显着的病毒信号（默认得分> = 0.95），但VirSorter2仍然对其应用末端修剪步骤（每端10％的基因）以找到最佳病毒片段（最长默认情况下在峰值分数的 95% 以内）。同样，末端修剪步骤修剪的“完整”序列不应被解释为原病毒，因为对分数影响较小的基因，例如未知基因或宿主和病毒共有的基因，可以被修剪掉。如果希望不修剪完整序列（以 ||full 结尾）并将其留给 checkV 等专用工具，可以使用--keep-original-seq选项。

run 帮助文档

-w, --working-dir PATH        输出目录
  -d, --db-dir PATH             数据库目录
  -i, --seqfile PATH            fa或fq格式的序列文件（可以通过gzip或bz2压缩
  -l, --label TEXT              将标签作为前缀添加到输出文件中;
                                这在使用不同的筛选选项重新运行分类时很有用
  --include-groups TEXT         包括病毒组的分类器（逗号分隔，中间没有空格）;
                                可用的选项是：dsDNAphage，NCLDV，RNA，ssDNA，lavidaviridae 
                                [默认：dsDNAphage，ssDNA]
  -j, --jobs INTEGER            并行允许的最大作业数。 [默认值：160]
  --min-score FLOAT             被识别为病毒的最低评分 [默认：0.5]
  --min-length INTEGER          所需的最短序列;所有短于此数值的序列都将被删除 
                                [默认值：0]
  --keep-original-seq           保留原始序列而不是修剪;
                                默认情况下，已识别病毒序列两端的未翻译区域被修剪;
                                对环状序列进行修饰以消除两端之间的重叠，
                                并针对分裂成两端的基因进行调整 [默认值：否]
  --exclude-lt2gene             短序列（少于2个基因）不进行任何评分，
                                但那些具有标志性基因的序列默认为病毒;
                                使用此选项可排除它们 [默认：否]
  --prep-for-dramv              为DRAMv程序准备文件[默认：否]
  --high-confidence-only        仅输出高置信度病毒序列;
                                这相当于筛选具有以下标准的最终病毒 
                                score.tsv中：（max_score > = 0.9）或
                                （max_score > = 0.7并且hallmark标志基因> = 1）[默认值：否]
  --hallmark-required           所有病毒序列上都需要标志性基因[默认：否]
  --hallmark-required-on-short  需要在短病毒序列上有标记基因
                                （长度截止为"短"由template-config.yaml
                                文件中的"MIN_SIZE_ALLOWED_WO_HALLMARK_GENE"设置，
                                默认为3kbp）;这可以以合理的灵敏度成本减少误报[默认：否]
  --viral-gene-required         需要对病毒基因进行注释，去除没有注释基因的假定病毒
                                序列; 这可以以合理的敏感度成本减少误报 [默认值：否]
  --viral-gene-enrich-off       关闭调用全序列病毒的病毒基因比细胞基因更多的要求;
                                当仅使用VirSorter2产生具有从其他工具识别的病毒序列
                                的DRAMv输入时，或者当使用checkV修剪的病毒序列时，
                                这很有用[默认值：否]
  --seqname-suffix-off          关闭向序列名称添加后缀(||full, ||{i}_partial,
                                ||lt2gene)  请注意，这可能会导致来自同一重叠群的
                                部分序列具有相同的名称;当您确定每个重叠群中最大值
                                有一个部分序列时，可以使用此选项[默认值：否] 
  --provirus-off                在对完整重叠群进行分类后关闭提取原噬菌体。
                                和--max-orf-per-seq一起使用, 可以显著加快运行速度 
                                [默认值：否]
  --max-orf-per-seq INTEGER     用于计算分类特征的ORF的最大数量;
                                此选项只能在 --provirus-off 模式下使用;
                                如果序列中的ORF的超过最大限制，
                                则将其数值采样到此以减少计算 [默认值： -1]
                                可以设置ORF数量阈值
  --tmpdir TEXT                 中间文件的目录名称
  --rm-tmpdir                   删除中间文件目录（--tmpdir）; 
                                每次运行都会创建 100 多个中间文件，因此，
                                在处理 100 多个样本时，建议使用此文件，
                                以避免超出用户的文件[默认值：否]
  --verbose                     显示详细规则输出[默认值：否]
  --profile TEXT                snakemake 配置文件，例如用于集群执行。
                                如直接添加--forceall选项可用于强制重新运行
  -n, --dryrun                  检查要运行的规则和要生成的文件[默认值：否]
  --use-conda-off               停止使用此软件包附带的 conda envs （vs2.yaml），
                                并使用当前系统中安装的内容;仅当您想使用自己喜欢的
                                版本自行安装依赖项时才有用;此选项适用于开发版本 [默认：否]
  -h, --help                    显示此消息并退出
  
  all                           默认值，这意味着运行整个管道，包括 1）预处理，2）注释（特征提取）和 3）分类。
                                主要的计算瓶颈是注释步骤，占用了大约 95% 的 CPU 时间。
                                如果您只想使用不同的分数截止值 (--min-score) 重新运行，
                                则使用classify参数可以跳过注释步骤，仅重新运行分类步骤

参考文献：

GitHub - jiarong/VirSorter2: customizable pipeline to identify viral sequences from (meta)genomic data