PICRUSt：16S预测宏基因组-扩增子分析锦上添花
已有 4091 次阅读
|个人分类:|系统分类:
写在前面16S分析能获得的信息比较有限，一般找到差异OTU，就很难再深入分析了。如何把差异OTU与细菌自身的基因组功能建立联系呢？很多人在这方面做出了努力。PICRUSt就是让16S扩增子分析锦上添花的工具，可以基于OTU表预测宏基因组的功能基因组成。该软件2013年发表在Nature Biotechnology上，截止17年9月11日引用1098次，由Curtis Huttenhower团队开发。最新本地版本1.1.2于17年8月30日发布。在线版最新为1.1.1于4月25日发布。已经有很多人写过PICRUSt的原理介绍，比如。但推荐的方法再好，很多小伙伴实际操作中还是有一定困难的。“宏基因组”本着只分享干货的原则，推出扩增子预测宏基因组的中文实操指南。将介绍一系列相关软件的详细使用。本部分练习所需文件位于百度网盘，链接： 密码：y33d。制作PICRUSt要求OTU表大多数人还要用de novo方法建的OTU表，而PICRUSt需要以greengene13.5版本为参考序列的OTU，这就需要自己动手了，可以用linux服务器，也可以用qiime虚拟机搞定(一般公司测序都会给分析好，没有就找他们要）。创建基于Greengene13.5的close reference OTU表的方法# 进入测序数据的文件夹
cd ~/example_PE250
# 下载13.5版参考序列，QIIME中13.8版本此软件不支持
wget ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_5_otus.tar.gz
tar xvzf gg_13_5_otus.tar.gz
# Usearch产生OTU表
./usearch10 -usearch_global temp/seqs_usearch.fa -db gg_13_5_otus/rep_set/97_otus.fasta -otutabout result/gg135_otu_table.txt -biomout result/gg135_table.biom -strand plus -id 0.97 -threads 20
# 转换文本为biom标准格式
biom convert -i result/gg135_otu_table.txt -o result/gg135_otu_table.biom --table-type=&OTU table& --to-json
在线版操作Picurst在线服务器 准备数据和上传先本地生成close reference OTU(gg135_otu_table.txt 或 gg135_table.biom)，或用在线测序数据 先从左侧下部get data — upload file — choose local file — type为Auto-dected或picurst；再选择start，完成后点close拷贝数标准化biom格式选择第一项，文本格式选择第二项，提交即可(biom格式报错；文本的也报错。可能是服务器暂时不可用，换个时间一般会好)。预测宏基因组左侧工具栏 PICRUSt — Predict Metagenome,选择刚生成的Normalized By Copy Number on data xx，点击execute即可。按功能类别分类汇总左侧工具栏 PICRUSt — Categorize by function, 输入选择Predict Metagenome on data xx, 级别使用默认的3，输出默认为BIOM，建议改为txt，点击Execute。就会生成一个表下载，这就是转换为KO的表，和OTU类似，但是功能描述。Linux本地版使用安装picrustPICRUST 依赖关系包括: python 2.7, PyCogent 1.5.3, biom 2.1.4, numpy 1.5.1, 这些都是qiime1.9.1所依赖的，因此只要安装过qiime这都不是问题了(all installed in qiime)，详见 &# 下载并安装，管理员权限非常顺利，也方便所有人使用
git clone git://github.com/picrust/picrust.git picrust
cd picrust
sudo pip install .
# 下载数据库
download_picrust_files.py
运行picrust# 按拷备数标准化OTU表， -c指定数据库
normalize_by_copy_number.py -i result/gg135_table.biom -o normalized_otus.biom # 用usearch的biom结果会报错
normalize_by_copy_number.py -i result/gg135_otu_table.biom -o result/normalized_otus.biom -c picrust/picrust/data/16S_13_5_precalculated.tab.gz
# 预测宏基因组
predict_metagenomes.py -i result/normalized_otus.biom -o result/metagenome_predictions.biom -c picrust/picrust/data/ko_13_5_precalculated.tab.gz
# 转换结果为txt查看
biom convert -i result/metagenome_predictions.biom -o result/metagenome_predictions.txt --table-type=&OTU table& --to-tsv
# 调整为标准表格
sed -i '/# Const/d;s/#OTU //g' result/metagenome_predictions.txt
# 预览KO表，和OTU表类似，只是OTU变为了KO
less result/metagenome_predictions.txt
# 统计KO数据6910个
wc -l result/metagenome_predictions.txt
# 按功能级别分类汇总
# -c指输出类型，有KEGG_Pathways, COG_Category，RFAM三种，-l是级别，分4级
categorize_by_function.py -i result/metagenome_predictions.biom -c KEGG_Pathways -l 3 -o result/metagenome_predictions.L3.biom
# 转换结果为txt查看
biom convert -i result/metagenome_predictions.L3.biom -o result/metagenome_predictions.L3.txt --table-type=&OTU table& --to-tsv
# 调整为标准表格
sed -i '/# Const/d;s/#OTU //g' result/metagenome_predictions.L3.txt
# 预览KO表，和OTU表类似，只是OTU变为了KO
less -S result/metagenome_predictions.L3.txt
# 统计KO数据L3级别有329个
wc -l result/metagenome_predictions.L3.txt
# 查找指定通路的丰度情况(可选，用于具体查看某个KO)
metagenome_contributions.py -i result/normalized_otus.biom -l K0,K00004 -o result/ko_metagenome_contributions.tab
最终结果数据格式展示KO表的部分结果格式ID & & &KO4 & & KO6 & & KO5 & & KO3 & & WT8 & & KO9 & &
K0 & & 0.0 & & 0.0 & & 0.0 & & 0.0 & & 0.0 & &
K0 & & 0.0 & & 0.0 & & 0.0 & & 0.0 & & 0.0 & &
K0 & &29.0 & &31.0 & &24.0 & &2.0 & & 26.0 &
K0 & & 0.0 & & 1.0 & & 0.0 & & 1.0 & & 0.0 & &
K0.0 31.0 85.0 90576.0
K0 & & 11.0 & &4.0 & & 4.0 & & 0.0 & & 2.0 & &
K0.0 &26.0 &34.0 &2795.0
展示KO level3级的部分结果格式ID & & &KO4 & & KO6 & & KO5 & & KO3 & & WT8 & & KO9 & &
1,1,1-Trichloro-2,2-bis(4-chlorophenyl)ethane (DDT) degr
ABC transporters & & & & & & &
Adherens junction & & & 142.0 & 154.0 & 123.0 & 119.0 &
Adipocytokine signaling pathway 34.0 76865.0
African trypanosomiasis 42.0 64.0
Alanine, aspartate and glutamate metabolism & &
Aldosterone-regulated sodium reabsorption & & & 85.0 & &
Alzheimer's disease & & 46.0 64.0
Amino acid metabolism &
& & & & & & & &
Amino acid related enzymes & & & & & & &
Amino sugar and nucleotide sugar metabolism & &
Aminoacyl-tRNA biosynthesis & &
Aminobenzoate degradation & & &
宏基因组功能分类(KO)表统计与可视化我们现在获得了KO表或L3级别的功能表，它们同OTU表类似。接下来，可以用官方推荐的STAMP、HUManN、LEfSe和GraPhlAn进行分析。也可以采用之前我介绍的edgeR或DEseq等R包，对KO表按分组进行统计，找到差异用热图或柱状图展示即可。详见
&我们可以使用上文中的代码，对KO像OTU一样分析，看样品间的相关性聚类、差异KO功能分析并热图可视化。KO与OTU相关分析比较图1. 原始基于OTU相对丰度样品相关性的热图，WT， OE，KO(这里KO是组名，表示基因敲除)明显分为三组图2. 基于宏基因组功能(KO)相对丰度样品相关性热图，OE还是明显的一组，KO和OE的差别非常不明显，而且分成了两个混杂的类群。一种可能是这个基因的突变体引起了菌群结构变化，但在功能基因层面真的和WT差异不大。但更可能的原因是：PICRUSt在本研究中预测的宏基因组数据不够准确全面，偏差到了无法区分WT与突变体，因为原始注释中有很多OTU并不在greengene数据中可以匹配，这些PICRUSt是无法分析的。差异功能KO统计与热图展示但OE与WT差异还是明显的，我用上文热图中的代码套用在KO表进一步统计差异KO，并展示。图3. 差异KO(Level3, P & 5e-7)详细的绘图代码参考上文热图，也可以在百度云result目录中下载DAKO_egr.r文件，读取OTU替换为KO的表，绘图主要变化只有两处如下：# 筛选差异KO的阈值，用默认的0.05会有168个差异，热图展示会看不清字。调整至5e-7减少至50个左右只是为方便展示
de_lrt = decideTestsDGE(lrt, adjust.method=&fdr&, p.value=0.0000005)
# 热图绘制，注意KO的名称较长，调整margins第二个值为17才显示完整
heatmap.2(sub_norm, scale=&row&, Colv=TRUE, Rowv=TRUE,dendrogram=&both&, col=rev(colorRampPalette(brewer.pal(11, &RdYlGn&))(256)), cexCol=1,keysize=1,density.info=&none&,main=NULL,trace=&none&,margins = c(3, 17))
关于差异KO的展示，还常用泡泡图、柱状图结合不同分类级分组展示，这都是绘图的范畴了，我会在将来也会进行详细讲解的。想要自己画的，可以先学习陈同博士的(见公众号菜单-文章目录-R语言绘图部分)。接触过泛基因组人都知道，即使16S相同，由于细菌序列的大量水平转移，功能基因也是相差极大。所以本软件的结果，仅供描述。富人做宏基因组，穷人16S预测功能锦上添花。ReferencePredictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Langille, M. G.I.; Zaneveld, J.; Caporaso, J. G.; McDonald, D.; Knights, D.; a Reyes, J.; Clemente, J. C.; Burkepile, D. E.; Vega Thurber, R. L.; Knight, R.; Beiko, R. G.; and Huttenhower, C. Nature Biotechnology, 1-10. 8 2013. 原文 主页 1.1.1服务器 写在后面为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内五十位PI，五百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读学习解决问题思路，仍末解决推荐生信技能树-微生物组版块() 发贴，并转发链接入群，问题及解答方便检索，造福后人。学习16S扩增子、宏基因组思路和分析实战，快关注“宏基因组”，干货第一时间推送。点击阅读原文，跳转最新文章目录阅读
转载本文请联系原作者获取授权，同时请注明本文来自刘永鑫科学网博客。链接地址：
上一篇：下一篇：
当前推荐数：0
推荐到博客首页
评论 ( 个评论)
扫一扫，分享此博文
作者的精选博文
作者的其他最新博文
热门博文导读
Powered by
Copyright &我来说两句
点赞& 条评论&&&&&人次参与
查看更多评论 > >
生物通微信公众号
生物通新浪微博
热搜：||||
知名企业招聘
版权所有 生物通
Copyright©
eBiotrade.com, All Rights Reserved
联系信箱：16S预测细菌组表型-bugbase：革兰氏阴阳、生物膜、致病力、移...
已有 1149 次阅读
|个人分类:|系统分类:
基于16S OTU表预测细菌表型数据库，同时可进行组间差异分析BugBase Predicts Organism Level Microbiome PhenotypesTonya Ward, Jake Larson, Jeremy Meulemans, Ben Hillmann, Joshua Lynch, Dimitri Sidiropoulos, John Spear, Greg Caporaso, Ran Blekhman, Rob Knight, Ryan Fink, Dan KnightsbioRxiv 133462; doi: 此文今天5月2日发布在预印本杂志bioRxiv上面，还没有正式发表。短短几个月己被引用3次，阅读2138次，全文下载789次。待正式发表发，必将成为本领域功能预测的又一神器。居然有33条Tweets的转发和留言，看看大家的评价：官网：老司机直接Parse Data上传数据开分析，新人一定要先点下方Documentation，仔细阅读使用说明。使用说明BugBase是一款分析微生物组样品表型的工具，此网站可以基于OTU表和Mapping files，预测大量信息和比较，包括以下七方面：革兰氏阳性 Gram Positive革兰氏阴性 Gram Negative生物膜形成 Biofilm Forming致病潜力 Pathogenic Potential移动元件含量 Mobile Element Containing氧的利用 Oxygen Utilizing氧化胁迫耐受 Oxidative Stress Tolerant同时BugBase还会按组进行分类统计与可视化。输入文件要求OTU表BIOM 1.0格式16S以GreenGenes 13.5 为Reference宏基因组以IMG为参考小于15 mb (本地版无限制)Maping File制表符分隔第一行必须以#SampleID起始第一行全为列标题第一列必须为SampleID只允许使用字母、数字、下划线和连字符不允许包含空格，逗号、引号、括号不要包含机密信息其实就是符合QIIME的标准即可，如下图示例分组列名：指定分组信息，如上面的BODY_SITE，可产生如下结果：还有很多其它功能，大家自己读帮助文档吧。测试数据分析可以点主页右边的Downloads，下载OTU表和Mapping file，再到 Parse Data中上传这两个文件，Column Header添“BODY_SITE”，点Parse Data很快生成了结果链接，打包下载所有结果。本地化bugbase：安装和运行对一直用服务器的小伙伴，肯定希望数据库有本地版，可以整合入自己的分析pipeline，随时搞，随便搞，一条命令搞定不求人。更不怕网站无法访问或无人维护。这个数据库还真可以本地方，满足muggle和geek的所有需求，安装代码如下：源代码和数据库： # 下载安装程序
cd ~/software
wget https://github.com/knights-lab/BugBase/archive/master.zip
mv master.zip BugBase.zip
unzip BugBase.zip
mv BugBase-master/ BugBase
# 此程序运行必须定义下面环境变量，根据实际目录修改
export BUGBASE_PATH=/home/user/software/BugBase
export PATH=$PATH:/home/user/software/BugBase/bin
安装依赖R包并显示帮助run.bugbase.r -h &# 安装了所有依赖包
# 以上R包如果已经安装，此步可跳过
# 每次运行都会重复安装10多个包近半小时
显示程序的使用帮助如下：内容较多，只显示前四个参数，完整信息见附录1。所有依赖包列表和安装代码见附录2Usage: /mnt/bai/yongxin/software/BugBase/bin/run.bugbase.r [options]
-i OTUTABLE, --otutable=OTUTABLE
otu table to plot [default NULL]
-c MAPCOLUMN, --mapcolumn=MAPCOLUMN
column of mapping file to plot [default NULL]
-m MAPPINGFILE, --mappingfile=MAPPINGFILE
mapping file to plot [default NULL]
-o OUTPUT, --output=OUTPUT
output directory [default .]
运行演示数据# 运行演示数据
run.bugbase.r -i $BUGBASE_PATH/doc/data/HMP_s15.txt -m $BUGBASE_PATH/doc/data/HMP_map.txt -c HMPBODYSUBSITE -o output
运行中会显示运行内容如下[1] &Loading Inputs...&
[1] &16S copy number normalizing OTU table...&
[1] &Predicting phenotypes...&
[1] &313 OTUs from the input table matched the 203452 available database OTUs&
[1] &Plotting thresholds...&
[1] &Plotting predictions...&
[1] &Plotting OTU contributions...&
[1] &BugBase analysis complete&
输入参数文件详解OTU表 -i $BUGBASE_PATH/doc/data/HMP_s15.txt```less -S $BUGBASE_PATH/doc/data/HMP_s15.txtConstructed from biom fileOTU ID SRS280 & & & &SRS993.0 & & 0.0 & & 0.0 & & 0.0 & & 0.0 & & 0.0.0 & &0.0 & & 0.0 & & 16.0 & &1.0 & & 12.0.0 & & 0.0 & & 0.0 & & 10.0 & &3.0 & & 3.0 2. mapping file `-m` $BUGBASE_PATH/doc/data/HMP_map.txt
sed -i ‘s/^M/\n/g’ BUGBASE_PATH/doc/data/HMP_map.txtless -S $BUGBASE_PATH/doc/data/HMP_map.txtSampleID & & & BarcodeSequence SAMP_COLLECT_DEVICE & & TITLE & RUN_PREFIX & & &AGE & & COMMON_NAME & & BODY_SITESRS681 & & & &TCAGAATACGTTC & stool_specimen_collection_kit & HMP_production_phase_2715 & & & SRR040712SRS389 & & & &TCAGAGTCGAC & & stool_specimen_collection_kit & HMP_production_phase_199 & & & &SRR040857SRS892 & & & &TCAGAACCGGATAC &stool_specimen_collection_kit & HMP_production_phase_1681 & & & SRR040994 3. 分组列名 `-c` HMPBODYSUBSITE
### 输出文件详解
输出结果在线分析直接生成压缩包下载。本地在Ouput目录中有四个目录：
- `normalized_otus`目录中有`16s_normalized_otus.txt`文件，为标准化的OTU表
less -S output/normalized_otus/16s_normalized_otus.txt文件如下： & &SRS681 & & & &SRS389 & & & &SRS892
& & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 2 & & & 11.75 & 3 & & & 3 & & & 7 & & & 0 & & & 0 & & & 3 & & & 0 & & & 8 & & & 3 & & & 27 & & &0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0 & & & 0- `otu_contributions`目录中
主要有9种表型或功能预测结果表`contributing_otus.txt`，和9种表型按实验组比较的结果堆叠柱状图和物种颜色方案图例PDF版
less -S output/otu_contributions/contributing_otus.txt预测结果前4行内容如下 & &Aerobic Anaerobic & & & Contains_Mobile_Elements & & & &Facultatively_Anaerobic Forms_Biofilms &Gram_Negative & Gram_Positive & Potentially_Pathogenic &Stress_Tolerant
963239 &FALSE & TRUE & &TRUE & &FALSE & FALSE & FALSE & TRUE & &FALSE & TRUE4431292 FALSE & TRUE & &TRUE & &FALSE & FALSE & FALSE & TRUE & &TRUE & &TRUE4480529 FALSE & TRUE & &TRUE & &FALSE & FALSE & FALSE & TRUE & &FALSE & TRUE
图0. 注释物种门水平图例
图1. 各组需氧菌相对丰度
图2 各组厌氧菌相对丰度
图3. 革兰氏阴性菌相对丰度
一共有9种，不再一一列举
- `predicted_phenotypes`目录中
主要有9种表型或功能预测结果表`predictions.txt`，和9种表型按实验组比较箱线图，和相关组间统计信息。
less -S output/predicted_phenotypes/predictions.txt & &Aerobic Anaerobic & & & Contains_Mobile_Elements & & & &Facultatively_Anaerobic Forms_Biofilms &Gram_Negative & Gram_Posi
SRS681 & & & &0 & & & 0.892 & & & 0.441 & & & 0 & & & 0 & & & 0.405 & & & 0.2178269SRS389 & & & &0 & & & 1 & & & 0.371 & & & 0 & & & 0 & & & 0.15 & & & &0.85 & & & &0SRS892 & & & &0 & & & 0.02 & & & &0.425 & & & 0.88292 & & 0 & & & 0.566
图4. 箱线图展示9种表型中的移动元件含量
- `thresholds`目录中包含分析使用的阈值`thresholds_used.txt`和不同阈值下的结果`variances.txt`数据，以及9表型在不同阈值下相对丰度变化
图5. 折线图展示9种表型中的生物膜形成菌在不同阈值下相对丰度变化
### Reference
1. BugBase Predicts Organism Level Microbiome Phenotypes
Tonya Ward, Jake Larson, Jeremy Meulemans, Ben Hillmann, Joshua Lynch, Dimitri Sidiropoulos, John Spear, Greg Caporaso, Ran Blekhman, Rob Knight, Ryan Fink, Dan Knights
bioRxiv 133462; doi: https://doi.org/10.
#### 附录1. 程序参数详解
Usage: /mnt/bai/yongxin/software/BugBase/bin/run.bugbase.r [options]Options:-v, —verbosePrint extra output [default]-i OTUTABLE, —otutable=OTUTABLEotu table to plot [default NULL]-c MAPCOLUMN, —mapcolumn=MAPCOLUMNcolumn of mapping file to plot [default NULL]-m MAPPINGFILE, —mappingfile=MAPPINGFILEmapping file to plot [default NULL]-o OUTPUT, —output=OUTPUToutput directory [default .]-t TAXALEVEL, —taxalevel=TAXALEVELtaxa level to plot otu contributions by, default is & & & & &2 (phylum) [default NULL]-p PHENOTYPE, —phenotype=PHENOTYPEspecific traits (phenotypes) to predict, separated by & & & & &commas, no spaces [default NULL]-x, —predictonly output the prediction table, do not make plots & & & & &[default FALSE]-T THRESHOLD, —threshold=THRESHOLDthreshold to use, must be between 0 and 1 & & & & &[default NULL]-g GROUPS, —groups=GROUPStreatment groups of samples, separated by commas, no spaces & & & & &[default NULL]-u USERTABLE, —usertable=USERTABLEuser define trait table, absolute file path required [ & & & & &default NULL]-z, —continuousplot continuous data [default FALSE]-k, —kegguse kegg pathway table [default FALSE]-C, —covuse coefficient of variance instead of variance [default FALSE]-l, —clr_transuse centered log-ratio transformation instead of relative abundance [default FALSE]-a, —allplot all samples without a mapping file (this outputs no & & & & &statistics) [default FALSE]-w, —shotgunData is metagenomic shotgun data & & & & &(picked against RefSeq database) [default FALSE]-h, —helpShow this help message and exit
#### 附录2. 安装所有依赖R包并测试
运行RR安装13个依赖包，环境为Ubuntu16.04 + R3.4.1，祝你成功source(““)biocLite(c(“optparse”, “beeswarm”, “RColorBrewer”, “reshape2”, “plyr”, “grid”, “gridExtra”, “ggplot2”, “RJSONIO”, “biom”, “Matrix”, “labeling”, “digest”))测试是否安装成功library(“optparse”)library(“beeswarm”)library(“RColorBrewer”)library(“reshape2”)library(“plyr”)library(“grid”)library(“gridExtra”)library(“ggplot2”)library(“RJSONIO”)library(“biom”)library(“Matrix”)library(“labeling”)library(“digest”)q()
#### 附录3. 常见错误
1. otu表格式错误
[1] “Loading Inputs…”Error in load.inputs(otu_table, map, mapcolumn, groups) :Error: otu table must be either .txt or .biom (json)Execution halted解决方法：要求扩展名符合要求，如文本表必须为.txt结尾，更正即可
2. 输出目录己存在
Error: Output directory already existsExecution halted```程序防止覆盖结果，手动删除输出目录即可，或换个不存在的输出目录附录4. 革兰氏细菌分类解决思路百度搜索“革兰氏 分类 数据库”没有结果 文章中提到了细菌分类，但作者说是自己整理的数据库，且不可以给我google搜索“Gram stain database”主要相关结果：
biostar和我同样的问题，大家给了如下多条建议 临床最终参考细菌数据库，其中革兰氏阴、阳；只有一些名字，多个网页，需要整理 此文给出了wiki和分类手册的链接 有列表链接，比较好整理
1957年分类手册，太老了吧 NCBI有文件lproks_0.txt，染色为菌的一种表型但完全没能快速找到针对性的整理数据库或软件。后来在在群中提问：“有基于16S鉴定革兰氏分类的数据库吗？”，两个群里都有人回答了bugbase。使用发现专门解决我的问题，只是没有发表的软件，怪不得查不到，但也显示了大家知识的力量。一句话解决了可能几天都无法解决的问题。猜你喜欢一文读懂：
必备技能：
扩增子分析：
&科研团队经验： &
&系列教程：
& 生物科普
& 写在后面为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外100+ PI，800+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”点击阅读原文，跳转最新文章目录阅读
转载本文请联系原作者获取授权，同时请注明本文来自刘永鑫科学网博客。链接地址：
上一篇：下一篇：
当前推荐数：0
推荐到博客首页
评论 ( 个评论)
扫一扫，分享此博文
作者的精选博文
作者的其他最新博文
热门博文导读
Powered by
Copyright &生物信息学-宏基因组学方向
PICRUSt：预测宏基因组功能—16S扩增子分析锦上添花
16S分析能获得的信息比较有限，一般找到差异OTU，就很难再深入分析了。
如何把差异OTU与细菌自身的基因组功能建立联系呢？很多人在这方面做出了努力。
PICRUSt就是让16S扩增子分析锦上添花的工具，可以基于OTU表预测宏基因组的功能基因组成。该软件2013年发表在Nature Biotechnology上，截止17年9月11日引用1098次，由Curtis Huttenhower团队开发。最新本地版本1.1.2于17年8月30日发布。在线版最新为1.1.1于4月25日发布。
已经有很多人写过PICRUSt的原理介绍，比如。
但推荐的方法再好，很多小伙伴实际操作中还是有一定困难的。“宏基因组”本着只分享干货的原则，推出扩增子预测宏基因组的中文实操指南。将介绍一系列相关软件的详细使用。
本部分练习所需文件位于百度网盘，链接： 密码：y33d。
制作PICRUSt要求OTU表
大多数人还要用de novo方法建的OTU表，而PICRUSt需要以greengene13.5版本为参考序列的OTU，这就需要自己动手了，可以用linux服务器，也可以用qiime虚拟机搞定(一般公司测序都会给分析好，没有就找他们要）。
创建基于Greengene13.5的close reference OTU表的方法
# 进入测序数据的文件夹
cd ~/example_PE250
# 下载13.5版参考序列，QIIME中13.8版本此软件不支持
wget ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_5_otus.tar.gz
tar xvzf gg_13_5_otus.tar.gz
# Usearch产生OTU表
./usearch10 -usearch_global temp/seqs_usearch.fa -db gg_13_5_otus/rep_set/97_otus.fasta -otutabout result/gg135_otu_table.txt -biomout result/gg135_table.biom -strand plus -id 0.97 -threads 20
# 转换文本为biom标准格式
biom convert -i result/gg135_otu_table.txt -o result/gg135_otu_table.biom --table-type="OTU table" --to-json
在线版操作
Picurst在线服务器
1. 准备数据和上传
先本地生成close reference OTU(gg135_otu_table.txt 或 gg135_table.biom)，或用在线测序数据
先从左侧下部get data – upload file – choose local file – type为Auto-dected或picurst；再选择start，完成后点close
2. 拷贝数标准化
biom格式选择第一项，文本格式选择第二项，提交即可(biom格式报错；文本的也报错。可能是服务器暂时不可用，换个时间一般会好)。
3. 预测宏基因组
左侧工具栏 PICRUSt – Predict Metagenome,
选择刚生成的Normalized By Copy Number on data xx，点击execute即可。
4. 按功能类别分类汇总
左侧工具栏 PICRUSt – Categorize by function, 输入选择Predict Metagenome on data xx, 级别使用默认的3，输出默认为BIOM，建议改为txt，点击Execute。就会生成一个表下载，这就是转换为KO的表，和OTU类似，但是功能描述。
Linux本地版使用
安装picrust
依赖关系包括: python 2.7, PyCogent 1.5.3, biom 2.1.4, numpy 1.5.1, 这些都是qiime1.9.1所依赖的，因此只要安装过qiime这都不是问题了(all installed in qiime)，详见
# 下载并安装，管理员权限非常顺利，也方便所有人使用
git clone git:
cd picrust
sudo pip install .
# 下载数据库
download_picrust_files.py
运行picrust
normalize_by_copy_number.py -i result/gg135_table.biom -o normalized_otus.biom
normalize_by_copy_number.py -i result/gg135_otu_table.biom -o result/normalized_otus.biom -c picrust/picrust/data/16S_13_5_precalculated.tab.gz
predict_metagenomes.py -i result/normalized_otus.biom -o result/metagenome_predictions.biom -c picrust/picrust/data/ko_13_5_precalculated.tab.gz
biom convert -i result/metagenome_predictions.biom -o result/metagenome_predictions.txt
sed -i '/# Const/d;s/#OTU //g' result/metagenome_predictions.txt
less result/metagenome_predictions.txt
wc -l result/metagenome_predictions.txt
categorize_by_function.py -i result/metagenome_predictions.biom -c KEGG_Pathways -l 3 -o result/metagenome_predictions.L3.biom
biom convert -i result/metagenome_predictions.L3.biom -o result/metagenome_predictions.L3.txt
sed -i '/# Const/d;s/#OTU //g' result/metagenome_predictions.L3.txt
less -S result/metagenome_predictions.L3.txt
wc -l result/metagenome_predictions.L3.txt
metagenome_contributions.py -i result/normalized_otus.biom -l K0,K00004 -o result/ko_metagenome_contributions.tab
最终结果数据格式
展示KO表的部分结果格式
K0.0 31.0 85.0 90576.0
展示KO level3级的部分结果格式
1,1,1-Trichloro-2,2-bis(4-chlorophenyl)ethane (DDT) degr
ABC transporters
Adherens junction
Adipocytokine signaling pathway 90482.0 65434.0 76865.0
African trypanosomiasis 15284.0 10342.0 13778.0 13364.0
Alanine, aspartate and glutamate metabolism
Aldosterone-regulated sodium reabsorption
Alzheimer's disease
Amino acid metabolism
Amino acid related enzymes
Amino sugar and nucleotide sugar metabolism
Aminoacyl-tRNA biosynthesis
Aminobenzoate degradation
宏基因组功能分类(KO)表统计与可视化
我们现在获得了KO表或L3级别的功能表，它们同OTU表类似。接下来，可以用官方推荐的STAMP、HUManN、LEfSe和GraPhlAn进行分析。
也可以采用之前我介绍的edgeR或DEseq等R包，对KO表按分组进行统计，找到差异用热图或柱状图展示即可。详见
我们可以使用上文中的代码，对KO像OTU一样分析，看样品间的相关性聚类、差异KO功能分析并热图可视化。
KO与OTU相关分析比较
图1. 原始基于OTU相对丰度样品相关性的热图，WT， OE，KO(这里KO是组名，表示基因敲除)明显分为三组
图2. 基于宏基因组功能(KO)相对丰度样品相关性热图，OE还是明显的一组，KO和OE的差别非常不明显，而且分成了两个混杂的类群。
一种可能是这个基因的突变体引起了菌群结构变化，但在功能基因层面真的和WT差异不大。但更可能的原因是：PICRUSt在本研究中预测的宏基因组数据不够准确全面，偏差到了无法区分WT与突变体，因为原始注释中有很多OTU并不在greengene数据中可以匹配，这些PICRUSt是无法分析的。
差异功能KO统计与热图展示
但OE与WT差异还是明显的，我用上文热图中的代码套用在KO表进一步统计差异KO，并展示。
图3. 差异KO(Level3, P & 5e-7)
详细的绘图代码参考上文热图，也可以在百度云result目录中下载DAKO_egr.r文件，读取OTU替换为KO的表，绘图主要变化只有两处如下：
de_lrt = decideTestsDGE(lrt, adjust.method="fdr", p.value=0.0000005)
heatmap.2(sub_norm, scale="row", Colv=TRUE, Rowv=TRUE,dendrogram="both", col=rev(colorRampPalette(brewer.pal(11, "RdYlGn"))(256)), cexCol=1,keysize=1,density.info="none",main=NULL,trace="none",margins = c(3, 17))
关于差异KO的展示，还常用泡泡图、柱状图结合不同分类级分组展示，这都是绘图的范畴了，我会在将来也会进行详细讲解的。想要自己画的，可以先学习陈同博士的(见公众号菜单-文章目录-R语言绘图部分)。
接触过泛基因组人都知道，即使16S相同，由于细菌序列的大量水平转移，功能基因也是相差极大。所以本软件的结果，仅供描述。富人做宏基因组，穷人16S预测功能锦上添花。
Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Langille, M. G.I.; Zaneveld, J.; Caporaso, J. G.; McDonald, D.; Knights, D.; a Reyes, J.; Clemente, J. C.; Burkepile, D. E.; Vega Thurber, R. L.; Knight, R.; Beiko, R. G.; and Huttenhower, C. Nature Biotechnology, 1-10. 8 2013.
1.1.1服务器
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内五十位PI，五百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读学习解决问题思路，仍末解决推荐生信技能树-微生物组版块() 发贴，并转发链接入群，问题及解答方便检索，造福后人。
学习16S扩增子、宏基因组思路和分析实战，快关注“宏基因组”，干货第一时间推送。
点击阅读原文，跳转最新文章目录阅读
没有更多推荐了，