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=" ")

#讀取OTU表，這裡我選擇的是整個otu表格，但是一般沒有必要全部做差異的啊，相對丰度高的做做就可以了

otu_table =read.delim("otu_table.txt", row.names= 1,sep=" ",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

#組間比較,統計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

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

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

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

這裡還是對其中幾個參數做一個解釋吧：

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

#發現有很多極大值，特別紅的那些

wt2

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

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

#還是不行

wt2

#我們把差異極大的點修改一下c小一些，另外圖形需要窄一些，來適應期刊

wt2[wt2>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

wt2

wt2[wt2>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)

最終成圖欣賞：

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！