本文为R包microbiomeViz作者李陈浩原创，博文链接见文末。作者授权发布于本平台，刘永鑫对本文进行测试，有修改。

为啥写这个

平日经常会分析shotgun宏基因组的数据，我们的pipeline使用MetaPhlAn，Kraken等profiler。这种数据经常会产生一个表格，如下

download.file("https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/output/SRS014459-Stool_profile.txt", 'SRS014459-Stool_profile.txt')
knitr::kable(head(read.table('SRS014459-Stool_profile.txt')))

第一列是分类信息注释，第二列是相对丰度（百分比）。在做这种图可视化方面，目前个人见过最强大的是GraPhlAn:

官网上相关的教程很详细，但是问题是，这个完全封闭的python程序，想要hack，还真的是挺难得。Krona可能是另一个选择，但是同样还是会有同样的问题。最近发布的R包Metacoder，画出的图个人真心不是很喜欢：

跟Y叔讨论了一下用ggtree实现像GraPhlAn那样图的可能性，得到了肯定的答复，于是开始自己造轮子。

MicrobiomeViz–千里之行，始于足下
其实可以写一个简单的函数，但是还是想做一个拓展性更强的东西，所以就有了这个包（不断完善中）： https://github.com/lch14forever/microbiomeViz

使用实战

让我们产生lefse调用graphlan绘制的物种树标记差异物种的Cladogram

输入数据为metaphlan2结果合并的矩阵。如何生成详见：MetaPhlAn2一条命令获得宏基因组物种组成

ID      BM_SRS013506    BM_SRS015374    BM_SRS015646    BM_SRS017687    BM_SRS019221    BM_SRS019329    BM_SRS020336    BM_SRS022145    BM_SRS022532    
k__Bacteria     100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   100.0   
k__Bacteria|p__Actinobacteria   1.33609 2.90435 0.45117 6.7964  14.08966        2.30709 7.30108 0.53534 3.57207 8.47622 7.07037 17.30722        30.62601
k__Bacteria|p__Actinobacteria|c__Actinobacteria 1.33609 2.90435 0.45117 6.7964  14.08966        2.30709 7.30108 0.53534 3.57207 8.47622 7.07037 17.30722

包安装和加载

# microbiomeViz需要 R 3.5 以上，依赖包安装
source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")
devtools::install_github("lch14forever/microbiomeViz")
library(microbiomeViz)

物种相对丰对矩阵绘制物种树

# 加载测试数据
df <- read.table("http://bailab.genetics.ac.cn/markdown/R/microbiomeViz/merged_abundance_table.txt", head=TRUE, stringsAsFactors = FALSE)

## 计算均值用于呈现结点大小
dat <- data.frame(V1=df[,1], V2=rowMeans(df[,-1]), stringsAsFactors = FALSE)

# 用物种和丰度生成树骨架
tr <- parseMetaphlanTSV(dat, node.size.offset=2, node.size.scale=0.8)
p <- tree.backbone(tr, size=0.5)
p

差异物种注释

# 读取需要颜色标注的差异物种列表，本质上是两列和颜色对应表
lefse_lists = data.frame(node=c('s__Haemophilus_parainfluenzae','p__Proteobacteria',
                                'f__Veillonellaceae','o__Selenomonadales',
                                'c__Negativicutes', 's__Streptococcus_parasanguinis',

                                'p__Firmicutes','f__Streptococcaceae',
                                'g__Streptococcus','o__Lactobacillales',
                                'c__Bacilli','s__Streptococcus_mitis'),
                         color=c(rep('darkgreen',6), rep('red','6')),
                         stringsAsFactors = FALSE
)


# 注释树

p <- clade.anno(p, lefse_lists, alpha=0.3)
p

简单几行代码，美图大功告成。

Reference

http://lchblogs.netlify.com/post/2018-01-18-r-metagenomeviz/

http://lchblogs.netlify.com/post/2018-04-20-r-microbiomeviz_example/

猜你喜欢

• 10000+：肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板 Shell  R Perl
• 生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

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

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

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA