生信分析实战：用MetaPhlAn4处理完测序数据后，这些结果文件怎么用？（附常用脚本）-编程阁

MetaPhlAn4结果文件深度解析：从数据提取到高级可视化的完整指南

当你第一次拿到MetaPhlAn4生成的.txt结果文件时，可能会被那些看似晦涩的clade_name和relative_abundance搞得一头雾水。别担心，这篇文章将带你从零开始理解这些数据，并教你如何将它们转化为直观的图表和有价值的科学发现。

1. 理解MetaPhlAn4的输出结构

MetaPhlAn4生成的.txt文件实际上是一个结构化的微生物群落组成报告。让我们先来看一个典型的输出示例：

#SampleID Metaphlan_Analysis clade_name NCBI_tax_id relative_abundance additional_species k__Bacteria 2 99.92666 k__Archaea 2157 0.07334 k__Bacteria|p__Proteobacteria 2|1224 89.13284 k__Bacteria|p__Actinobacteria 2|201174 8.93442

每一列的含义如下：

clade_name：微生物的分类信息，采用层级标记法（k__表示界，p__表示门，c__表示纲，o__表示目，f__表示科，g__表示属，s__表示种）
NCBI_tax_id：对应分类单元的NCBI分类ID
relative_abundance：该分类单元在样本中的相对丰度（百分比）
additional_species：当MetaPhlAn无法精确到种水平时，可能包含的其他物种信息

关键点：MetaPhlAn4使用"|"符号连接不同分类层级，这种结构化的命名方式为后续的数据提取和处理提供了便利。

2. 结果文件的预处理技巧

在进行分析前，我们通常需要将多个样本的结果合并并进行预处理。MetaPhlAn4自带的merge_metaphlan_tables.py脚本可以轻松完成这项工作：

# 激活conda环境（假设环境名为metaphlan4） conda activate metaphlan4 # 合并所有样本结果 merge_metaphlan_tables.py sample1.txt sample2.txt sample3.txt > merged_abundance.txt

合并后的文件是一个矩阵，行是分类单元，列是不同样本，单元格内是对应的相对丰度值。这种格式非常适合后续的统计分析和可视化。

2.1 按分类层级提取数据

在实际分析中，我们经常需要关注特定分类层级（如门、属或种）的微生物组成。以下是一组实用的命令行工具组合，可以快速提取不同层级的丰度表：

# 提取种水平数据 grep -E '(s__)|(clade_name)' merged_abundance.txt | grep -v 't__' | sed 's/^.*s__//g' | awk '{$2=null;print}' | sed 's/\ \ /\ /g' | sed 's/\ /\t/g' > species_abundance.txt # 提取属水平数据 grep -E '(g__)|(clade_name)' merged_abundance.txt | grep -v 's__' | sed 's/^.*g__//g' | awk '{$2=null;print}' | sed 's/\ \ /\ /g' | sed 's/\ /\t/g' > genus_abundance.txt

提示：这些命令利用了grep、sed和awk的组合，通过模式匹配和文本处理提取特定层级的分类信息。可以根据需要修改命令中的层级标记（如将g__改为p__来提取门水平数据）。

3. 在R中进行统计分析与可视化

将处理好的数据导入R后，我们可以进行更深入的分析和可视化。以下是使用R语言进行微生物群落分析的完整流程。

3.1 数据导入与预处理

# 加载必要的包 library(tidyverse) library(phyloseq) library(vegan) library(ggplot2) # 读取数据 abundance_data <- read.table("species_abundance.txt", header=TRUE, sep="\t", row.names=1) # 转换为phyloseq对象（假设有元数据文件metadata.csv） metadata <- read.csv("metadata.csv", row.names=1) OTU <- otu_table(abundance_data, taxa_are_rows=TRUE) physeq <- phyloseq(OTU, sample_data(metadata))

3.2 α多样性分析

α多样性反映单个样本内的微生物多样性，常用的指标包括Shannon指数和Simpson指数：

# 计算α多样性 alpha_div <- estimate_richness(physeq, measures=c("Shannon", "Simpson")) # 将结果与元数据合并 alpha_div <- cbind(alpha_div, sample_data(physeq)) # 绘制箱线图 ggplot(alpha_div, aes(x=Group, y=Shannon, fill=Group)) + geom_boxplot() + theme_minimal() + labs(title="Shannon Diversity Index by Group")

3.3 β多样性分析

β多样性反映样本间的微生物组成差异，常用主坐标分析(PCoA)展示：

# 计算Bray-Curtis距离 dist_matrix <- phyloseq::distance(physeq, method="bray") # PCoA分析 pcoa <- ordinate(physeq, method="PCoA", distance=dist_matrix) # 可视化 plot_ordination(physeq, pcoa, color="Group") + geom_point(size=3) + stat_ellipse() + theme_minimal()

3.4 堆叠柱状图绘制

堆叠柱状图是展示微生物群落组成的直观方式：

# 按门水平聚合数据 physeq_phylum <- tax_glom(physeq, taxrank="Phylum") # 转换为相对丰度 physeq_phylum_rel <- transform_sample_counts(physeq_phylum, function(x) x/sum(x)) # 提取前10个最丰富的门 top_phyla <- names(sort(taxa_sums(physeq_phylum_rel), decreasing=TRUE)[1:10]) physeq_top <- prune_taxa(top_phyla, physeq_phylum_rel) # 绘制堆叠柱状图 plot_bar(physeq_top, fill="Phylum") + facet_wrap(~Group, scales="free_x") + theme(axis.text.x=element_blank()) + labs(y="Relative Abundance", title="Top 10 Phyla Composition")

4. 在Python中的高级分析

对于习惯使用Python的研究者，以下是如何利用Python进行MetaPhlAn4结果分析的示例。

4.1 数据导入与基本分析

import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt # 读取数据 abundance = pd.read_csv("merged_abundance.txt", sep="\t", index_col=0) # 数据预处理：过滤低丰度分类单元 abundance = abundance[abundance.mean(axis=1) > 0.01] # 保留平均丰度>1%的分类单元 # 计算样本间的相关性 corr_matrix = abundance.corr() sns.clustermap(corr_matrix, cmap="vlag", center=0) plt.title("Sample Correlation Heatmap") plt.show()

4.2 热图可视化

热图可以直观展示不同样本中微生物的分布模式：

# 标准化数据（Z-score） abundance_zscore = abundance.apply(lambda x: (x-x.mean())/x.std(), axis=1) # 绘制热图 plt.figure(figsize=(10, 8)) sns.heatmap(abundance_zscore, cmap="RdBu_r", center=0) plt.title("Microbial Abundance Z-scores") plt.show()

4.3 机器学习应用

我们可以使用机器学习方法识别不同组别间的特征微生物：

from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report # 假设metadata包含分组信息 X = abundance.T y = metadata['Group'] # 划分训练测试集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3) # 训练随机森林模型 clf = RandomForestClassifier(n_estimators=100) clf.fit(X_train, y_train) # 评估模型 y_pred = clf.predict(X_test) print(classification_report(y_test, y_pred)) # 查看特征重要性 feature_importance = pd.Series(clf.feature_importances_, index=abundance.index) top_features = feature_importance.sort_values(ascending=False)[:10] top_features.plot(kind='barh') plt.title("Top 10 Important Microbial Features") plt.show()

5. 与功能预测工具的联用分析

MetaPhlAn4的物种组成分析结果可以与其他工具联用，预测微生物群落的功能潜力。PICRUSt2是最常用的功能预测工具之一。

5.1 PICRUSt2分析流程

# 将MetaPhlAn4结果转换为PICRUSt2输入格式 metaphlan_to_stamp.pl -f merged_abundance.txt > metaphlan4.spf # 运行PICRUSt2 picrust2_pipeline.py -i metaphlan4.spf -o picrust2_out \ --processes 4 --stratified

5.2 功能分析结果可视化

在R中可视化PICRUSt2预测的通路差异：

# 读取PICRUSt2结果 pathway <- read.table("picrust2_out/pathways_out/path_abun_unstrat.tsv", header=TRUE, row.names=1, sep="\t") # 差异分析（假设有Control和Case两组） pathway_filtered <- pathway[rowMeans(pathway) > 0.001, ] pvals <- apply(pathway_filtered, 1, function(x) wilcox.test(x ~ metadata$Group)$p.value) sig_pathways <- pathway_filtered[p.adjust(pvals, method="fdr") < 0.05, ] # 绘制热图 pheatmap::pheatmap(log10(sig_pathways+1e-5), annotation_col=data.frame(Group=metadata$Group), show_rownames=FALSE, main="Differentially Abundant Pathways")

6. 高级可视化技巧

除了基本的堆叠柱状图和热图外，还有一些高级可视化方法可以更生动地展示微生物群落数据。

6.1 网络分析

微生物共现网络可以揭示物种间的相互作用关系：

library(igraph) library(psych) # 计算物种间的Spearman相关性 cor_matrix <- corr.test(t(abundance_data), method="spearman")$r # 只保留显著相关（p<0.01） cor_matrix[cor_matrix < 0.7] <- 0 # 构建网络图 net <- graph_from_adjacency_matrix(cor_matrix, weighted=TRUE, mode="undirected") # 简化网络 net <- simplify(net, remove.multiple=TRUE, remove.loops=TRUE) # 设置可视化参数 V(net)$size <- log(degree(net)+1)*2 E(net)$width <- abs(E(net)$weight)*3 # 绘制网络图 plot(net, layout=layout_with_fr, vertex.label.cex=0.6, vertex.color="lightblue", main="Microbial Co-occurrence Network")

6.2 三元相图

三元相图可以展示三个主要分类单元在不同样本中的相对比例：

library(ggtern) # 选择三个最丰富的门 top3_phyla <- names(sort(rowMeans(abundance_data), decreasing=TRUE)[1:3]) tern_data <- as.data.frame(t(abundance_data[top3_phyla, ])) # 绘制三元图 ggtern(tern_data, aes(x=get(top3_phyla[1]), y=get(top3_phyla[2]), z=get(top3_phyla[3]), color=metadata$Group)) + geom_point(size=3) + theme_rgbw() + labs(title="Ternary Plot of Top 3 Phyla")

7. 实战案例分析

让我们通过一个实际案例来整合前面学到的技术。假设我们有一组肠道微生物组数据，包含健康对照组和疾病组各15个样本。

7.1 差异物种分析

首先识别两组间显著差异的微生物：

library(DESeq2) # 创建DESeq2对象 dds <- DESeqDataSetFromMatrix(countData=round(abundance_data*1e6), colData=metadata, design=~Group) # 运行差异分析 dds <- DESeq(dds) res <- results(dds, contrast=c("Group", "Disease", "Control")) # 提取显著结果 sig_species <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1) write.csv(as.data.frame(sig_species), "differential_species.csv")

7.2 构建预测模型

利用差异物种构建疾病预测模型：

from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score, roc_curve # 准备数据 X = abundance.loc[sig_species.index].T y = [1 if g == "Disease" else 0 for g in metadata['Group']] # 训练逻辑回归模型 model = LogisticRegression(penalty='l1', solver='liblinear') model.fit(X, y) # 评估模型 prob = model.predict_proba(X)[:,1] print("AUC:", roc_auc_score(y, prob)) # 绘制ROC曲线 fpr, tpr, _ = roc_curve(y, prob) plt.plot(fpr, tpr) plt.plot([0,1], [0,1], 'k--') plt.xlabel("False Positive Rate") plt.ylabel("True Positive Rate") plt.title("ROC Curve for Disease Prediction") plt.show()