library(limma)#加载差异分析包limma
#将分组文件加载到环境中，分组信息第一列为样本名，第二列为分组信息如“high”“low”
targets<-read.csv("group.csv")
#将表达矩阵加载到环境中，行为基因，列为样本，这里应该注意去除重复项。
eset<-read.csv("expreset-basal1.csv",row.names = "symbol")
targets$Target=gsub("_",".",targets$Target)##若数据中存在特殊符号，将"_"替换成“.”，也可以不替换
##该数据集中实际存在不符合R的命名原则，所以在没个分类前加一个“F”，具体自己定
targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))
colnames(targets)=c("FileName","Target")#更改列名，为了和limma包中的一致
lev<-unique(targets$Target)##使用unique（）函数进行去重
f <- factor(targets$Target, levels=lev)
design <- model.matrix(~0+f)
colnames(design) <- lev
cont.wt <- makeContrasts("high-low",
+ levels=design)
fit <- lmFit(eset, design)#前面矩阵的row.name=“symbol”
fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
tT=topTable(fit2, adjust="BH",sort.by="logFC",n=Inf)
tT = subset(tT, select=c("adj.P.Val","P.Value","logFC"))
colnames(tT)=c("FDR","P.Value","logFC")
write.csv(tT,"DEGbasal.csv")
#最后的tT就是得到的差异基因，其中包含基因，P.Value和logFC
#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文件如下所示
