扩增子统计绘图3热图：样品相关分析，差异OTU/ASV

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

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

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

本课程中所需的测序数据、实验设计；和课程分析生成的中间文件，均可以直去百度云下载。后台回复'扩增子'获取。

热图展示样品相关性

1. # 运行前，请在Rstudio中菜单栏选择“Session - Set work directory -- Choose directory”，弹窗选择之前分析目录中的result文件夹
2. # 读入实验设计
3. design = read.table('design.txt', header=T, row.names= 1, sep='\t')
4. # 读取OTU表
5. otu_table = read.delim('otu_table.txt', row.names= 1,  header=T, sep='\t')
6. # 过滤数据并排序
7. idx = rownames(design) %in% colnames(otu_table)
8. sub_design = design[idx,]
9. count = otu_table[, rownames(sub_design)]
10. # 转换原始数据为百分比
11. norm = t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100
12. # 计算所有样品间相关系数
13. sim=cor(norm,method='pearson')
14. # 使用热图可视化，并保存为8x8英寸的PDF
15. library('gplots')
16. library('RColorBrewer')
17. pdf(file=paste('heat_cor_samples.pdf', sep=''), height = 8, width = 8)
18. # 想预览，跳过上面Pdf行直接运行heatmap.2
19. heatmap.2(sim, Rowv=TRUE, Colv=TRUE, dendrogram='both', trace='none', margins=c(6,6), col=rev(colorRampPalette(brewer.pal(11, 'RdYlGn'))(256)),density.info='none')
20. dev.off()

图1. 热图展示所有样品基于相对丰度的Pearson相关系数矩阵。我们可以看到样品明显分成了三类，KO，OE，WT，表明该基因的过表达和基因敲除对菌群均有影响，其中过表达到WT差异明显。其中KO3与WT聚在了一起，表明其野生型相似，我能想到三种可能：过表达的基因被沉默而回复成与野生型相似；该份材料的种子是混入的WT；可能该材料的标WT串成了KO3。

edgeR统计组间差异OTU

# 使用edgeR统计组间差异OTU，以OE vs WT为例library(edgeR)# create DGE listd = DGEList(counts=count, group=sub_design$genotype)d = calcNormFactors(d)# 生成实验设计矩阵design.mat = model.matrix(~ 0 + d$samples$group)colnames(design.mat)=levels(genotypes)d2 = estimateGLMCommonDisp(d, design.mat)d2 = estimateGLMTagwiseDisp(d2, design.mat)fit = glmFit(d2, design.mat)# 设置比较组BvsA <- makeContrasts(contrasts = 'OE-WT', levels=design.mat)# 组间比较,统计Fold change, Pvaluelrt = glmLRT(fit,contrast=BvsA)# FDR检验，控制假阳性率小于5%de_lrt = decideTestsDGE(lrt, adjust.method='fdr', p.value=0.05)# 导出计算结果x=lrt$tablex$sig=de_lrtenriched = row.names(subset(x,sig==1))depleted = row.names(subset(x,sig==-1))

热图展示差异OTU

1. ## 热图展示差异OTU
2. pair_group = subset(sub_design, genotype %in% c('OE', 'WT'))
3. # Sig OTU in two genotype
4. DE=c(enriched,depleted)
5. sub_norm = as.matrix(norm[DE, rownames(pair_group)])
6. #colnames(sub_norm)=gsub('DM','KO',colnames(sub_norm),perl=TRUE) # rename samples ID
7. pdf(file=paste('heat_otu_OEvsWT_sig.pdf', sep=''), height = 8, width = 8)
8. # scale in row, dendrogram only in row, not cluster in column
9. heatmap.2(sub_norm, scale='row', Colv=FALSE, Rowv=FALSE,dendrogram='none', col=rev(colorRampPalette(brewer.pal(11, 'RdYlGn'))(256)), cexCol=1,keysize=1,density.info='none',main=NULL,trace='none')
10. dev.off()

图中可到OTU95在OE中高表达，而其它OTU均在OE中下调；表达该基因的表达，主要来拟制一些OTU。这些OTU的编号较大，代表其丰度较高。比如OTU_1，就是聚类结果中最高丰度的OTU。