microbiomeViz:绘制lefse结果中Cladogram「建议收藏」

microbiomeViz:绘制lefse结果中Cladogram「建议收藏」平日经常会分析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

大家好,又见面了,我是你们的朋友全栈君。

平日经常会分析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')))
V1 V2
k__Bacteria 100.00000
kBacteria|pFirmicutes 64.91753
kBacteria|pBacteroidetes 35.08247
kBacteria|pFirmicutes|c__Clostridia 64.91753
kBacteria|pBacteroidetes|c__Bacteroidia 35.08247
kBacteria|pFirmicutes|cClostridia|oClostridiales 64.91753

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

Image

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

Image

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

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

使用实战

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

Image

输入数据为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

Image

差异物种注释

# 读取需要颜色标注的差异物种列表,本质上是两列和颜色对应表
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

Image

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

Reference

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

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

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/144822.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • redis集群主从复制原理_主从关系紫音

    redis集群主从复制原理_主从关系紫音Redis主从复制主从复制简介主从复制的概念主从复制的作用主从复制工作流程阶段一:建立连接阶段主从连接(slave连接master)第一种方式第二种方式第三种方式授权访问阶段二:数据同步阶段工作流程数据同步阶段master说明数据同步阶段slave说明阶段三:命令传播阶段命令传播阶段的部分复制服务器的运行id复制缓冲区复制缓冲区内部工作原理复制缓冲区主从服务器复制偏移量(offset)数据同步+命令传播阶段工作流程心跳机制心跳阶段注意事项主从复制常见问题引发频繁的全量复制1引发频繁的全量复制2频繁的网络中

  • 特殊多位数乘法口算算法

    特殊多位数乘法口算算法本文转自:我爱口算网,但是本人有一定更正,因此转载请注明出处一、关于9的数学口算技巧(两位数乘法)关于9的口诀:1×9=9  2×9=18  3×9=27    4×9=365×9=45  6×9=54  7×9=63    8×9=729×9=81上面的口诀小朋友们已经会了吗?小学

  • 从mssql (sqlserver2000)中导出数据到mysql 中用load data

    从mssql (sqlserver2000)中导出数据到mysql 中用load data

  • python3.6写一个http接口服务,给别人调用1

    python3.6写一个http接口服务,给别人调用1一、python3.6写一个http接口服务,给别人调用1首先推荐tornado,Tornado是一个Pythonweb框架和异步网络库,最初在FriendFeed开发。通过使用无阻塞网络I/O,Tornado可以扩展到数万个开放连接,使其成为长轮询、WebSocket和其他需要与每个用户建立长时间连接的应用程序的理想选择。简易而且本地win10能够跑起来。二、Torna…

  • elk面试题_百家公司运维面试题汇总

    elk面试题_百家公司运维面试题汇总备注:这一我在去年国庆节期间,整理的整个19年,学员的面试遇到的问题,整理出来之后发给后期的学员,让他们做参考和学习,看看公司会面试哪些问题。前言小的时候,哭着哭着就笑了;长大后笑着笑着就哭了,这是一种人生经历,当你经历的越多,你越发现世界不像童话里那么美好。真正值得在乎的东西,不会越来越多,只会越来越少,所以珍惜你当下的每一寸时光。现在的每一份努力,都会变成倍增的回收,在公众面前表现出来。距…

  • leetcode 最长有效括号_字符指针赋值为什么不能加大括号

    leetcode 最长有效括号_字符指针赋值为什么不能加大括号给你一个只包含 ‘(’ 和 ‘)’ 的字符串,找出最长有效(格式正确且连续)括号子串的长度。示例 1:输入:s = “(()”输出:2解释:最长有效括号子串是 “()”示例 2:输入:s = “)()())”输出:4解释:最长有效括号子串是 “()()”示例 3:输入:s = “”输出:0题解括号匹配:(看作+1,)看作-1,所有满足条件的括号应该是前缀和>=0,并且总和==0class Solution {public: const int INF =

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号