禁止商业或二改转载,仅供自学使用,侵权必究,如需截取部分内容请后台联系作者!
介绍
批次效应是实验中由于样本处理和测序技术变异引起的非生物学差异,可能干扰研究结果。这种效应难以完全消除,但可通过方法如PCA、PCoA等降维技术进行评估和降低其影响。在多研究或多平台数据整合时,需要特别注意消除数据差异。微生物组学领域由于数据的稀疏性,常规的校正方法如sva::ComBat和limma::removeBatchEffect不适用。
**MMUPHin(Meta-analysis via Mixed Models Utilizing Public Health Information)**是由哈佛大学Huttenhover实验室开发的微生物组数据统计分析包,特别适合处理与公共健康相关的多研究数据。它采用混合模型原理来处理批次效应,有助于减少技术变异对生物学信号的干扰。
MMUPHin处理批次效应的原理:
-
混合模型:MMUPHin使用混合模型来纳入批次效应。在这种模型中,批次效应被视为随机效应,它们与研究中的固定效应(例如治疗组与对照组)分开处理。
-
元分析方法:MMUPHin利用元分析的技术,允许来自不同研究的数据共同分析。元分析是一种统计方法,它综合并分析多个研究的结果,以获得更广泛、更全面的结论。
-
数据整合:通过混合模型和元分析方法,MMUPHin能够在考虑批次效应的同时,整合多个研究的数据,提高分析的统计能力和结论的泛化性。
-
校正批次效应:通过在模型中包括批次效应作为一个变量,MMUPHin可以校正这些非生物学差异,从而使研究结果更加可靠和准确。
加载依赖包
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
:提供可视化降维结果
-
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 获取提取码
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格式的数据对象,先对数据进行转换
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进行降维分析,并进行可视化。
# 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
结果: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分析时,研究者需要仔细考虑变量的顺序和它们之间的潜在关系。
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%(对肠道结构解释的变异)要远高于“study_condition”的R2,这表明后续基于“study_condition”的差异研究会受到“studyID”的影响,因此需要做批次校正。
批次效应校正
批次效应是一种在生物样本处理和数据分析过程中可能引入的变异,它能够对实验结果产生显著影响。尽管无法彻底消除批次效应,但可以采取措施降低其对研究结果的干扰。在处理批次效应的过程中,采用线性模型进行校正可能会对生物学意义的解释造成一定程度的影响,这可能是正面的也可能是负面的。
MMUPHin是一种先进的统计分析工具,它通过使用混合模型的方法来处理批次效应。在混合模型框架内,批次效应被归类为随机效应,这允许它们与实验设计中的固定效应(如不同治疗组别)区分开来。这种分离有助于更准确地评估固定效应对生物学问题的影响,同时控制批次效应可能带来的变异。通过这种方式,MMUPHin能够有效地调整批次效应,提高研究结果的可靠性和解释力。
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进行降维分析,并进行可视化。
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
结果:相比校正前,studyID簇在PCoA上没有明显的区分现象。
校正后的解释度
采用PERMANOVA
分析方法(MicrobiomeAnalysis::run_PERMANOVA
函数可做该分析)评估"study_condition"和"studyID"两个分组分别对整体肠道微生物的影响大小。
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”的总体变异没有太大变化。
荟萃分析
荟萃分析的主要目标是通过综合多个独立研究的结果,来得出更为全面和可靠的结论。在生物统计学领域,已有多种工具和方法被开发用于执行此类分析,其中meta
包是广泛使用的工具之一。MMUPHin软件包进一步扩展了这一功能,提供了lm_meta
函数,专门用于进行荟萃分析。该函数的分析流程通常包括以下步骤:
- 差异物种分析:首先,使用Maaslin2工具对每个独立研究中对照组(control)和结直肠癌(CRC)组之间的差异物种进行识别和分析。
- 混合模型汇总:随后,通过MMUPHin的混合模型方法,将不同研究的结果进行汇总。这种方法允许研究者纳入批次效应作为随机效应,同时对固定效应进行评估。
- 协变量校正:在分析过程中,可以引入诸如年龄、性别和体重指数(BMI)等人口统计学协变量,以校正可能的混杂因素,从而提高分析结果的准确性。
- 效应量估计:在统计模型中,
coef
参数代表效应量(Effect Size),它量化了对照组与CRC组之间在生物学指标上的差异。
计算效应值
使用MMUPHin::lm_meta
函数,计算每个species在不同研究的效应值。
if(!dir.exists("./MMUPHin_lm_meta")) {
dir.create("./MMUPHin_lm_meta", recursive = TRUE)
}
fit_lm_meta <- MMUPHin::lm_meta(
feature_abd = crc_prof,
batch = "studyID",
exposure = "study_condition",
covariates = c("gender", "age", "BMI"),
data = crc_meta,
control = list(verbose = FALSE,
output = "./MMUPHin_lm_meta"))
head(fit_lm_meta$meta_fits)
feature | exposure | coef | stderr | pval | k | |
---|---|---|---|---|---|---|
s__Faecalibacterium_prausnitzii | s__Faecalibacterium_prausnitzii | CRC | -0.028315399 | 0.015894907 | 0.07484495 | 5 |
s__Streptococcus_salivarius | s__Streptococcus_salivarius | CRC | -0.014906775 | 0.008080789 | 0.06507861 | 5 |
s__Ruminococcus_sp_5_1_39BFAA | s__Ruminococcus_sp_5_1_39BFAA | CRC | -0.031334358 | 0.013904021 | 0.02422019 | 5 |
s__Subdoligranulum_unclassified | s__Subdoligranulum_unclassified | CRC | -0.001222739 | 0.017139834 | 0.94312793 | 5 |
s__Bacteroides_stercoris | s__Bacteroides_stercoris | CRC | 0.005573002 | 0.008397279 | 0.50690311 | 5 |
s__Bifidobacterium_longum | s__Bifidobacterium_longum | CRC | -0.020595000 | 0.008893154 | 0.02056775 | 5 |
结果:5个研究的CRC荟萃分析的结果。需要关注qval.fdr和coef两个指标。
- feature:表示正在分析的物种特征。
- exposure:表示研究中考虑的暴露因素,这里是study_condition因素(control vs CRC)。
- coef:即系数,表示在控制了其他变量后,暴露因素与特征之间关系的估计效应量。正值可能表示随着暴露因素的增加,特征的水平也增加;反之亦然。
- stderr:标准误差(Standard Error),表示coef估计的精确度。较小的标准误差意味着估计更可靠。
- pval:P值,用于检验coef是否在统计上显著不同于0,即暴露因素是否真的影响了特征。
- k:表示参与荟萃分析的研究数量,这里总计5个CRC研究。
- tau2:表示随机效应方差,是不同研究效应大小不一致性的度量。
- stderr.tau2:tau2的标准误差。
- pval.tau2:tau2的P值,用于检验不同研究间效应大小的不一致性是否显著。
- I2:I平方(I-squared),表示研究间变异性的比例,用于评估荟萃分析中异质性的强度。
- H2:H平方(H-squared),是另一种衡量异质性的方法,与I平方类似,但考虑了研究的权重。
- weight:各研究的权重,通常与样本大小和效应大小的精确度有关。在荟萃分析中,权重较大的研究对总体估计的影响更大。
- weight_XXX:特定研究的权重,例如"weight_FengQ_2015"是FengQ在2015年的研究在荟萃分析中的权重。
- pval.bonf:Bonferroni校正后的P值,用于多重比较校正,减少I型错误的风险。
- qval.fdr:False Discovery Rate(FDR)校正后的q值,同样用于处理多重比较问题。
可视化结果
根据qval.fdr和coef两个指标获得差异物种。
signif_res <- fit_lm_meta$meta_fits %>%
dplyr::filter(qval.fdr < 0.05) %>%
dplyr::arrange(coef) %>%
dplyr::mutate(feature = factor(feature, levels = feature),
mycolor = ifelse(coef > 0, "CRC", "control"))
head(signif_res)
feature | exposure | coef | stderr | pval | k | |
---|---|---|---|---|---|---|
s__Faecalibacterium_prausnitzii | s__Faecalibacterium_prausnitzii | CRC | -0.028315399 | 0.015894907 | 0.07484495 | 5 |
s__Streptococcus_salivarius | s__Streptococcus_salivarius | CRC | -0.014906775 | 0.008080789 | 0.06507861 | 5 |
s__Ruminococcus_sp_5_1_39BFAA | s__Ruminococcus_sp_5_1_39BFAA | CRC | -0.031334358 | 0.013904021 | 0.02422019 | 5 |
s__Subdoligranulum_unclassified | s__Subdoligranulum_unclassified | CRC | -0.001222739 | 0.017139834 | 0.94312793 | 5 |
s__Bacteroides_stercoris | s__Bacteroides_stercoris | CRC | 0.005573002 | 0.008397279 | 0.50690311 | 5 |
s__Bifidobacterium_longum | s__Bifidobacterium_longum | CRC | -0.020595000 | 0.008893154 | 0.02056775 | 5 |
采用棒棒糖图(Lollipop Chart)可视化富集在不同组的差异物种。
pl <- ggplot(data = signif_res, aes(x = feature, y = coef, color = mycolor)) +
geom_segment(aes(x = feature, xend = feature, y = 0, yend = coef), color = "grey") +
geom_hline(yintercept = 0, linewidth = 0.1, linetype = 2) +
geom_point(size = 4) +
labs(x = "", y = "Maaslin2 (coef)") +
coord_flip() +
theme_light() +
theme(
#legend.position = "none",
panel.grid.major.x = element_blank(),
panel.border = element_blank(),
axis.ticks.x = element_blank()
)
pl
结果:5个研究荟萃分析筛选到的差异物种富集情况。Ruminococcus torques富集在CRC分组,这和Aging characteristics of colorectal cancer based on gut microbiota研究结果一致。
Eighty differential bacteria were screened out from the old group. Faecalibacterium prausnitzii, Eubacterium rectale, and Roseburia faecis were abundant in healthy individuals. Bacteroides vulgatus, Escherichia_coli, and Ruminococcus torques were enriched in CRC samples.
总结
MMUPHin是一款先进的统计分析工具,专为微生物组数据的批次效应校正而设计。它提供了一个强大的函数MMUPHin::adjust_batch
,能够对来自不同数据源的批次效应进行有效校正,确保分析结果的准确性和可靠性。
在评估批次效应的过程中,MMUPHin结合了MicrobiomeAnalysis
包中的多个函数,以实现全面和深入的分析。具体来说,MicrobiomeAnalysis::run_ord
函数用于执行基于距离的排序分析,帮助我们理解样本间的相对位置关系;MicrobiomeAnalysis::plot_ord
函数则用于绘制排序分析的结果图,直观展示样本间的分布和差异;而MicrobiomeAnalysis::run_PERMANOVA
函数则用于评估样本间差异的统计显著性,进一步验证批次效应的影响。
在完成批次效应的评估和校正后,MMUPHin进一步提供了MMUPHin::lm_meta
函数,用于进行荟萃分析。该函数能够整合多个研究的数据,通过线性模型综合评估不同研究间的效应大小和一致性,从而得到更为全面和可靠的结论。
参考
- MMUPHin