#DGE for microarray by limma
library("edgeR")
library('WGCNA')
library('gplots')
library('limma')
setwd("C:/Users/lenovo/Desktop/Diabetic children/DEG（logFC=0.5,P=0.01）/DEG")
foldChange=0.5 
padj=0.01
rawexprSet=read.csv("express-counts2.csv",header=TRUE,row.names=1,check.names = FALSE)  
dim(rawexprSet)
exprSet=log2(rawexprSet)
par(mfrow=c(1,2))
boxplot(data.frame(exprSet),col="blue") 
exprSet[1:5,1:5]
group <- read.csv("datTraits.csv",header=TRUE,row.names=1,check.names = FALSE)
group <- group[,1]
design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(exprSet)

fit <- lmFit(exprSet,design)
cont.matrix<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)  
##eBayes() with trend=TRUE
tempOutput = topTable(fit2,coef=1,n=Inf,adjust="BH") 

nrDEG = na.omit(tempOutput) 

allDiff <- nrDEG
diff=allDiff
write.csv(diff, "limmaOut.csv")
diffSig = diff[(diff$P.Value < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
#write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
write.csv(diffSig, "diffSig.csv")
diffUp = diff[(diff$P.Value < padj & (diff$logFC>foldChange)),]
#write.table(diffUp, file="up.xls",sep="\t",quote=F)
write.csv(diffUp, "diffUp.csv")
diffDown = diff[(diff$P.Value < padj & (diff$logFC<(-foldChange))),]
#write.table(diffDown, file="down.xls",sep="\t",quote=F)
write.csv(diffDown, "diffDown.csv")

#volcano
pdf(file="vol.pdf") 
xMax=max(-log10(allDiff$P.Value))+1
yMax=5
plot(-log10(allDiff$P.Value), allDiff$logFC, xlab="-log10(P.Value)",ylab="logFC",
     main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
diffSub=allDiff[allDiff$P.Value<padj & allDiff$logFC>foldChange,]
points(-log10(diffSub$P.Value), diffSub$logFC, pch=20, col="red",cex=0.4)
diffSub=allDiff[allDiff$P.Value<padj & allDiff$logFC<(-foldChange),]
points(-log10(diffSub$P.Value), diffSub$logFC, pch=20, col="green",cex=0.4)
abline(h=0,lty=2,lwd=3)
dev.off()

#heatmap
diffSig=read.csv("diffSig.csv",header=TRUE,row.names=1,check.names = FALSE)  
rownames1=rownames(exprSet)
exprSet2=cbind(rownames1,exprSet)
rownames2=rownames(diffSig)
diffSig2=cbind(rownames2,diffSig)
a=merge(exprSet2,diffSig2,by.x = "rownames1",by.y = "rownames2")
rownames(a)=a$Var.2
exprSet=as.matrix(a[,2:37])

heatmapData <- exprSet
orderdiffsig = diffSig[order(diffSig$logFC,decreasing=T),]
rownames(orderdiffsig)=orderdiffsig[,1]
orderdiffsig=orderdiffsig[,-1]
heatmapData = heatmapData[rownames(orderdiffsig),]
hmExp=log10(heatmapData+0.001)
mergedColors = as.character(labels2colors(group,colorSeq=c("red","yellow")))
hmMat=as.matrix(t(hmExp))
pdf(file="heatmap.pdf",width=45,height=30)
par(oma=c(10,3,3,17))
heatmap.2(hmMat,col='greenred',key=T,keysize = 0.7,
          trace="none", cexCol=0.6, cexRow=2.5,dendrogram="col",Rowv=F,
          RowSideColors = mergedColors,scale="col")
dev.off()