16s分析之差异展示（热图）
前两天我向大家推了16s做差异分析的两个包（没有看的请点击下面链接）：
差异做出来了如何展示，也是一个值得思考的问题，所以今天我们就尝试一下热图，看看效果：
#首先安装这个包
#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)
最终成图欣赏：