#mRNA表达矩阵与GROUP文件样式，heatmap样式见文章最后
library(limma)
mRNA <- read.table("表达矩阵.txt",sep = "\t",header = T,comment.char = "!",encoding = "UTF-8")
#mRNA数据框行名为基因名，列命为样本名称
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顺序，实验组要在前面否则影响上下调结果
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结果
#将表达矩阵与差异分析结果合并
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)) #标注样本的分组信息
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()
#表达矩阵与GROUP文件如下所示
