SIAMCAT: Statistical Inference of Associations between Microbial Communities And host phenoTypes
SIAMCAT: 微生物群落和宿主表型之间关联的统计推断
Konrad Zych, Jakob Wirbel, and Georg Zeller1*
1EMBL Heidelberg
目录
SIAMCAT: Statistical Inference of Associations between Microbial Communities And host phenoTypes
1.关于本手册
本小册子旨在为SIAMCAT的主要功能提供一个简短的教程。其他工作流程的例子或更详细的教程可以在其他小册子中找到（ BioConductor page).
SIAMCAT是由Peer Bork和Georg Zeller小组在EMBL主持的计算微生物组分析工具套件的一部分。欲了解更多信息，请访问t EMBL-microbiome tools.
2.介绍
微生物组和宿主表型之间的关联最好由能够从微生物组组成预测宿主状态的定量模型来描述。 SIAMCAT可以对数百个样本的数十万个微生物分类群、基因家族或代谢途径的数据进行预测。 SIAMCAT产生图形输出，方便评估输入数据的质量和统计关联，用于模型诊断和推理，揭示最具预测性的微生物生物标志物。
3.快速开始
我们使用了SIAMCAT软件包中的一个示例数据集。作为示例数据集，我们使用了Zeller等人的出版物中的数据，该出版物证明了粪便样本中的微生物物种有可能将结直肠癌（CRC）患者与健康对照组区分开来。
library("SIAMCAT")
data("feat_crc_zeller", package="SIAMCAT")
data("meta_crc_zeller", package="SIAMCAT")首先，SIAMCAT需要一个特征矩阵（可以是矩阵、数据框架或phyloseq-otu_table），其中包含不同样品（列）的不同特征值（在行）。例如，这里包括的特征矩阵包含用mOTU分析器计算的141个样品的细菌物种的相对丰度。 mOTU profiler
feat.crc.zeller[1:3, 1:3]## CCIS27304052ST-3-0 CCIS15794887ST-4-0
## UNMAPPED 0.589839 0.7142157
## Methanoculleus marisnigri [h:1] 0.000000 0.0000000
## Methanococcoides burtonii [h:10] 0.000000 0.0000000
## CCIS74726977ST-3-0
## UNMAPPED 0.7818674
## Methanoculleus marisnigri [h:1] 0.0000000
## Methanococcoides burtonii [h:10] 0.0000000dim(feat.crc.zeller)## [1] 1754 141请注意，SIAMCAT的工作对象是相对丰度。其他类型的数据（如计数）也可以工作，但不是软件包的所有功能都能产生有意义的输出。
其次，我们在另一个数据框架中也有关于样本的元数据。
head(meta.crc.zeller)## Age BMI Gender AJCC_stage FOBT Group
## CCIS27304052ST-3-0 52 20 F -1 Negative CTR
## CCIS15794887ST-4-0 37 18 F -1 Negative CTR
## CCIS74726977ST-3-0 66 24 M -1 Negative CTR
## CCIS16561622ST-4-0 54 26 M -1 Negative CTR
## CCIS79210440ST-3-0 65 30 M -1 Positive CTR
## CCIS82507866ST-3-0 57 24 M -1 Negative CTR为了告诉SIAMCAT，哪些样本是癌症病例，哪些是健康对照，我们可以从元数据中的Group列构建一个标签对象。
label.crc.zeller <- create.label(meta=meta.crc.zeller,
label='Group', case='CRC')## Label used as case:
## CRC
## Label used as control:
## CTR## + finished create.label.from.metadata in 0.003 s现在我们有了创建SIAMCAT对象的所有材料。请看一下关于输入格式的手册 vignette about input formats ，了解更多关于支持的格式和创建SIAMCAT对象的其他方法。
sc.obj <- siamcat(feat=feat.crc.zeller,
label=label.crc.zeller,
meta=meta.crc.zeller)关于SIAMCAT对象的一些信息可以通过phyloseq的show函数访问（SIAMCAT建立在phyloseq数据结构上）。
show(sc.obj)## siamcat-class object
## label() Label object: 88 CTR and 53 CRC samples
##
## contains phyloseq-class experiment-level object @phyloseq:
## phyloseq@otu_table() OTU Table: [ 1754 taxa and 141 samples ]
## phyloseq@sam_data() Sample Data: [ 141 samples by 6 sample variables ]由于我们目前在数据集中有相当多的微生物物种，我们可以使用函数filter.features进行无监督的特征选择。
sc.obj <- filter.features(sc.obj,
filter.method = 'abundance',
cutoff = 0.001)## Features successfully filtered4.关联分析
微生物物种和标签之间的关联可以用check.association函数来测试。该函数使用非参数Wilcoxon检验和关联的不同效应大小（如AUC或折合变化）计算每个物种的显著性。
sc.obj <- check.associations(
sc.obj,
sort.by = 'fc',
alpha = 0.05,
mult.corr = "fdr",
detect.lim = 10 ^-6,
plot.type = "quantile.box",
panels = c("fc", "prevalence", "auroc"))该函数产生一个pdf文件作为输出，因为该图被优化为横向DIN-A4布局，但也可用于在一个活动的图形设备上绘图，例如在RStudio中。然后产生的绘图看起来像这样。
5.协变量检测
由于感兴趣的主要表型之外的许多生物和技术因素会影响微生物组的组成，简单的关联研究可能会受到其他变量的混杂影响，从而导致虚假的结果。check.confounders功能提供了测试相关元数据变量的潜在混杂影响的选项。在SIAMCAT对象中没有存储任何信息，但不同的分析是可视化的，并保存在一个合并的pdf文件中，用于定性解释。
sc.obj <- check.confounders(
sc.obj,
fn.plot = 'confounder_plots.pdf',
meta.in = NULL,
feature.type = 'filtered'
)条件熵检查的主要作用是在后续检查中去除无意义的变量。条件熵量化了一个变量（行）与另一个变量（列）各自包含的独特信息。分享完全相同信息的相同变量和派生变量的值将为零。在这个例子中，标签是由AJCC分期确定的组别变量派生出来的，所以两者都被排除。
条件熵图
为了更好地量化元数据变量对单个微生物特征的潜在混杂影响，check.confounder绘制了标签解释的方差与每个单个特征的元数据变量解释的方差的对比。左上角有许多特征的变量可能会混淆标签的关联。
解释的方差图
6.模型构建
SIAMCAT的一个优势是其多功能但易于使用的界面，用于在微生物物种的基础上构建机器学习模型。 SIAMCAT包含数据归一化、将数据分割成交叉验证褶皱、训练模型以及基于交叉验证实例和训练模型进行预测的功能。
6.1数据标准化
数据标准化是通过normalize.features函数进行的。在这里，我们使用log.unit方法，但也有其他一些方法和自定义选项（请查看文档）。
sc.obj <- normalize.features(
sc.obj,
norm.method = "log.unit",
norm.param = list(
log.n0 = 1e-06,
n.p = 2,
norm.margin = 1
)
)## Features normalized successfully.6.2准备交叉验证
交叉验证折叠的准备是机器学习的关键步骤。 SIAMCAT极大地简化了交叉验证方案的设置，包括对样本进行分层或根据元数据保持样本不可分割。在这个小例子中，我们选择了一个重复了两次的5倍交叉验证方案。数据拆分将被保存在SIAMCAT对象的data_split槽中。
sc.obj <- create.data.split(
sc.obj,
num.folds = 5,
num.resample = 2
)## Features splitted for cross-validation successfully.6.3训练模型
实际的模型训练是通过函数train.model进行的。同样，有多种自定义选项，从机器学习方法到模型选择的措施，或用于超参数调整的可定制参数集。
sc.obj <- train.model(
sc.obj,
method = "lasso"
)这些模型被保存在SIAMCAT对象的model_list槽中。模型的建立是通过mlr R包进行的。所有的模型都可以很容易地被访问。
# get information about the model type
model_type(sc.obj)## [1] "lasso"# access the models
models <- models(sc.obj)
models[[1]]## Model for learner.id=classif.cvglmnet; learner.class=classif.cvglmnet
## Trained on: task.id = data; obs = 112; features = 207
## Hyperparameters: nlambda=100,alpha=16.4预测
使用数据分叉和上一步训练的模型，我们可以使用函数make.predictions，以便在数据分叉的测试实例上应用模型。预测结果将被保存在SIAMCAT对象的pred_matrix槽中。
sc.obj <- make.predictions(sc.obj)
pred_matrix <- pred_matrix(sc.obj)head(pred_matrix)## CV_rep1 CV_rep2
## CCIS27304052ST-3-0 0.10011768 0.16309522
## CCIS15794887ST-4-0 0.18626022 0.27379882
## CCIS74726977ST-3-0 0.47634142 0.39204057
## CCIS16561622ST-4-0 0.31825180 0.20415432
## CCIS79210440ST-3-0 0.54648880 0.07606544
## CCIS82507866ST-3-0 0.04439485 0.092449687.模型评估及解释
在最后一部分，我们想知道模型的表现如何，以及模型中选择了哪些微生物物种。为了做到这一点，我们首先使用函数evaluate.predictions来计算预测结果与真实数据的吻合程度。这个函数计算了每个重新采样的交叉验证运行的接收操作特征（ROC）曲线下的面积（AU-ROC）和精确召回（PR）曲线。
sc.obj <- evaluate.predictions(sc.obj)## Evaluated predictions successfully.7.1图片评估
为了绘制评估结果，我们可以使用函数model.evaluation.plot，它产生一个pdf文件，显示不同重样运行的ROC和PR曲线，以及平均ROC和PR曲线。
model.evaluation.plot(sc.obj)
7.2图片解释
SIAMCAT产生的最后一幅图是模型解释图，由model.interpretation.plot函数创建。该图以柱状图的形式显示了所选择的顶级特征的
- 模型的权重（以及它们的稳健程度），以柱状图的形式显示。
- 热图，显示所选特征的Z-cores或折叠变化，以及列表显示每个模型的权重比例，这些权重是由最优先选择的特征所决定的。
此外，元数据的分布也显示在下面的热图中。
该函数再次产生一个为景观DIN-A4绘图区域优化的pdf文件。
model.interpretation.plot(
sc.obj,
fn.plot = 'interpretation.pdf',
consens.thres = 0.5,
norm.models = TRUE,
limits = c(-3, 3),
heatmap.type = 'zscore',
)结果图片如图
8.会话信息
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] SIAMCAT_1.10.0 phyloseq_1.34.0 mlr_2.18.0 ParamHelpers_1.14
## [5] BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] Biobase_2.50.0 jsonlite_1.7.1 splines_4.0.3
## [4] foreach_1.5.1 LiblineaR_2.10-8 PRROC_1.3.1
## [7] BiocManager_1.30.10 stats4_4.0.3 progress_1.2.2
## [10] yaml_2.2.1 corrplot_0.84 pillar_1.4.6
## [13] backports_1.1.10 lattice_0.20-41 glue_1.4.2
## [16] pROC_1.16.2 digest_0.6.27 RColorBrewer_1.1-2
## [19] XVector_0.30.0 checkmate_2.0.0 colorspace_1.4-1
## [22] htmltools_0.5.0 Matrix_1.2-18 plyr_1.8.6
## [25] XML_3.99-0.5 pkgconfig_2.0.3 bookdown_0.21
## [28] zlibbioc_1.36.0 purrr_0.3.4 scales_1.1.1
## [31] parallelMap_1.5.0 tibble_3.0.4 beanplot_1.2
## [34] mgcv_1.8-33 generics_0.0.2 IRanges_2.24.0
## [37] ggplot2_3.3.2 infotheo_1.2.0 ellipsis_0.3.1
## [40] BiocGenerics_0.36.0 survival_3.2-7 magrittr_1.5
## [43] crayon_1.3.4 evaluate_0.14 nlme_3.1-150
## [46] MASS_7.3-53 vegan_2.5-6 prettyunits_1.1.1
## [49] tools_4.0.3 data.table_1.13.2 hms_0.5.3
## [52] matrixStats_0.57.0 gridBase_0.4-7 lifecycle_0.2.0
## [55] BBmisc_1.11 stringr_1.4.0 Rhdf5lib_1.12.0
## [58] S4Vectors_0.28.0 munsell_0.5.0 glmnet_4.0-2
## [61] cluster_2.1.0 Biostrings_2.58.0 ade4_1.7-15
## [64] compiler_4.0.3 rlang_0.4.8 rhdf5_2.34.0
## [67] grid_4.0.3 iterators_1.0.13 rhdf5filters_1.2.0
## [70] biomformat_1.18.0 igraph_1.2.6 rmarkdown_2.5
## [73] gtable_0.3.0 codetools_0.2-16 multtest_2.46.0
## [76] reshape2_1.4.4 R6_2.4.1 gridExtra_2.3
## [79] knitr_1.30 dplyr_1.0.2 fastmatch_1.1-0
## [82] permute_0.9-5 shape_1.4.5 ape_5.4-1
## [85] stringi_1.5.3 parallel_4.0.3 Rcpp_1.0.5
## [88] vctrs_0.3.4 tidyselect_1.1.0 xfun_0.18
扫一扫
2186
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
