功能分析|clusterProfile包用于宏基因组学和代谢组学等多组学数据中的功能富集分析

摘要:不定时分享科研干货,设置🌟就不会错过有用的内容啦,需要的小伙伴可以及时去文末领取今日份code~

不定时分享科研干货,设置🌟就不会错过有用的内容啦,需要的小伙伴可以及时去文末领取今日份code~

科研方法分享

Sharing of scientific research methods

Y叔团队2024年7月在Nature protocols发表“使用clusterProfiler 来表征多组学数据”,该文章通过3个例子包括代谢组学和宏基因组学数据中的功能富集分析、毛竹中与耐寒性有关的转录因子分析、和单细胞转录组细胞类型注释,详细展示了怎么使用clusterProfiler包进行功能富集分析。小编对文章和代码进行学习和测试,发现这3个例子基本涵盖了功能富集分析的绝大多少内容,值得大家学习与操作,因此小编将分三期对这三个例子进行解读,期望我们的解读能让大家快速上手自己数据从而进行个性化分析。本节带给大家第一块内容,即代谢组学和宏基因组学数据中的功能富集分析。

💡输出结果展示

图片来源,Xu, S., Hu, E., Cai, Y. et al. Nat Protoc 19, 3292–3320 (2024).

比较疾病的功能谱图

在疾病研究中,整合多组学数据是深入理解生物系统复杂性和探究疾病成因的关键宏基因组学主要关注微生物的遗传组成,而代谢组学则集中于宿主-微生物相互作用的代谢产物。将这两个领域相结合提供了一种独特的视角,使我们能够全面考察疾病的复杂性及其微妙的亚型。在这一框架下,我们利用 clusterProfiler 对两种炎症性肠病(IBD)亚型的差异代谢物和差异微生物基因进行功能富集分析,并对结果进行了比较。这种整合分析可以为不同疾病条件下微生物群落的代谢潜力和功能能力提供重要见解。

ORA和GSEA算法基本步骤

主要采用 clusterProfiler ORA(在三个示例中均使用)和 GSEA 算法(在第二个示例中使用)来探索下游生物通路调控和细胞类型识别。

对于ORA分析,主要步骤包括:

(1)通过差异分析获得感兴趣的特征列表,例如差异基因或代谢物

(2)选择合适的注释数据库,例如 GO 或 KEGG 通路数据库

(3)使用 clusterProfiler 提供的富集功能(例如 enrichGO、enrichKEGG 或 enricher)使用 ORA 算法进行分析

对于GSEA,目标是获取具有生物学意义的特征的排序列表,例如按倍数变化值、T统计量或组间分析的校正P值排序的特征。然后,利用 clusterProfiler 中的 gseKEGG、gseGO 或 GSEA 等函数进行 GSEA 分析。对于涉及多个组的场景,可以使用 compareCluster 函数进行比较分析。

__

小知识点1 GSEA 或 ORA 的选择

本知识点介绍了如何在 GSEA 和 ORA 分析之间进行选择,用于指导研究者根据输入数据和分析目标进行选择。在这两种富集方法之间进行选择时,我们通常根据输入数据的特征和分析目标做出决定。对于输入数据ORA算法只需要一个感兴趣的基因列表,其中可能包括差异基因、单细胞标记基因或来自免疫沉淀相互作用的蛋白质。而GSEA需要一个完整的、排序的基因列表,例如按log2倍数变化排序的基因列表关于分析目标ORA采用 Fisher 精确检验来检查一组基因是否显著存在于一个基因集中,这是富集分析中的一种典型方法GSEA评估一个基因集是否倾向于聚集在排序列表的顶部或底部,确定该基因集在特定条件下是被激活还是被抑制。GSEA 的优势在于它能够发现微妙但一致的途径,而不仅仅依赖于差异基因。

01 示例1基本步骤详解

在本示例中,Y叔使用了预处理后的宏基因组KEGG基因丰度信息表和代谢组代谢物特征丰度表。这两种组学特征丰度表如何从原始测序数据中获得并不是本示例的重点,Y叔也提供了一个通用的处理流程,在小知识2中进行了详细阐述。在KEGG基因丰度表中,行名表示KEGG基因ID(KO),列名表示样本。在代谢物特征丰度表中,行名为KEGG化合物(COMPOUND),列名同样为样本名。此外,还需要准备一个样本分组信息表,其中行名为样本名,确保其与特征丰度表的列名一致。如果数据中的行和列反向排列,可以使用R中的t函数对特征丰度表进行转置。

在此示例的步骤3中,Y叔采用了MicrobiotaProcess包提供的mp_diff_analysis方法进行差异分析,该方法结合了Wilcoxon秩和检验和线性判别分析(LDA)来识别差异特征。也可以考虑其他差异分析方法,例如用于代谢组数据的参数检验(如t检验或方差检验),或者直接使用Wilcoxon秩和检验。对于识别宏基因组功能基因的差异特征,可以直接使用Wilcoxon秩和检验,或者采用edgeR或ALDEx2等工具。

接下来,我们使用 enrichKEGGcompareCluster函数,分别指定 organism 参数为 ko 和 cpd,对宏基因组和代谢组的差异特征进行富集分析。如果您使用的是其他常见数据库中的特征ID,例如酶委员会(EC)编号、Human Metabolome Database (HMDB)Small Molecule Pathway Database (SMPDB),可以使用自研的MicrobiomeProfiler包。该包提供了enrichKOenrichHMDBenrichSMPDB函数,作为enrichKEGG的替代方法。

____

小知识点2 代谢组学和宏基因组数据预处理

本部分主要描述了宏基因组和代谢组学数据的上游处理流程。用户可以参考以下步骤,当然也可以查阅其他已发表的相关协议以获取更多指导。

宏基因组数据处理

1.原始数据质控:使用fastp进行原始序列的初步质量控制。

2.宿主 DNA 去除:通过Bowtie2将序列比对至人类基因组,去除宿主(人类)DNA。

3.分类学注释:利用MetaPhlAn2进行分类学注释,仅保留相对丰度大于 0.1%的物种。4.功能分析:通过HUMAnN2进行功能分析,并根据KEGG KO 编号EC 编号组织基因丰度数据。

代谢组学数据处理

1. 样本分析:使用液相色谱-串联质谱技术(LC-MS/MS)分析样本代谢物,涵盖极性代谢物脂类游离脂肪酸胆汁酸

2. 高灵敏度检测:采用高灵敏度质谱仪,可检测已知和新型代谢物

3. 数据处理:包括以下步骤:噪声去除、峰检测、同位素簇检测、保留时间对齐、

通过数据库搜索进行代谢物鉴定、可使用工具 Genedata Expressionist完成上述步骤。

________

02 示例数据及R代码

⚠️!!!本次测试所用到的数据均来自Y叔的原文,对数据有什么疑问可以阅读原文。

🔑宏基因组示例数据和R代码

🌟准备数据

mg.meta.csv 是样本分组信息表,用于分组信息提取

mg.expr.csv是基因丰度信息表,用于kegg富集分析

🌟完整代码

⚠️已经安装过相关包的将包安装命令注释掉,避免耗费时间太久

⚠️代码中定义了DA函数需要注意

# 清除全局环境中的所有对象rm(list = ls)#安装和加载相关R包(需要安装的删除注释标注“#”把安装命令释放出来)pkgs "enrichplot", "ggfun", "ggplot2", "ggrepel", "ggsc", "MicrobiotaProcess", "Seurat")install.packages("BiocManager")BiocManager::install(pkgs)#加载包library(MicrobiotaProcess)library(clusterProfiler)library(ggplot2)library(enrichplot)宏基因组数据差异分析###################1、读取宏基因组数据meta_mg metagenome row.names = 1, check.name = FALSE)#读取基因丰度信息#2、定义执行差异分析的函数DA meta, # 样本元数据表 abundance = 'Abundance', # 默认指定丰度列的列名为 "Abundance" group_colname = 'Diagnosis', # 默认分组列的列名为 "Diagnosis" force = TRUE, # 是否强制执行差异分析 relative = FALSE, # 是否进行相对丰度转换subset_group, # 子集分组,指定要保留的组 diff_group, # 差异组,指定要比较的目标组 filter.p = 'pvalue', # 差异分析的 p 值筛选条件 ...) { # 允许传入其他参数 # 构造差异分析的分组列名 sign_group_colname # 将表达数据转换为 MPSE 对象 mpse # 将元数据与 MPSE 对象合并 mpse left_join(meta, by = "Sample") # 筛选样本中属于 subset_group 的组 mpse |> dplyr::filter(!!as.symbol(group_colname) %in% subset_group) |> # 使用 mp_diff_analysis 函数进行差异分析 mp_diff_analysis( .abundance = !!as.symbol(abundance), # 使用指定的丰度列 .group = !!as.symbol(group_colname), # 使用指定的分组列 force = force, # 是否强制执行分析 relative = relative, # 是否相对丰度转换 filter.p = filter.p, # 设置 p 值过滤条件 ... # 传入额外参数 ) |> # 提取差异分析得到的特征 mp_extract_feature |> # 筛选出属于 diff_group 的差异特征 dplyr::filter(!!as.symbol(sign_group_colname) == diff_group) |> # 提取 OTU 信息作为结果 dplyr::pull(OTU) |> suppressMessages # 隐藏冗余消息#鉴定差异表达基因#3、定义一个包含分组信息的列表,键为组名,值为组标记groups # 使用 lapply 遍历 groups 列表,对每个分组执行 DA 函数de_gene DA( expr = metagenome, # 传入 metagenome 数据作为表达矩阵 meta = meta_mg, # 传入元数据表 meta_mg subset_group = c(x, 'Control'), # 筛选出当前组(x)和对照组(Control) diff_group = x # 将当前组(x)指定为差异组 )})#差异基因的功能分析#4、执行基因功能富集分析,并比较不同分组的结果gene_enrich_result geneClusters = de_gene, # 输入包含差异基因的列表,按分组组织 fun = "enrichKEGG", # 使用 KEGG 数据库进行功能富集分析 organism = "ko" # 指定分析的物种为 KEGG 的 KO 编号(通用物种标识))#5、宏基因组富集结果可视化# 使用 dotplot 可视化基因富集分析结果p1 gene_enrich_result, # 富集分析结果对象 facet = "intersect", # 按交集(intersect)分面显示 showCategory = 10, # 每组显示前 10 个富集功能 split = "intersect", # 根据交集分组 label_format = 60 # 标签的最大字符宽度,超出部分换行) + ggtitle("Functional enrichment of intestinal genes") + # 设置图标题 theme( plot.title = element_text(hjust = 1) # 调整标题位置,右对齐 )p1#6、输出结果ggsave("microbiome genes.pdf",plot= p1, height = 6, width =8)write.csv(gene_enrich_result, file = "metagenome_results.csv")

🌟输出结果

microbiome genes.pdf 基因富集分析图

metagenome_results.csv 基因富集结果总计

🔑代谢组组示例数据和R代码

🌟准备数据

metabolism_meta.csv 是样本分组信息表,用于分组信息提取

metabolism_expr.csv是基因丰度信息表,用于富集分析

#安装和加载相关R包(需要安装的删除注释标注“#”把安装命令释放出来)#pkgs # "enrichplot", "ggfun", "ggplot2", "ggrepel",# "ggsc", "MicrobiotaProcess", "Seurat")#install.packages("BiocManager")#BiocManager::install(pkgs)#加载包library(MicrobiotaProcess)library(clusterProfiler)library(ggplot2)library(enrichplot)#定义执行差异分析的函数DA meta, # 样本元数据表 abundance = 'Abundance', # 默认指定丰度列的列名为 "Abundance" group_colname = 'Diagnosis', # 默认分组列的列名为 "Diagnosis" force = TRUE, # 是否强制执行差异分析 relative = FALSE, # 是否进行相对丰度转换 subset_group, # 子集分组,指定要保留的组 diff_group, # 差异组,指定要比较的目标组 filter.p = 'pvalue', # 差异分析的 p 值筛选条件 ...) { # 允许传入其他参数 # 构造差异分析的分组列名 sign_group_colname # 将表达数据转换为 MPSE 对象 mpse # 将元数据与 MPSE 对象合并 mpse left_join(meta, by = "Sample") # 筛选样本中属于 subset_group 的组 mpse |> dplyr::filter(!!as.symbol(group_colname) %in% subset_group) |> # 使用 mp_diff_analysis 函数进行差异分析 mp_diff_analysis( .abundance = !!as.symbol(abundance), # 使用指定的丰度列 .group = !!as.symbol(group_colname), # 使用指定的分组列 force = force, # 是否强制执行分析 relative = relative, # 是否相对丰度转换 filter.p = filter.p, # 设置 p 值过滤条件 ... # 传入额外参数 ) |> # 提取差异分析得到的特征 mp_extract_feature |> # 筛选出属于 diff_group 的差异特征 dplyr::filter(!!as.symbol(sign_group_colname) == diff_group) |> # 提取 OTU 信息作为结果 dplyr::pull(OTU) |> suppressMessages # 隐藏冗余消息}代谢组数据差异分析###################1、读取代谢组分组和丰度数据meta_mb metabolism row.names = 1, check.name = FALSE)#读取代谢表达丰度数据#2、鉴定差异代谢物#!!!注意只做代谢物分析的DA函数需要定义# 定义疾病分组groups # 对每个分组进行差异代谢物分析de_cpd DA( # 调用 DA 函数 expr = metabolism, # 输入代谢物表达矩阵 meta = meta_mb, # 样本元数据表 subset_group = c(x, 'Control'), # 筛选子集组,包含当前分组和对照组 diff_group = x # 当前分组为差异比较目标 )})#3、差异代谢物富集分析cpd_enrich_result geneClusters = de_cpd, # 输入每个分组的差异代谢物列表 fun = "enrichKEGG", # 调用 KEGG 富集分析函数 organism = "cpd" # 设定分析对象为 KEGG 化合物数据库)#4、代谢组结果可视化p2 cpd_enrich_result, # 富集分析结果对象 facet = 'intersect', # 根据交集分面显示每组数据的富集结果 showCategory = 10, # 显示最显著的 10 个富集通路 split = "intersect", # 按交集类型拆分显示 label_format = 60 # 设置标签长度,超过 60 个字符会换行) + ggtitle("Functional enrichment of chemical compounds") + # 设置图表标题 theme( plot.title = element_text(hjust = 1) # 图表标题右对齐 )p2#5、输出结果ggsave("metabolism.pdf",plot= p2, height = 6, width =8)write.csv(cpd_enrich_result,file="metabolism_results.csv")

🌟输出结果

metabolism.pdf 代谢物富集分析图

Metabolism_results.csv 代谢物分析结果总计

参考文献

[1] Xu, S., Hu, E., Cai, Y. et al. Using clusterProfiler to characterize multiomics data. Nat Protoc 19, 3292–3320 (2024).

往期相关内容

功能富集分析,我们之前有过几期介绍,大家可以参考下述内容:

扩增子分析|基于R语言ggpicrust2包进行PICRUSt2预测功能的分析和可视化(通路PCA,热图和差异图)

富集分析|从物种到功能:基因富集分析之GO、KEGG数据库以及过表达分析(ORA)和基因集富集分析(GSEA)

富集分析|如何解析微生物基因功能?Y叔经典R包“clusterProfiler 4.0:用于解释组学数据的通用富集工具”文章导读

📚测试代码和实例数据领取

!!!有需要的小伙伴可以去评论区自取今天的测试代码和实例数据。

📌示例代码中提供了数据和代码,已经小编测试,可直接运行。

⚠️本账号会不时更新代码,代码失效可联系小编~

以上就是clusterProfiler包用于代谢组学和宏基因组学数据中的功能富集分析的全部内容。

💬有任何问题可以后台联系,看文献,学方法,期待大家一起成长!

来源:微生物组

相关推荐