数据分析：宏基因组的荟萃分析之MMUPHin

2024-08-07 16:25:39 浏览数 (1)

禁止商业或二改转载，仅供自学使用，侵权必究，如需截取部分内容请后台联系作者!

数据分析：宏基因组的荟萃分析之MMUPHin​mp.weixin.qq.com/s/-XXLLkKYqlgx7j3NxAKY8w

介绍

批次效应是实验中由于样本处理和测序技术变异引起的非生物学差异，可能干扰研究结果。这种效应难以完全消除，但可通过方法如PCA、PCoA等降维技术进行评估和降低其影响。在多研究或多平台数据整合时，需要特别注意消除数据差异。微生物组学领域由于数据的稀疏性，常规的校正方法如sva::ComBat和limma::removeBatchEffect不适用。

MMUPHin（Meta-analysis via Mixed Models Utilizing Public Health Information）是由哈佛大学Huttenhover实验室开发的微生物组数据统计分析包，特别适合处理与公共健康相关的多研究数据。它采用混合模型原理来处理批次效应，有助于减少技术变异对生物学信号的干扰。

MMUPHin处理批次效应的原理：

  • 混合模型：MMUPHin使用混合模型来纳入批次效应。在这种模型中，批次效应被视为随机效应，它们与研究中的固定效应（例如治疗组与对照组）分开处理。
  • 元分析方法：MMUPHin利用元分析的技术，允许来自不同研究的数据共同分析。元分析是一种统计方法，它综合并分析多个研究的结果，以获得更广泛、更全面的结论。
  • 数据整合：通过混合模型和元分析方法，MMUPHin能够在考虑批次效应的同时，整合多个研究的数据，提高分析的统计能力和结论的泛化性。
  • 校正批次效应：通过在模型中包括批次效应作为一个变量，MMUPHin可以校正这些非生物学差异，从而使研究结果更加可靠和准确。

加载依赖包

代码语言：javascript复制
library(tidyverse)
​
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("MMUPHin")
​
library(MMUPHin)
library(phyloseq)
  • 安装MicrobiomeAnalysis包后，使用MicrobiomeAnalysis::run_ord和MicrobiomeAnalysis::plot_ord可视化数据。
  • MicrobiomeAnalysis::run_ord：提供多种降维分析函数
  • MicrobiomeAnalysis::plot_ord：提供可视化降维结果
代码语言：javascript复制
if (!requireNamespace(c("remotes", "devtools"), quietly=TRUE)) {
  install.packages(c("devtools", "remotes"))
}
remotes::install_github("HuaZou/MicrobiomeAnalysis")

导入数据

数据来源是 curatedMetagenomicData，分别来自五个研究的CRC癌症肠道微生物相对丰度数据。数据可从以下链接下载：

  • 百度云盘链接: https://pan.baidu.com/s/1ckqx-z1_WZRWJD41rCGsbQ
  • 提取码: 请关注WX公zhong号 生信学习者 后台发送 mmuphin 获取提取码
代码语言：javascript复制
crc_prof <- read.csv("CRC_profile.csv", row.names = 1)
crc_meta <- read.csv("CRC_metadata.csv", row.names = 1)
​
head(crc_meta)

subjectID

body_site

study_condition

disease

FengQ_2015.metaphlan_bugs_list.stool:SID31004

SID31004

stool

CRC

CRC;fatty_liver;hypertension

FengQ_2015.metaphlan_bugs_list.stool:SID31009

SID31009

stool

control

fatty_liver;hypertension

FengQ_2015.metaphlan_bugs_list.stool:SID31021

SID31021

stool

control

NA

FengQ_2015.metaphlan_bugs_list.stool:SID31071

SID31071

stool

control

fatty_liver

FengQ_2015.metaphlan_bugs_list.stool:SID31112

SID31112

stool

control

fatty_liver

FengQ_2015.metaphlan_bugs_list.stool:SID31129

SID31129

stool

control

fatty_liver;T2D;hypertension

数据预处理

run_ord函数的输入数据是phyloseq格式的数据对象，先对数据进行转换

代码语言：javascript复制
crc_meta$studyID <- gsub(".metaphlan_bugs_list.stool", "", crc_meta$studyID)
crc_meta$study_condition <- factor(crc_meta$study_condition, 
                                   levels = c("control", "CRC"))
​
colnames(crc_prof) <- gsub("ZellerG_2014.metaphlan_bugs_list.stool.|YuJ_2015.metaphlan_bugs_list.stool.|VogtmannE_2016.metaphlan_bugs_list.stool.|HanniganGD_2017.metaphlan_bugs_list.stool.|FengQ_2015.metaphlan_bugs_list.stool.", "", colnames(crc_prof))
rownames(crc_meta) <- gsub("ZellerG_2014.metaphlan_bugs_list.stool.|YuJ_2015.metaphlan_bugs_list.stool.|VogtmannE_2016.metaphlan_bugs_list.stool.|HanniganGD_2017.metaphlan_bugs_list.stool.|FengQ_2015.metaphlan_bugs_list.stool.", "", rownames(crc_meta))
​
​
rownames(crc_meta) <- make.names(rownames(crc_meta))
colnames(crc_prof) <- make.names(colnames(crc_prof))
rownames(crc_prof) <- make.names(rownames(crc_prof))
tax_tab <- data.frame(Species = rownames(crc_prof))
​
rownames(tax_tab) <- tax_tab$Species
​
crc_phy <- phyloseq::phyloseq(
  sample_data(crc_meta),
  otu_table(crc_prof, taxa_are_rows = TRUE),
  tax_table(as.matrix(tax_tab)))
​
crc_phy
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 484 taxa and 551 samples ]
sample_data() Sample Data:       [ 551 samples by 28 sample variables ]
tax_table()   Taxonomy Table:    [ 484 taxa by 1 taxonomic ranks ]

降维分析

使用PCoA进行降维分析，并进行可视化。

代码语言：javascript复制
# PCoA降维分析
ord_result <- MicrobiomeAnalysis::run_ord(
  object = crc_phy,
  variable = "study_condition",
  method = "PCoA")
​
# 可视化降维结果
PCoA_pl_raw <- MicrobiomeAnalysis::plot_ord(
  reslist = ord_result,
  variable = "studyID",
  #variable_color = c("red", "blue"),
  var_shape = "study_condition",
  ellipse_type = "ellipse_groups")
​
PCoA_pl_raw

添加图片注释，不超过 140 字（可选）

结果：5个CRC研究的健康和CRC分布情况。图中五个研究的样本存在不同的簇，这可能导致不同研究的批次效应影响到生物学效应（Control和CRC之间的差异物种识别）。

评估批次效应

采用PERMANOVA分析方法（MicrobiomeAnalysis::run_PERMANOVA函数可做该分析）评估"study_condition"和"studyID"两个分组分别对整体肠道微生物的影响大小。

PERMANOVA（Permutational Multivariate Analysis of Variance）是一种用于分析和比较多元数据集中不同群组间差异的统计方法，特别适用于生态和环境科学领域。它基于排列方法来检验数据中组间和组内的方差分配是否有显著差异，不要求数据的正态分布或方差齐性。 PERMANOVA的核心原理是利用任何形式的距离度量来比较不同群组间的平均距离。它使用排列测试来评估群组间的差异，通过计算观测值之间的距离，然后进行一定次数的随机置换，来得到p值，以判断不同组之间是否显著不同。 结果解释方面，PERMANOVA的输出主要包括F值、R²值和p值。F值反映了各组之间差异的大小，F值越大则差异越明显；R²值表示差异的解释程度，R²值越高则解释程度越高；P值为显著性检验结果，通常设定P值小于0.05，则认为差异显著。 需要注意的是，虽然PERMANOVA对数据分布没有严格要求，但它会对于群组间离散度（variability）的差异敏感。如果检测出显著性差异，这可能是由于群组间位置（location）、离散度或两者的组合所导致的。因此，在解释PERMANOVA结果时，建议同时评估数据离散度的影响，例如使用PERMDISP进行额外的检验。 此外，PERMANOVA分析时，变量的顺序可能会影响结果，因为它采用的是序贯检验方式，即先评估第一个变量解释的差异比例，再评估后续变量解释的剩余总体差异的比例。在多因素分析中，这种顺序效应可能会导致对不同因素重要性的解释产生偏差。因此，在进行PERMANOVA分析时，研究者需要仔细考虑变量的顺序和它们之间的潜在关系。

代码语言：javascript复制
MicrobiomeAnalysis::run_PERMANOVA(
  crc_phy,
  variables = c("study_condition", "studyID"),
  mode = "one",
  method = "bray")

SumsOfSample

Df

SumsOfSqs

MeanSqs

F.Model

R2

Pr(>F)

AdjustedPvalue

study_condition

551

1

1.542097

1.542097

5.111122

0.009224003

0.001

0.001

studyID

551

4

13.359955

3.339989

11.855396

0.079912133

0.001

0.001

结果：PERMANOVA的Pr(>F) < 0.05表明“study_condition”和“studyID”均与整体肠道微生物结构有显著差异。但与此同时，“studyID”的R2是8%(对肠道结构解释9.1%的变异)要远高于“study_condition”的R2，这表明后续基于“study_condition”的差异研究会受到“studyID”的影响，因此需要做批次校正。

批次效应校正

批次效应是一种在生物样本处理和数据分析过程中可能引入的变异，它能够对实验结果产生显著影响。尽管无法彻底消除批次效应，但可以采取措施降低其对研究结果的干扰。在处理批次效应的过程中，采用线性模型进行校正可能会对生物学意义的解释造成一定程度的影响，这可能是正面的也可能是负面的。

MMUPHin是一种先进的统计分析工具，它通过使用混合模型的方法来处理批次效应。在混合模型框架内，批次效应被归类为随机效应，这允许它们与实验设计中的固定效应（如不同治疗组别）区分开来。这种分离有助于更准确地评估固定效应对生物学问题的影响，同时控制批次效应可能带来的变异。通过这种方式，MMUPHin能够有效地调整批次效应，提高研究结果的可靠性和解释力。

代码语言：javascript复制
fit_adjust_batch <- MMUPHin::adjust_batch(
  feature_abd = crc_prof,
  batch = "studyID",
  covariates = "study_condition",
  data = crc_meta,
  control = list(verbose = FALSE))
​
crc_prof_adj <- fit_adjust_batch$feature_abd_adj
​
head(crc_prof_adj[, 1:6])
                                 SID31004    SID31009     SID31021    SID31071     SID31112     SID31129
s__Faecalibacterium_prausnitzii 0.10120482 0.050692229 0.0870712280 0.306298806 9.138322e-03 0.0412915623
s__Streptococcus_salivarius     0.06044265 0.001458136 0.0018225717 0.003949364 4.243000e-03 0.0006752066
s__Ruminococcus_sp_5_1_39BFAA   0.02374596 0.016513626 0.0169870545 0.037752728 0.000000e 00 0.0826998372
s__Subdoligranulum_unclassified 0.03265566 0.039447880 0.0753412636 0.050567518 3.055722e-02 0.0071549150
s__Bacteroides_stercoris        0.31510103 0.000000000 0.0004372288 0.003220664 5.676591e-06 0.0000000000
s__Bifidobacterium_longum       0.01595582 0.016304677 0.0260266972 0.044504423 1.641349e-05 0.0036925262

校正后批次效应

使用PCoA进行降维分析，并进行可视化。

代码语言：javascript复制
rownames(crc_meta) <- make.names(rownames(crc_meta))
colnames(crc_prof_adj) <- make.names(colnames(crc_prof_adj))
rownames(crc_prof_adj) <- make.names(rownames(crc_prof_adj))
tax_tab_adj <- data.frame(Species = rownames(crc_prof_adj))
rownames(tax_tab_adj) <- tax_tab_adj$Species
​
crc_phy_adj <- phyloseq::phyloseq(
  sample_data(crc_meta),
  otu_table(crc_prof_adj, taxa_are_rows = TRUE),
  tax_table(as.matrix(tax_tab_adj)))
​
# 运行
ord_result_adj <- MicrobiomeAnalysis::run_ord(
  object = crc_phy_adj,
  variable = "study_condition",
  method = "PCoA")
​
PCoA_pl_adj <- MicrobiomeAnalysis::plot_ord(
  reslist = ord_result_adj,
  variable = "studyID",
  # variable_color = c("red", "blue"),
  var_shape = "study_condition",
  ellipse_type = "ellipse_groups")
​
PCoA_pl_adj

添加图片注释，不超过 140 字（可选）

结果：相比校正前，studyID簇在PCoA上没有明显的区分现象。

校正后的解释度

采用PERMANOVA分析方法（MicrobiomeAnalysis::run_PERMANOVA函数可做该分析）评估"study_condition"和"studyID"两个分组分别对整体肠道微生物的影响大小。

代码语言：javascript复制
MicrobiomeAnalysis::run_PERMANOVA(
  crc_phy_adj,
  variables = c("study_condition", "studyID"),
  mode = "one",
  method = "bray")

SumsOfSample

Df

SumsOfSqs

MeanSqs

F.Model

R2

Pr(>F)

AdjustedPvalue

study_condition

551

1

1.498815

1.498815

5.118839

0.00923780

0.001

0.001

studyID

551

4

4.937974

1.234493

4.284743

0.03043471

0.001

0.001

结果：相比校正前，“studyID”的解释肠道结构总体变异度从8%降低到了3%。与此同时，“study_condition”的总体变异没有太大变化。

更多内容见：

数据分析：宏基因组的荟萃分析之MMUPHin

总结

MMUPHin是一款先进的统计分析工具，专为微生物组数据的批次效应校正而设计。它提供了一个强大的函数MMUPHin::adjust_batch，能够对来自不同数据源的批次效应进行有效校正，确保分析结果的准确性和可靠性。

在评估批次效应的过程中，MMUPHin结合了MicrobiomeAnalysis包中的多个函数，以实现全面和深入的分析。具体来说，MicrobiomeAnalysis::run_ord函数用于执行基于距离的排序分析，帮助我们理解样本间的相对位置关系；MicrobiomeAnalysis::plot_ord函数则用于绘制排序分析的结果图，直观展示样本间的分布和差异；而MicrobiomeAnalysis::run_PERMANOVA函数则用于评估样本间差异的统计显著性，进一步验证批次效应的影响。

在完成批次效应的评估和校正后，MMUPHin进一步提供了MMUPHin::lm_meta函数，用于进行荟萃分析。该函数能够整合多个研究的数据，通过线性模型综合评估不同研究间的效应大小和一致性，从而得到更为全面和可靠的结论。

参考

  • MMUPHin
数据分析

0 人点赞