logFoldChange=2
adjustP=0.05
library(limma)
setwd("C:\\Users\\Administrator\\Desktop\\07.diff")
rt=read.table("normalize.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
#differential
modType=c(rep("normal",9),rep("OA",10),rep("normal",10),rep("OA",10),rep("normal",10))
design <- model.matrix(~0+factor(modType))
colnames(design) <- c("OA","normal")
fit <- lmFit(rt,design)
cont.matrix<-makeContrasts(OA-normal,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
allDiff=topTable(fit2,adjust='fdr',number=200000)
write.table(allDiff,file="limmaTab.xls",sep="\t",quote=F,row.names=F)
#write table
diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffSig,file="diff.xls",sep="\t",quote=F,row.names=F)
diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F)
diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]
write.table(diffDown,file="down.xls",sep="\t",quote=F,row.names=F)
#write expression level of diff gene
hmExp=rt[as.vector(diffSig[,1]),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="diffExp.txt",sep="\t",quote=F,col.names=F)
#volcano
pdf(file="vol.pdf")
xMax=max(-log10(allDiff$adj.P.Val))
yMax=max(abs(allDiff$logFC))
plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC",
     main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.8)
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange)
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8)
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange))
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="green",cex=0.8)
abline(h=0,lty=2,lwd=3)
dev.off()
setwd("C:\\Users\\Administrator\\Desktop\\OAS\\08.pheatmap")  rt=read.table("diffExp.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
#rt=log2(rt+1)
#rt[rt>15]=15
library(pheatmap)
Geo=c(rep("GSE12021",19),rep("GSE55235",20),rep("GSE55457",20))
Type=c(rep("N",9),rep("OA",10),rep("N",10),rep("OA",10),rep("N",10),rep("OA",10))   names(Geo)=colnames(rt)
ann=cbind(Geo,Type)
ann=as.data.frame(ann)
tiff(file="heatmap.tiff",
     width = 20,            
     height =25,           
     units ="cm",
     compression="lzw",
     bg="white",
     res=500)
pheatmap(rt, annotation=ann, 
         color = colorRampPalette(c("green", "black", "red"))(50),
         cluster_cols =F,
         #         scale="row",
         fontsize_row=3,
         fontsize_col=5)
dev.off()