扩增子统计绘图7三元图

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

写在前面

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

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

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

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

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

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

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

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

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

1箱线图：Alpha多样性
2散点图：Beta多样性，PCoA, CCA
3热图：差异菌、OTU及功能
4曼哈顿图：差异OTU或Taxonomy
5火山图：差异OTU数量及变化规律
6韦恩图：比较组间共有和特有OTU或分类单元
本节需要在”3热图：差异菌、OTU及功能”和”6韦恩图”基础上继续运行

7三元图

三元图有两种用法，常用的本质上是维恩图的一种变形，但维恩图只是数字比较单调，三元图类型上是散点图，可以用点大小和颜色代表丰度、显著性等信息，来进一步丰富图片信息。而且中国的文化中有事不过三，一而再、再而三等文化；三角形成最稳定也给人稳重、信认之感，三元图的美观和实用自然必不可少。

三元图主要分两种：展示两组共有和特有显著富集OTU，展示三种特异富集OTU。详见7三元图：美的不要不要的，再多用也不过分。本文主要以绘制比较常用的两组共有和特有显著富集OTU的三角图。另一种只要分析思路清楚，大家在此基础上很容易修改出来，只是代码量是需要加倍的(6次两两比较+三次取交集)。

加载三元图的配色方案和自定义函数

# 定义常用颜色 Defined color with transparent alpha = .7 c_yellow =          rgb(255 / 255, 255 / 255,   0 / 255, alpha) c_blue =            rgb(  0 / 255, 000 / 255, 255 / 255, alpha) c_orange =          rgb(255 / 255,  69 / 255,   0 / 255, alpha) c_green =           rgb(  50/ 255, 220 / 255,  50 / 255, alpha) c_dark_green =      rgb( 50 / 255, 200 / 255, 100 / 255, alpha) c_very_dark_green = rgb( 50 / 255, 150 / 255, 100 / 255, alpha) c_sea_green =       rgb( 46 / 255, 129 / 255,  90 / 255, alpha) c_black =           rgb(  0 / 255,   0 / 255,   0 / 255, alpha) c_grey =            rgb(180 / 255, 180 / 255,  180 / 255, alpha) c_dark_brown =      rgb(101 / 255,  67 / 255,  33 / 255, alpha) c_red =             rgb(200 / 255,   0 / 255,   0 / 255, alpha) c_dark_red =        rgb(255 / 255, 130 / 255,   0 / 255, alpha) # 三元图函数，无须理解直接调用即可 Function of ternary plot tern_e=function (x, scale = 1, dimnames = NULL, dimnames_position = c("corner",                                                                       "edge", "none"), dimnames_color = "black", id = NULL, id_color = "black",                  coordinates = FALSE, grid = TRUE, grid_color = "gray", labels = c("inside",                                                                                    "outside", "none"), labels_color = "darkgray", border = "black",                  bg = "white", pch = 19, cex = 1, prop_size = FALSE, col = "red",                  main = "ternary plot", newpage = TRUE, pop = TRUE, ...) {  labels = match.arg(labels)  if (grid == TRUE)    grid = "dotted"  if (coordinates)    id = paste("(", round(x[, 1] * scale, 1), ",", round(x[,                                                            2] * scale, 1), ",", round(x[, 3] * scale, 1), ")",                sep = "")  dimnames_position = match.arg(dimnames_position)  if (is.null(dimnames) && dimnames_position != "none")    dimnames = colnames(x)  if (is.logical(prop_size) && prop_size)    prop_size = 3  if (ncol(x) != 3)    stop("Need a matrix with 3 columns")  if (any(x < 0))    stop("X must be non-negative")  s = rowSums(x)  if (any(s <= 0))    stop("each row of X must have a positive sum")  x = x/s  top = sqrt(3)/2  if (newpage)    grid.newpage()  xlim = c(-0.03, 1.03)  ylim = c(-1, top)  pushViewport(viewport(width = unit(1, "snpc")))  if (!is.null(main))    grid.text(main, y = 0.9, gp = gpar(fontsize = 18, fontstyle = 1))  pushViewport(viewport(width = 0.8, height = 0.8, xscale = xlim,                        yscale = ylim, name = "plot"))  eps = 0.01  grid.polygon(c(0, 0.5, 1), c(0, top, 0), gp = gpar(fill = bg,                                                     col = border), ...)  if (dimnames_position == "corner") {    grid.text(x = c(0, 1, 0.5), y = c(-0.02, -0.02, top +                                        0.02), label = dimnames, gp = gpar(fontsize = 12))  }  if (dimnames_position == "edge") {    shift = eps * if (labels == "outside")      8    else 0    grid.text(x = 0.25 - 2 * eps - shift, y = 0.5 * top +                shift, label = dimnames[2], rot = 60, gp = gpar(col = dimnames_color))    grid.text(x = 0.75 + 3 * eps + shift, y = 0.5 * top +                shift, label = dimnames[1], rot = -60, gp = gpar(col = dimnames_color))    grid.text(x = 0.5, y = -0.02 - shift, label = dimnames[3],              gp = gpar(col = dimnames_color))  }  if (is.character(grid))    for (i in 1:4 * 0.2) {      grid.lines(c(1 - i, (1 - i)/2), c(0, 1 - i) * top,                 gp = gpar(lty = grid, col = grid_color))      grid.lines(c(1 - i, 1 - i + i/2), c(0, i) * top,                 gp = gpar(lty = grid, col = grid_color))      grid.lines(c(i/2, 1 - i + i/2), c(i, i) * top, gp = gpar(lty = grid,                                                               col = grid_color))      if (labels == "inside") {        grid.text(x = (1 - i) * 3/4 - eps, y = (1 - i)/2 *                    top, label = i * scale, gp = gpar(col = labels_color),                  rot = 120)        grid.text(x = 1 - i + i/4 + eps, y = i/2 * top -                    eps, label = (1 - i) * scale, gp = gpar(col = labels_color),                  rot = -120)        grid.text(x = 0.5, y = i * top + eps, label = i *                    scale, gp = gpar(col = labels_color))      }      if (labels == "outside") {        grid.text(x = (1 - i)/2 - 6 * eps, y = (1 - i) *                    top, label = (1 - i) * scale, gp = gpar(col = labels_color))        grid.text(x = 1 - (1 - i)/2 + 3 * eps, y = (1 -                                                      i) * top + 5 * eps, label = i * scale, rot = -120,                  gp = gpar(col = labels_color))        grid.text(x = i + eps, y = -0.05, label = (1 -                                                     i) * scale, vjust = 1, rot = 120, gp = gpar(col = labels_color))      }    }  xp = x[, 2] + x[, 3]/2  yp = x[, 3] * top  size = unit(if (prop_size)    #emiel inserted this code. x are proportions per row.  x*s is original data matrix. s = rowsums of original data matrix (x*s)    prop_size * rowSums(x*x*s) / max(  rowSums(x*x*s) )    #prop_size * rowSums(    (x*s) * ((x*s)/s)) / max(  rowSums(    (x*s) * ((x*s)/s)) )    else cex, "lines")  grid.points(xp, yp, pch = pch, gp = gpar(col = col), default.units = "snpc",              size = size, ...)  if (!is.null(id))    grid.text(x = xp, y = unit(yp - 0.015, "snpc") - 0.5 *                size, label = as.character(id), gp = gpar(col = id_color,                                                          cex = cex))  if (pop)    popViewport(2)  else upViewport(2) }

绘制三组比较三元图，WT对照为顶点

# merge group to mean ## 按样品名合并实验组与转置的OTU mat_t2 = merge(sub_design[c("genotype")], t(norm), by="row.names")[,-1] ## 按实验设计求组平均值 mat_mean = aggregate(mat_t2[,-1], by=mat_t2[1], FUN=mean) # mean # 重新转载并去除组名 per3=t(mat_mean[,-1]) colnames(per3) = mat_mean$genotype per3=as.data.frame(per3[rowSums(per3)>0,]) # remove all 0 OTU #per3=per3[,tern] # reorder per3 as input color=c(c_green,c_orange,c_red,c_grey) # 两底角相对于顶点显著富集的OTU，分共有和特有，类似维恩图 per3$color=color[4] # set all default # 设置点默认颜色为灰 AvC = KO_enriched BvC = OE_enriched C = intersect(row.names(AvC), row.names(BvC)) A = setdiff(AvC, C) B = setdiff(BvC, C) if (length(A)>0){per3[A,]$color=color[1]} if (length(B)>0){per3[B,]$color=color[2]} if (length(C)>0){per3[C,]$color=color[3]} ## output pdf and png in 8x8 inches per3lg=log2(per3[,1:3]*100+1) # 对数变换，剩数字可以调整OTU千分比的差距，点大小更均匀 pdf(file=paste("ter_",tern[1],tern[2],tern[3],"venn.pdf", sep=""), height = 8, width = 8) tern_e(per3lg[,1:3], prop=T, col=per3$color, grid_color="black", labels_color="transparent", pch=19, main="Tenary Plot") dev.off()

此图展示KO和OE突变体中特异或共有富集的OTU，本数据因为是测试数据，过统计显著的每组只有一个显著富集的OTU，没有共有富集的OTU。
详细的图片讲解，可参考7三元图：美的不要不要的，再多用也不过分

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

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