宏病毒组量化指南：从丰度计算到差异挖掘

天意生信云

发布于 2025-12-24 15:34:42

1310

文章被收录于专栏：天意生信俱乐部天意生信俱乐部

在上一篇教程中，我们利用 iPHoP 成功为病毒找到了宿主。现在，我们手中已经掌握了病毒的身份（Taxonomy）和关系（Host-Virus link）。

接下来的核心问题是：“这些病毒在不同样本中发生了什么变化？”

病人肠道里的某种噬菌体是否比健康人更多？
海洋深层的病毒多样性是否低于表层？
哪些特定的“病毒-宿主”互作模式在处理组中显著增强？

为了回答这些问题，我们需要进行定量与差异分析。本篇将带你从原始测序数据出发，计算病毒丰度，并利用 R 语言进行生态统计。

第一部分：给病毒“数人头” —— 丰度计算

要比较差异，首先得知道“有多少”。在宏病毒组中，我们通常将测序得到的 Clean Reads 比对回组装好的病毒序列（vOTUs），根据比对上的 Reads 数量来计算丰度。

核心工具与原理

Bowtie2: 目前最常用的短序列比对工具，速度快，内存占用低。用于将 Reads 比对到 vOTUs 上。
SAMtools: 处理比对结果（SAM/BAM格式）的标准工具，用于排序和索引。
CoverM: （强烈推荐）这是一个专门为宏基因组设计的封装工具，能自动调用 Bowtie2 并直接输出标准化的丰度表（TPM, RPKM, Count），省去了很多手写脚本的麻烦。

环境部署

# 推荐创建一个新的环境用于比对
mamba create -n mapping -c bioconda -c conda-forge coverm bowtie2 samtools
conda activate mapping

实战操作：一步生成丰度表

前提条件：

vOTUs.fasta: 你的病毒序列集合（去冗余后的）。
SampleA_1.fq.gz, SampleA_2.fq.gz: 样本 A 的双端测序数据。
SampleB_1.fq.gz, SampleB_2.fq.gz: 样本 B 的双端测序数据。

传统方法（Bowtie2 + Samtools）：

# 1. 建库
bowtie2-build vOTUs.fasta vOTUs_index

# 2. 比对 (以 SampleA 为例)
bowtie2 -x vOTUs_index -1 SampleA_1.fq.gz -2 SampleA_2.fq.gz -S SampleA.sam -p 16

# 3. 格式转换与排序
samtools view -bS SampleA.sam | samtools sort -o SampleA.sorted.bam

# 4. 统计 Reads 数 (idxstats)
samtools index SampleA.sorted.bam
samtools idxstats SampleA.sorted.bam > SampleA.counts.txt

进阶方法（CoverM，推荐）：

CoverM 能够自动完成建库、比对、过滤和计算，一步到位。

# 计算 TPM (Transcripts Per Million) 和 Raw Counts
# --min-read-percent-identity 95: 要求 Read 必须 95% 相似才算比对上（严格质控）
# --min-covered-fraction 0.75: 要求病毒基因组至少被覆盖 75% 才算存在（防假阳性关键参数！）

coverm contig \
    --reference vOTUs.fasta \
    --coupled SampleA_1.fq.gz SampleA_2.fq.gz SampleB_1.fq.gz SampleB_2.fq.gz \
    --min-read-percent-identity 95 \
    --min-covered-fraction 0.75 \
    --methods tpm count \
    --threads 16 \
    --output-file vOTUs_abundance_table.tsv

结果解析与质控

运行后，你会得到 vOTUs_abundance_table.tsv。

数据示例：

如何解读与“排雷”：

TPM vs Count:
- TPM (Transcripts Per Million): 归一化后的相对丰度，消除了测序深度和基因长度的影响。适合做样本间比较（如 PCoA、多样性）。
- Count (Raw Reads): 原始比对数。适合做差异分析（DESeq2 要求输入整数）。
结果是否正常？
- 0值过多（Zero-inflation）：宏病毒组数据非常稀疏，如果表格里有 70%-90% 的 0 是正常的。
- 比对率过低：查看 Bowtie2 的日志。如果 Overall alignment rate < 1%，说明测序数据里几乎没有你组装出来的病毒，或者宿主污染极高。正常环境样本通常在 5% - 50% 之间（取决于是否富集过病毒颗粒）。

第二部分：组间差异分析 —— 统计与可视化

拿到了丰度表，我们转战 R 语言。这一步的目标是回答：“组间群落结构是否有显著差异？” 以及 “哪个病毒导致了这种差异？”

软件工具

R 语言: 数据分析平台。
vegan: 生态学分析的神器（多样性计算）。
DESeq2: 基于负二项分布的差异分析工具（原用于转录组，现广泛用于宏基因组/宏病毒组）。

准备数据

需要两个文件：

counts.txt: 行为 vOTU，列为样本的原始 Count 表。
group.txt: 样本分组信息（如 Control vs Treatment）。

# R 代码开始
library(vegan)
library(DESeq2)
library(ggplot2)
library(pheatmap)

# 1. 读取数据
counts <- read.table("vOTUs_abundance_table_counts_only.tsv", header=T, row.names=1, sep="\t")
metadata <- read.table("group.txt", header=T, row.names=1, sep="\t")

# 确保样本名称一致
counts <- counts[, rownames(metadata)]

多样性分析 (Diversity)

α-多样性 (Alpha Diversity)

反映样本内部的病毒丰富度。

# 计算 Shannon 指数
alpha_div <- diversity(t(counts), index = "shannon")

# 结合分组画箱线图
plot_data <- data.frame(Sample = names(alpha_div), Shannon = alpha_div, Group = metadata$Group)
ggplot(plot_data, aes(x=Group, y=Shannon, fill=Group)) + 
  geom_boxplot() + 
  geom_jitter(width=0.2) +
  theme_bw() +
  ggtitle("Viral Alpha Diversity")

结果怎么看？如果 Treatment 组的 Shannon 指数显著低于 Control 组，说明处理导致了病毒群落变得单一（可能是某种优势病毒爆发）。

β-多样性 (Beta Diversity)

反映样本间的群落结构差异。

# 计算 Bray-Curtis 距离矩阵 (通常建议先对 count 做标准化，这里简化演示)
dist_mat <- vegdist(t(counts), method = "bray")

# PCoA 分析
pcoa <- cmdscale(dist_mat, k=2, eig=TRUE)
pcoa_points <- data.frame(pcoa$points)
colnames(pcoa_points) <- c("PCoA1", "PCoA2")
pcoa_points$Group <- metadata$Group

# 绘图
ggplot(pcoa_points, aes(x=PCoA1, y=PCoA2, color=Group)) +
  geom_point(size=3) +
  stat_ellipse() +  # 加置信椭圆
  theme_bw() +
  ggtitle("PCoA Analysis of Viral Communities")

结果是否正常？

正常：同组样本（颜色相同）聚在一起，不同组分开。
异常：所有点混在一起，或者按照“提取批次”聚类（说明有批次效应，需去除）。

寻找差异病毒 (Differential Analysis with DESeq2)

这是文章中最“值钱”的结果——找到标志性病毒（Biomarker）。

# 1. 构建 DESeq 对象
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~ Group)

# 2. 运行差异分析
dds <- DESeq(dds)

# 3. 提取结果 (假设比较 Treatment vs Control)
res <- results(dds, contrast=c("Group", "Treatment", "Control"))

# 4. 筛选显著差异的 vOTUs (阈值：Padj < 0.05 且 差异倍数 > 2)
diff_vOTUs <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

# 5. 保存结果
write.csv(as.data.frame(diff_vOTUs), "Differential_vOTUs.csv")

数据结果解析：生成的 CSV 文件中重点看以下两列：