|
前兩天我向大家推了16s做差異分析的兩個(gè)包（沒(méi)有看的請(qǐng)點(diǎn)擊下面鏈接）： 差異做出來(lái)了如何展示,，也是一個(gè)值得思考的問(wèn)題,，所以今天我們就嘗試一下熱圖，看看效果： #首先安裝這個(gè)包 #source("http:///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表,，這里我選擇的是整個(gè)otu表格，但是一般沒(méi)有必要全部做差異的啊，相對(duì)豐度高的做做就可以了 otu_table =read.delim("otu_table.txt", row.names= 1,sep="\t",header=T,check.names=F) #本步驟采用otu。table基于抽平后的qiime文件,，去掉第一行，和第二行首行的#符號(hào),，即可導(dǎo)入成功,。 # 過(guò)濾數(shù)據(jù)并排序 idx =rownames(design) %in% colnames(otu_table) sub_design =design[idx,] count =otu_table[, rownames(sub_design)] # 轉(zhuǎn)換原始數(shù)據(jù)為百分比， 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是一個(gè)數(shù)據(jù)框,，$samples是第二個(gè)數(shù)據(jù)框 d =calcNormFactors(d) # 生成實(shí)驗(yàn)設(shè)計(jì)矩陣 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) # 設(shè)置比較組 BvsA <-makeContrasts(contrasts = "GC1-GF1",levels=design.mat)# # 組間比較,統(tǒng)計(jì)Fold change,Pvalue lrt =glmLRT(fit,contrast=BvsA) # FDR檢驗(yàn)，控制假陽(yáng)性率小于5% de_lrt =decideTestsDGE(lrt, adjust.method="fdr", p.value=0.00005) # 導(dǎo)出計(jì)算結(jié)果 x=lrt$table head(x) x$sig=de_lrt enriched =row.names(subset(x,sig==1)) depleted =row.names(subset(x,sig==-1)) ## 提取做差異的兩個(gè)組的分組信息 pair_group =subset(sub_design, SampleType %in% c("GC1", "GF1")) # 合并上升和下降的otu行名 DE=c(enriched,depleted) #選擇norm文件,，就是行前面計(jì)算的相對(duì)豐度文件中的差異otu,，列選擇對(duì)應(yīng)的組名 wt<-norm[DE,rownames(pair_group)] str(wt) # 繪圖代碼、預(yù)覽,、保存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)) 這里還是對(duì)其中幾個(gè)參數(shù)做一個(gè)解釋吧：
dev.off() #圖形太丑,，放棄使用 下面開始使用pheatmap library(pheatmap) #默認(rèn)參數(shù)出圖,，發(fā)現(xiàn)圖形還可以，添加的灰色邊框也很好看 pheatmap(wt ) 但是問(wèn)題也是有的,，數(shù)據(jù)顏色不是很分明 下面我們修改顏色和數(shù)據(jù),，方便展示： colorRampPalette 函數(shù)支持自定義的創(chuàng)建一系列的顏色梯度； #這里我設(shè)置藏青色,，白色,，和紅色為漸變色，當(dāng)然可以修改 color =colorRampPalette(c("navy", "white","firebrick3"))(30) # cluster_rows橫向無(wú)序聚類,；fontsize=6字體大小我調(diào)節(jié)為6 pheatmap(wt,fontsize=6,cluster_rows= FALSE, color =colorRampPalette(c("navy", "white","firebrick3"))(60) ) #這里我需要對(duì)數(shù)據(jù)進(jìn)行變換,，因?yàn)樯厦嬉部吹搅耍业臄?shù)據(jù)都集中在藍(lán)色的區(qū)域,，我對(duì)讓他們差異縮小一些 wt2<--log10(wt+0.00001) #發(fā)現(xiàn)有很多極大值,，特別紅的那些 wt2<-sqrt(wt) pheatmap(wt2,fontsize=6,cluster_rows= FALSE, color =colorRampPalette(c("navy", "white","firebrick3"))(60) ) #還是不行 wt2<-sqrt(wt) #我們把差異極大的點(diǎn)修改一下c小一些，另外圖形需要窄一些,，來(lái)適應(yīng)期刊 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") #下面我們來(lái)學(xué)習(xí)一下分分隔,，橫向分組我們使用聚類,，縱向自己制定 #當(dāng)然我們不能亂分組，這里我們建立兩個(gè)分組文件 下面我減少了差異OTU的數(shù)目通過(guò)控制： # FDR檢驗(yàn),，控制假陽(yáng)性率小于5% de_lrt =decideTestsDGE(lrt, adjust.method="fdr", p.value=0.0000000000005) # 下面命令和上面相同,，再運(yùn)行一次 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 #我們支持分組,這里我隨便造一個(gè)分組，作為縱向分組信息 annotation_row =data.frame(OTUClass = factor(rep(c("wt1", "wt2","wt3", "wt4"), c(7, 7,7,7)))) rownames(annotation_row)= rownames(wt2) #我再造一個(gè)橫向分組信息 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 = "")) #運(yùn)行代碼 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) 出現(xiàn)錯(cuò)誤：
我們的分隔信息和分組信息明不匹配,，所以大家要特別注意 #修改分隔信息： 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) 最終成圖欣賞：
|
|
|
