扩增子统计绘图4曼哈顿图：差异OTU和Taxonomy

点击上方蓝色「宏基因组」关注我们！专业干货每日推送！

写在前面

优秀的作品都有三部分曲，如骇客帝国、教父、指环王等。

扩增子系列课程也分为三部曲：

第一部《扩增子图表解读》：加速大家对同行文章的解读能力。

第二部《扩增子分析解读》：学习数据分析的基本思路和流程。

第三部《扩增子统计绘图》：即是对结果进行可视和统计检验，达到出版级的图表结果。

《扩增子统计绘图》系列文章介绍

《扩增子统计绘图》是之前发布的《扩增子图表解读》和《扩增子分析解读》的进阶篇，是在大家可以看懂文献图表，并能开展标准扩增子分析的基础上，进行结果的统计与可视化。其章节设计与《扩增子图表解读》对应，为八节课八种常用图形(箱线图、散点图、热图、曼哈顿图、火山图、维恩图、三元图和网络图)，基本满足文章常用的图片种类需求。

也适合对公司标准化分析返回结果的进一步统计、可视化及美化，达到出版级别，冲击高分文章。

本部分练习所需文件位于百度网盘，链接：http://pan.baidu.com/s/1hs1PXcw 密码：y33d。

本节需要上一节热图基础上继续运行

曼哈顿图：差异OTU和Taxonomy

## Manhattan图展示差异OTU和Taxonomy # 读取taxonomy，并添加各列名称 taxonomy = read.delim("rep_seqs_tax.txt", row.names= 1,header=F, sep="\t") colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus","species","evalue") # 标记差异OTU类型 x$level = as.factor(ifelse(x$sig==1, "enriched",ifelse(x$sig==-1, "depleted","nosig"))) x$otu = rownames(x) # 转换Pvalue为负对数 x$neglogp = -log(x$PValue) # Taxonomy排序，并筛选OTU表中存在的 library(dplyr) taxonomy$id=rownames(taxonomy) taxonomy = arrange(taxonomy, phylum, class, order, family, genus, species) rownames(taxonomy) = taxonomy$id idx = rownames(taxonomy) %in% x$otu tax = taxonomy[idx, ] # subset taxonomy from used OTU # 手动筛选显著的组 x = x[rownames(tax), ] # reorder according to tax x$tax = gsub("p__","",tax$phylum,perl=TRUE) top_phylum=c("Bacteroidetes","Firmicutes","Planctomycetes","Proteobacteria","Verrucomicrobia") x[!(x$tax %in% top_phylum),]$tax = "Low Abundance" # no level can get value # 设置各类的level对应顺序 x$otu = factor(x$otu, levels=x$otu)   # set x order x$level = factor(x$level, levels=c("enriched","depleted","nosig")) levels(x$tax)=c(top_phylum,"Low Abundance") # 调整Y轴范围更美观 x[x$neglogp>15,]$neglogp  = 15 # Manhattan plot ## 添加显著阈值线 FDR = min(x$neglogp[x$level=="depleted"]) library(ggplot2) p = ggplot(x, aes(x=otu, y=neglogp, color=tax, size=logCPM, shape=level)) +  geom_point(alpha=.7) +  geom_hline(yintercept=FDR, linetype=2, color="lightgrey") +  scale_shape_manual(values=c(17, 25, 20))+  scale_size(breaks=c(5, 10, 15)) +  labs(x="OTU", y="-loge(P)") +  theme(axis.ticks.x=element_blank(),axis.text.x=element_blank(),legend.position="top") p ggsave(file=paste("man_otu.pdf", sep=""), p, width = 10, height = 3, useDingbats=F) ggsave(file=paste("man_otu.png", sep=""), p, width = 10, height = 3)

图1. 曼哈顿图展示差异OTU所在的门。看到差异OTU以下调为主(空心下三角形)，其中以拟杆菌门(Bacteroidetes)密集下调。

详细的图片讲解，可参考4曼哈顿图：差异OTU或Taxonomy  

想了解更多宏基因组、16S分析相关文章，

快关注“宏基因组”公众号，干货第一时间推送。

系统学习生物信息，快关注“生信宝典”，

那里有几千志同道合的小伙伴一起学习。