#source("http://bioconductor.org/biocLite.R") #biocLite("limma") logFoldChange=0.8 adjustP=0.05 library(limma) setwd("D:\\Bioinformatics\\WGCNA\\DEG") 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("AAA",7),rep("Normal",8),rep("AAA",14),rep("Normal",8),rep("AAA",49),rep("Normal",10)) design <- model.matrix(~0+factor(modType)) colnames(design) <- c("treat","con") fit <- lmFit(rt,design) cont.matrix<-makeContrasts(treat-con,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) #write table diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ] write.table(diffSig,file="diff.xls",sep="\t",quote=F) diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ] write.table(diffUp,file="up.xls",sep="\t",quote=F) diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ] write.table(diffDown,file="down.xls",sep="\t",quote=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.VallogFoldChange) points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8) diffSub=subset(allDiff, adj.P.Val