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

本网对Markdown排版支持较差，对格式不满意的用户请跳转至 或“宏基因组”公众号阅读；

写在前面

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

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

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

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

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

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

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

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

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

1箱线图：Alpha多样性
2散点图：Beta多样性，PCoA, CCA
3热图：差异菌、OTU及功能

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

曼哈顿图：差异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$ididx = rownames(taxonomy) %in% x$otutax = taxonomy[idx, ] # subset taxonomy from used OTU# 手动筛选显著的组x = x[rownames(tax), ] # reorder according to taxx$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 orderx$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")pggsave(file=paste("man_otu_", sampleA, "vs", sampleB, ".pdf", sep=""), p, width = 10, height = 3, useDingbats=F)ggsave(file=paste("man_otu_", sampleA, "vs", sampleB, ".png", sep=""), p, width = 10, height = 3)

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

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

想了解更多宏基因组、16S文献阅读和分析相关文章，快关注“宏基因组”公众号，干货第一时间推送。

系统学习生物信息，快关注“生信宝典”，那里有几千志同道合的小伙伴一起学习。