MIMOSA2：基于微生物组和代谢组数据的整合分析

MIMOSA2 升级自MIMOSA1。是 Borenstein 实验室(http://borensteinlab.com/ , 专注宏基因组系统 生物学)最新开发的工具。用于微生物群落和代谢组的整合分析，寻找微生物和代谢产物之间的关系。

先前Borenstein bab开发过Fishtaco工具(http://borensteinlab.com/software_fishtaco.html)，用于微生物群落和功能联合分析，[Fishtaco工具无论是扩增子测序数据还是宏基因组数据都可以使用，只是功能和微生物群落数据的对应关系不同，详细参见我自己写的Fishtaco教程，这也是目前唯一一份中文教程:参见Fishtaco](https://mp.weixin.qq.com/s/DBSFnTYGofRgEgGn7QGWjQ)。本次MIMOSA2用于整合微生物群落并在功能层面上耦合代谢数据，希望可以找到关键微生物和特定代谢产物的之间的关系。作者为此还开发了网页版本的工具：https://elbo-spice.gs.washington.edu/shiny/MIMOSA2shiny/ ，可以先尝鲜。

MIMOSA2有一个数据库，包含物种及其对应的代谢能力信息，我们的扩增子测序数据可以通过不同的方式映射至数据库。让我欣喜的是无论是基于去噪方式的ASV还是传统聚类方式的OTU，都可以很好地进行分析。这也是随着新的聚类方式发展的新式工具之一。

MIMOSA2通过微生物群落特征和代谢组特征关联，回答以下几个问题：

1. 代谢组数据是否和微生物组数据紧密相关；

2. 微生物群落差异是否可以解释代谢差异；

3. 何种微生物或者基因造成了代谢差异，或者贡献了代谢差异。

MIMOSA2的工作原理

工作原理：首先得到微生物组数据，使用参考数据库构建代谢模型；计算不同微生物对不同代谢物的潜在代谢潜能指标(CMP)；对每种代谢产物全部潜在贡献微生物CMP值进行汇总并做回归分析，评估微生物群落对指定代谢物的预测能力；最后筛选对特定代谢物影响加大的微生物作为关键微生物。

MIMOSA2分析的主要步骤

1. 构建群落和代谢模型，这里包括群落中物种对代谢通路的贡献程度；

2. 使用代谢模型计算每种微生物及其代谢得，评估在每个样本汇总每种微生物和代谢产物的影响。

3. 计算群落水平的代谢潜能得分，使用回归模型评估在不同样本中是否差异。

4. 可视化特征微生物对特定代谢的影响，并寻找关键微生物。

作者做了一份完整的介绍，参见网页：https://borenstein-lab.github.io/MIMOSA2shiny/analysis_description.html

上图展示了MIMOSA2分析所有可能的组合分析形式：使用微生物组和宏基因组数据作为输入，通过MIMOSA2可选构建两个模型，分别为：AGORA和EMBL；本次我为大家演示扩增子数据和代谢组数据如何构建模型并进行后续分析。从图中我们可以看到基于Greengenes数据库和SILVA数据库的物种注释结果，还有目前最为流行的非聚类方式得到的ASV(amplicon sequence variant)都可以用于分析。这里我们可以看到需要两个数据库：AGORAh和EMBL数据库。所以首先我们来下载数据库。

下面我们下载数据库，基于非聚类方式ASV的两种模型数据库和基于传统聚类(OTU)的两种模型数据库，注意这里基于OTU的两个模型适用于greengene和Silva数据库注释结果。download_reference_data为数据库下载命令，网站 https://borenstein-lab.github.io/MIMOSA2shiny/downloads.html# 可以手动下载数据库，但是好像我不能下载，显示拒绝访问。大家还是用下面命令好些。数据库默认下载到当前文件夹中，所以如果后使用多的话指定一个固定文件夹。四个数据库下载下来大小一共为500多M。

软件部署

安装R包

```{R warning=TRUE, include=FALSE}
rm(list=ls())

#安装包   这里我安装完成之后就将其注释掉

library(devtools)

install_github(“borenstein-lab/mimosa2”)

#如果安装失败，下载github包，进行本本地安装：

install.packages(“C:/Users/liulanlan/Desktop/mimosa2-master”, repos = NULL, type = “source”)

载入依赖关系

```{R}
suppressWarnings(suppressMessages(library(ggcorrplot)))
suppressWarnings(suppressMessages(library(igraph)))
suppressWarnings(suppressMessages(library(psych)))
suppressWarnings(suppressMessages(library(network)))
suppressWarnings(suppressMessages(library(ggplot2)))
suppressWarnings(suppressMessages(library(sna)))
suppressWarnings(suppressMessages(library(ergm)))
suppressWarnings(suppressMessages(library(ggrepel)))
suppressWarnings(suppressMessages(library("mimosa")))
suppressWarnings(suppressMessages(library(data.table)))

下载参考数据库

```{R include=FALSE}

基于非聚类方法

download_reference_data(seq_db = “Sequence variants (ASVs)”, target_db = “AGORA genomes and models”)

download_reference_data(seq_db = “Sequence variants (ASVs)”, target_db = “RefSeq/EMBL_GE Ms genomes and models”)

#基于参考数据库方法

download_reference_data(seq_db = “Greengenes 13_5 or 13_8 OTUs”, target_db = “RefSeq/EMBL_GEMs genomes and models”)

download_reference_data(seq_db = “Greengenes 13_5 or 13_8 OTUs”, target_db = “AGORA genomes and models”)

### 我们来看看MIMOSA2究竟做了什么？

![image](http://210.75.224.110/Note/WenTao/191207mimosa/fig4.png)

A. M为模拟的物质，M浓度的变化和物种1-3之间的关系如下图所示。

B. A-F为不同样本，物种1和2对于M物质具有相同的变化，这表明了这两种物质促进M物质形成，物种3消耗M物质。比较群落水平的代谢潜力(Compare total community-level metabolic potential, CMP)是最重要的指标，代表每个物种对该物质的影响得分。

C. 表示CMP进行回归并评估对给物质影响的显著性。这里给出R为评估的全部物种对该物质的解释度。

D. 给出了三个物种影响程度，这里最大的`Taxon 3`也就是所谓的支配物种，对M物质具有支配作用。

### 运行计算

这里一共有四种模式的运算：

#### 模式1：基于Greengenes OTU的结果和EMBL模型

第一种模式：使用基于聚类的结果，也就是OTU进行分析，我们在这里选择RefSeq/EMBL_GEMs genomes and models模型，其对应的数据库再之前下载得到为OTU_EMBL，这里我们指定数据库的相对路径，并指定到其中的data中。

参数文件：config_example1.txt，方便理解查查看，这也是我们唯一的输入，因为全部的参数都在这里指定好了。这里值得注意的参数是vsearch_path，我们需要下载vsearch并安装，这里的安装无论是再linux还是win下直接赋予可执行权限并放入环境变量中即可。本例子我运行在win环境，这里必须写全称：vsearch.exe，在linux下只需要写vsearch就好了。

参数文件(复制保存为txt即可调用)：

```{R eval=FALSE, include=FALSE}
file1    test_gg.txt #gg参考数据库聚类结果
file2    test_mets.txt#代谢物质结果，使用KEGG编号
ref_choices    RefSeq/EMBL_GEMs genomes and models#使用模型
file1_type    Greengenes 13_5 or 13_8 OTUs#定义输入微生物数据类型
simThreshold    0.99#相似度
data_prefix    ./database/OTU_EMBL/data/#定义参考数据库位置
vsearch_path    vsearch.exe#定义vsearch位置，此时在win下进行，所以指定为exe
rank_based    T
metType    KEGG Compound IDs#指定代谢产物ID类型
signifThreshold    1

运行第一种模式

# 下面我们选择一个例子进行运行，全部输入文件和参考数据库都备注再这里：config_example1.txt
#---运行-首先测试第一种模式：就gg数据库聚类结果使用RefSeq/EMBL_GEMs genomes and models模型分析
mimosa_results =  run_mimosa2("./config_example1.txt", make_plots = T, save_plots = T)
# unclass(mimosa_results)

模式2：基于Greengenes OTU和AGORA模型

第二种模式:使用基于聚类的结果，也就是OTU进行分析，我们在这里选择AGORA genomes and models模型，其对应的数据库再之前下载得到为OTU_AGORA，这里我们指定干数据库的相对路径，并指定到其中的data中，

参数文件：./test4.txt
```{R eval=FALSE, include=FALSE}
file1    test_gg.txt
file2    test_mets.txt
ref_choices    RefSeq/EMBL_GEMs genomes and models
file1_type    Greengenes 13_5 or 13_8 OTUs
simThreshold    0.99
data_prefix    ./database/OTU_AGORA/data/
vsearch_path    vsearch.exe
rank_based    T
metType    KEGG Compound IDs
signifThreshold    1

```{R}
#---运行=第二个测试基于OTU结果的RefSeq/EMBL_GEMs genomes and models模型分析
mimosa_results = run_mimosa2("./test4.txt", make_plots = T, save_plots = T)

模式3：基于ASV和EMBL模型

第三种模式: 使用基于非聚类的结果，也就是ASV进行分析，file1的文件类型需要更改为Sequence variants (ASVs)，我们在这里选择EMBL模型，其对应的数据库再之前下载得到为ASV_EMBL，这里我们指定干数据库的相对路径，并指定到其中的data中，
参数文件：test3.txt

```{R eval=FALSE, include=FALSE}
file1    test_seqs.txt
file2    test_mets.txt
ref_choices    EMBL_GEMs genomes and models
file1_type    Sequence variants (ASVs)
simThreshold    0.99
data_prefix    ./ASV_EMBL/data/
 vsearch_path    vsearch.exe
rank_based    T
metType    KEGG Compound IDs
signifThreshold    1

```{R}
#---运行=第二个测试基于非聚类结果
mimosa_results = run_mimosa2("./test3.txt", make_plots = T, save_plots = T)

模式4：基于ASV和AGORA模型

第四种模式：使用基于非聚类的结果，也就是ASV进行分析，file1的文件类型需要更改为Sequence variants (ASVs)，我们在这里选择AGORA genomes and models模型，其对应的数据库再之前下载得到为ASV_AGORA，这里我们指定干数据库的相对路径，并指定到其中的data中，
参数文件：test4.txt
```{R eval=FALSE, include=FALSE}
file1    test_seqs.txt
file2    test_mets.txt
ref_choices    AGORA genomes and models
file1_type    Sequence variants (ASVs)
simThreshold    0.99
data_prefix    ./database/ASV_AGORA/data/
vsearch_path    vsearch.exe
rank_based    T
metType    KEGG Compound IDs
signifThreshold    1

这里再运行出现错误，文件夹作者指定的和数据库中的对应不上，我修改了作者的两个函数，进行运行：

```{R}
#---运行=第二个测试基于非聚类结果的AGORA genomes and models模型分析
library(RColorBrewer)#调色板调用包
library("cowplot")
library(data.table)
# 主要是build_metabolic_model函数指定文件夹错误，所以我修改了该函数两处文件夹的位置
source("./build_metabolic_model-wt.R")
#载入修改后的主函数run_mimosa2-wt
source("./run_mimosa2-wt.R")
#当我调整之后，进行分析已经没有问题了
mimosa_results = run_mimosa2("./test2.txt", make_plots = T, save_plots = T)

结果说明

表格1: 微生物对代谢物的贡献表格

最重要的一列：varshare，代表了该模型汇总为微生物可以解释代谢物多少差异。仅仅包括贡献p值小于0.1代谢物。也就是说这些代谢物受到那些微生物影响，或者那些微生物影响了这些物质变化。找出对物质变化影响比较大的微生物，作为后续研究方向。

表格2：原始物种表格匹配到数据库后的新的物种表格

我们无论是基于ASV表格还是参考数据库聚类的OTU表，那种模型，都会重新匹配模拟到新的物种，这个文件展示相关信息(我们可以看到虽然参加分析的OTU表样品数量很多，但是仅展示了和代谢组共有的样本)

表格3：CMP值矩阵

CMP值矩阵，展示了每个样本中每种物种对每个物质的CMP评估值

表格4：模型统计摘要

模型统计总体简要统计：参与分析的样本数量，物种数量，参与模型的物种数量，代谢物质数量，参与模型的代谢物质数量，计算得到CMP值的代谢物质数量，匹配模型的代谢物质数量，显著的代谢物质数量···

表格5：输入参数文件，这里包含了我么的输入文件，模型选择，数据库指定等信息

图片文件展示

这里两种物质，对应两个结果。第一列图片展示的是多个样本涉及特定代谢物质的CMP值加和做的回归，R值展示为加粗数字，数值越大代表微生物群落的预测能力越准确。第二列是对应的微生物的贡献柱状图，柱状图右边的是促进该物质生成，左边的是消耗该物质。

手动输出选择：以上两种方式我们得到的默认图表如果是合并图例，如果想单独展示，这里结果中会包含图片文件我们可以手动输出，我在这里提供手动输出的文件：

#手动输出CMP拟合代谢物的图表
p = mimosa_results$CMPplots[[1]]
p+ theme_bw()

#手动导出微生物对代谢物的贡献图
p = mimosa_results$metContribPlots[[1]]
p+ theme_bw()

输出网络：

这里我们可以看到基于微生物和代谢物之间建立了一个网络，这是一个有向网络，varshare值在这里类似我们传统网络中的R值。构建网络的信息位于varshare值矩阵，varshare值作为边的权重文件，我们可以添加对微生物群落对代谢物预测能力的R值。

首先我们来构造边和节点文件，然后输出，可以选择的其他软件进行网络绘制：

# 构造边文件
networ = as.data.frame(mimosa_results$varShares)
# 构造一个向量用于承接正负关系
aa = rep(1,length(networ$VarShare))
for (i in 1:length(networ$VarShare)) {
  if (networ$VarShare[i]> 0) {
    aa[i] = 1
  }else{
    aa[i] = -1
  }
}
# aa

# 组合成边文件，添加了varshare值和方向。
edge = data.frame(compound = networ$compound,tax =networ$Species,value = networ$VarShare,direct= "directed",N_P = aa)
# head(edge)
#构造节点文件，代谢物和微生物共同组成，对分泌物的预测预测效果R值标记为一列
node = data.frame(ID= c(unique(networ$compound), unique(networ$Species) ),R =  c(unique(networ$Rsq),rep("NA",length(unique(networ$Species)))))

# head(node)
write.csv(edge,"./Gephi_edge.csv")
write.csv(node,"./Gephi_node.csv")

R 中网络可视化(可选)

#构造网络文件
G #转化为邻接矩阵，attr：填充权重
occor.r
#构造网络文件，network包提供
g # (summary(g))
#做没有权重的邻接矩阵
m #设置可视化layout
plotcord #添加表头
colnames(plotcord) = c("X1", "X2")
#添加节点名称
plotcord$elements #制作边文件
edglist edglist = as.data.frame(edglist)
# head(edglist)
#添加边的权重
set.edge.value(g,"weigt",occor.r)
edglist$weight = as.numeric(network::get.edge.attribute(g,"weigt"))
# dim(edglist)
#添加其他信息
edges # head(edges)
edges$weight = as.numeric(network::get.edge.attribute(g,"weigt"))

##这里将边权重根据正负分为两类
aaa = rep("a",length(edges$weight))
for (i in 1:length(edges$weight)) {
  if (edges$weight[i]> 0) {
    aaa[i] = "+"
  }
  if (edges$weight[i]< 0) {
    aaa[i] = "-"
  }
}
#添加到edges中
edges$wei_label = as.factor(aaa)
colnames(edges) # head(edges)

##plotcord这个表格添加注释
row.names(plotcord) = plotcord$elements
row.names(node) = node$ID
#合并之前制作的节点文件
index = merge(plotcord,node,by = "row.names",all = T)
# dim(index)
# head(index)
index$R = as.numeric(as.character(index$R))
#这里我为了突出代谢物，将全部微生物节点大小定义为代谢物的十分之一
index$R [is.na(as.numeric(as.character(index$R)))] = min(as.numeric(as.character(index$R)),na.rm = TRUE)/10
plotcord = index
#网络出图
pnet                                 data = edges, size = 0.5) +
  geom_point(aes(X1, X2,size=R),pch = 21,color = "black",fill = "red", data = plotcord) + scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  geom_text_repel(aes(X1, X2,label=elements),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

ggsave("./network.pdf", pnet, width = 12, height =8)

撰文：文涛 南京农业大学

责编：刘永鑫 中科院遗传发育所

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”