16s分析之差异展示（热图）

前两天我向大家推了16s做差异分析的两个包（没有看的请点击下面链接）：

1.16s分析之差异分析（DESeq2）

2.16s分析之差异OTU 挑选（edgeR）

---

差异做出来了如何展示，也是一个值得思考的问题，所以今天我们就尝试一下热图，看看效果：

#首先安装这个包

#source("http://bioconductor.org/biocLite.R")

#biocLite("edgeR")

#install.packages("rlang")

#install.packages("vegan")

library("gplots")

library("RColorBrewer")

library(edgeR)

library("ggplot2")

library("vegan")

setwd("E:/Shared_Folder/HG_usearch_HG")

# 读入mapping文件

design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")

# 读取OTU表，这里我选择的是整个otu表格，但是一般没有必要全部做差异的啊，相对丰度高的做做就可以了

otu_table =read.delim("otu_table.txt", row.names= 1,sep="\t",header=T,check.names=F)

#本步骤采用otu。table基于抽平后的qiime文件，去掉第一行，和第二行首行的#符号，即可导入成功。

# 过滤数据并排序

idx =rownames(design) %in% colnames(otu_table)

sub_design =design[idx,]

count =otu_table[, rownames(sub_design)]

# 转换原始数据为百分比，

norm =t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100

# create DGE list

d = DGEList(counts=count,group=sub_design$SampleType)#count是一个数据框，$samples是第二个数据框

d =calcNormFactors(d)

# 生成实验设计矩阵

design.mat =model.matrix(~ 0 + d$samples$group)

colnames(design.mat)=levels(design$SampleType)

colnames(design.mat)

d2 =estimateGLMCommonDisp(d, design.mat)

d2 =estimateGLMTagwiseDisp(d2, design.mat)

fit = glmFit(d2,design.mat)

# 设置比较组

BvsA <-makeContrasts(contrasts = "GC1-GF1",levels=design.mat)#

# 组间比较,统计Fold change,Pvalue

lrt =glmLRT(fit,contrast=BvsA)

# FDR检验，控制假阳性率小于5%

de_lrt =decideTestsDGE(lrt, adjust.method="fdr", p.value=0.00005)

# 导出计算结果

x=lrt$table

head(x)

x$sig=de_lrt

enriched =row.names(subset(x,sig==1))

depleted =row.names(subset(x,sig==-1))

## 提取做差异的两个组的分组信息

pair_group =subset(sub_design, SampleType %in% c("GC1", "GF1"))

# 合并上升和下降的otu行名

DE=c(enriched,depleted)

#选择norm文件，就是行前面计算的相对丰度文件中的差异otu，列选择对应的组名

wt<-norm[DE,rownames(pair_group)]

str(wt)

# 绘图代码、预览、保存PDF

tiff(file="heatmap_otu.tif",res = 300, compression = "none", width=720,height=540,units="mm")

pdf("heatmap_otu.pdf")

# scale in row,dendrogram only in row, not cluster in column

p<-heatmap.2(wt,scale="row", Colv=F, Rowv=FALSE,dendrogram="none",

col=rev(colorRampPalette(brewer.pal(11, "RdYlGn"))(256)),

cexCol=1,keysize=1,density.info="none",main=NULL,trace="none",margins= c(6, 17))

这里还是对其中几个参数做一个解释吧：

dendrogram="none"：就是不聚类，就是没有下图的这个：

dendrogram="row"对列进行聚类

dendrogram="col"对行进行聚类，这里行我们表示的样品，可以来一个聚类：

当然参数还有很多，我先不调整了，我要换包做

dev.off()

#图形太丑，放弃使用

下面开始使用pheatmap

library(pheatmap)

#默认参数出图，发现图形还可以，添加的灰色边框也很好看

pheatmap(wt )

但是问题也是有的，数据颜色不是很分明

下面我们修改颜色和数据，方便展示：

colorRampPalette 函数支持自定义的创建一系列的颜色梯度；

#这里我设置藏青色，白色，和红色为渐变色，当然可以修改

color =colorRampPalette(c("navy", "white","firebrick3"))(30)

# cluster_rows横向无序聚类；fontsize=6字体大小我调节为6

pheatmap(wt,fontsize=6,cluster_rows= FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60) )

#这里我需要对数据进行变换，因为上面也看到了，我的数据都集中在蓝色的区域，我对让他们差异缩小一些

wt2<--log10(wt+0.00001)

#发现有很多极大值，特别红的那些

wt2<-sqrt(wt)

pheatmap(wt2,fontsize=6,cluster_rows= FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60) )

#还是不行

wt2<-sqrt(wt)

#我们把差异极大的点修改一下c小一些，另外图形需要窄一些，来适应期刊

wt2[wt2>0.9]<-0.9

pheatmap(wt2,fontsize=6,cellwidth= 10, cellheight = 10,cluster_rows = FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60) )

#加保存和题目

pheatmap(wt2,fontsize=6,cellwidth= 10, cellheight = 10,cluster_rows = FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60),

,main = "OTU heatmap",filename= "otu.pdf")

#下面我们来学习一下分分隔，横向分组我们使用聚类，纵向自己制定

#当然我们不能乱分组，这里我们建立两个分组文件

下面我减少了差异OTU的数目通过控制：

# FDR检验，控制假阳性率小于5%

de_lrt =decideTestsDGE(lrt, adjust.method="fdr", p.value=0.0000000000005)

# 下面命令和上面相同，再运行一次

x=lrt$table

head(x)

x$sig=de_lrt

enriched =row.names(subset(x,sig==1))

depleted =row.names(subset(x,sig==-1))

pair_group =subset(sub_design, SampleType %in% c("GC1", "GF1"))

DE=c(enriched,depleted)

wt<-norm[DE,rownames(pair_group)]

wt2<-sqrt(wt)

wt2[wt2>0.9]<-0.9

#我们支持分组,这里我随便造一个分组，作为纵向分组信息

annotation_row =data.frame(OTUClass = factor(rep(c("wt1", "wt2","wt3", "wt4"), c(7, 7,7,7))))

rownames(annotation_row)= rownames(wt2)

#我再造一个横向分组信息

annotation_col =data.frame(breaks = factor(rep(c("GC1", "GF1"),c(9,9))))

rownames(annotation_col)= c(paste("GC1", 1:9, sep = "") ,paste("GF5",1:9, sep = ""))

#运行代码

pheatmap(wt2,fontsize=6,cellwidth= 10, cellheight = 10,cluster_rows = FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60) ,

cutree_col = 2,gaps_row = c(5,25,30),

annotation_col =annotation_col,annotation_row = annotation_row)

出现错误：

我们的分隔信息和分组信息明不匹配，所以大家要特别注意

#修改分隔信息：

pheatmap(wt2,fontsize=6,cellwidth= 10, cellheight = 10,cluster_rows = FALSE,

color =colorRampPalette(c("navy", "white","firebrick3"))(60) ,

cutree_col = 2,gaps_row = c(7,14,21),

annotation_col =annotation_col,annotation_row = annotation_row)

最终成图欣赏：