seo网站关键词优化多少钱,简单网页模板图片,m导航网站如何做淘宝客,wordpress 多站 列表MIMOSA2#xff1a;基于微生物组和代谢组数据的整合分析MIMOSA2 升级自MIMOSA1。是 Borenstein 实验室(http://borensteinlab.com/ , 专注宏基因组系统 生物学)最新开发的工具。用于微生物群落和代谢组的整合分析#xff0c;寻找微生物和代谢产物之间的关系。先前Borenstein … 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通过微生物群落特征和代谢组特征关联回答以下几个问题代谢组数据是否和微生物组数据紧密相关微生物群落差异是否可以解释代谢差异何种微生物或者基因造成了代谢差异或者贡献了代谢差异。MIMOSA2的工作原理工作原理首先得到微生物组数据使用参考数据库构建代谢模型计算不同微生物对不同代谢物的潜在代谢潜能指标(CMP)对每种代谢产物全部潜在贡献微生物CMP值进行汇总并做回归分析评估微生物群落对指定代谢物的预测能力最后筛选对特定代谢物影响加大的微生物作为关键微生物。MIMOSA2分析的主要步骤构建群落和代谢模型这里包括群落中物种对代谢通路的贡献程度使用代谢模型计算每种微生物及其代谢得评估在每个样本汇总每种微生物和代谢产物的影响。计算群落水平的代谢潜能得分使用回归模型评估在不同样本中是否差异。可视化特征微生物对特定代谢的影响并寻找关键微生物。作者做了一份完整的介绍参见网页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 warningTRUE, includeFALSE}rm(listls())#安装包 这里我安装完成之后就将其注释掉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 includeFALSE}基于非聚类方法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究竟做了什么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 evalFALSE, includeFALSE}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下进行所以指定为exerank_based TmetType 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 evalFALSE, includeFALSE}file1 test_gg.txtfile2 test_mets.txtref_choices RefSeq/EMBL_GEMs genomes and modelsfile1_type Greengenes 13_5 or 13_8 OTUssimThreshold 0.99data_prefix ./database/OTU_AGORA/data/vsearch_path vsearch.exerank_based TmetType KEGG Compound IDssignifThreshold 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 evalFALSE, includeFALSE}file1 test_seqs.txtfile2 test_mets.txtref_choices EMBL_GEMs genomes and modelsfile1_type Sequence variants (ASVs)simThreshold 0.99data_prefix ./ASV_EMBL/data/ vsearch_path vsearch.exerank_based TmetType KEGG Compound IDssignifThreshold 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 evalFALSE, includeFALSE}file1 test_seqs.txtfile2 test_mets.txtref_choices AGORA genomes and modelsfile1_type Sequence variants (ASVs)simThreshold 0.99data_prefix ./database/ASV_AGORA/data/vsearch_path vsearch.exerank_based TmetType KEGG Compound IDssignifThreshold 1这里再运行出现错误文件夹作者指定的和数据库中的对应不上我修改了作者的两个函数进行运行{R}#---运行第二个测试基于非聚类结果的AGORA genomes and models模型分析library(RColorBrewer)#调色板调用包library(cowplot)library(data.table)# 主要是build_metabolic_model函数指定文件夹错误所以我修改了该函数两处文件夹的位置source(./build_metabolic_model-wt.R)#载入修改后的主函数run_mimosa2-wtsource(./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表样品数量很多但是仅展示了和代谢组共有的样本)表格3CMP值矩阵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 #设置可视化layoutplotcord #添加表头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$elementsrow.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)/10plotcord index#网络出图pnet data edges, size 0.5) geom_point(aes(X1, X2,sizeR),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,labelelements),size4, 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())pnetggsave(./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扩增子、宏基因组科研思路和分析实战关注“宏基因组”
