#mRNA表達(dá)矩陣與GROUP文件樣式菲嘴，heatmap樣式見文章最后
library(limma)
?mRNA <- read.table("表達(dá)矩陣.txt",sep = "\t",header = T,comment.char = "!",encoding = "UTF-8")
#mRNA數(shù)據(jù)框行名為基因名想虎，列命為樣本名稱
?group <- read.table("GROUP.txt",header=T,sep = "\t",encoding = "UTF-8")
?group_CP <- group$treat?
?m_design<- model.matrix(~0+factor(group_CP))
?colnames(m_design) = levels(factor(group_CP))
?rownames(m_design)= group$ID
?contrast.matrix<-makeContrasts("P-C",levels=m_design) #注意P-C順序窝革，實(shí)驗(yàn)組要在前面否則影響上下調(diào)結(jié)果
?m_fit <- lmFit(mRNA,m_design)
?m_fit <- contrasts.fit(m_fit, contrast.matrix)
?m_fit <- eBayes(m_fit)
?m_genlist <- topTable(m_fit, coef = 1, n=Inf)? #limma結(jié)果
#將表達(dá)矩陣與差異分析結(jié)果合并
? ID_REF <- rownames(m_genlist)
? m_genlist <- data.frame(ID_REF,m_genlist)
? ID_REF <- rownames(mRNA)
? mRNA <- data.frame(ID_REF,mRNA)
? test <-merge(mRNA,m_genlist,by = "ID_REF")
? result <- subset(test,P.Value<0.05)
? row.names(result) <- result[,1]
#繪制熱圖
heatmap <- result[2:(nrow(group)+1)]
annotation <- data.frame(Factor = factor(group$treat)) #標(biāo)注樣本的分組信息
rownames(annotation) <- colnames(heatmap)
library(pheatmap)
filename <- paste("文件名",".pdf",sep="")
pdf(filename)
pheatmap(heatmap,
? ? ? ? annotation=annotation,
? ? ? ? annotation_legend = TRUE,
? ? ? ? main=filename ,
? ? ? ? scale = "row",
? ? ? ? show_rownames = F,
? ? ? ? color = colorRampPalette(c("green","black","red"))(100))
dev.off()
#表達(dá)矩陣與GROUP文件如下所示
