#source("http://bioconductor.org/biocLite.R")
#biocLite("limma")

logFoldChange=1
adjustP=0.05

library(limma)
setwd("C:\\Users\\lexb4\\Desktop\\geoLncRNA\\05.diff")
rt=read.table("biotype.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=avereps(rt)
rt=normalizeBetweenArrays(as.matrix(rt))
rt=log2(rt+1)

#differential
Type=c(rep("treat",13),rep("con",10))
design <- model.matrix(~0+factor(Type))
colnames(design) <- c("con","treat")
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)
type=sapply(strsplit(rownames(allDiff),"\\|"),"[",2)
protein=allDiff[type=="protein_coding",]
lncRNA=allDiff[type=="lncRNA",]
rownames(allDiff)=gsub("(.*?)\\|.*","\\1",rownames(allDiff))
rownames(protein)=gsub("(.*?)\\|.*","\\1",rownames(protein))
rownames(lncRNA)=gsub("(.*?)\\|.*","\\1",rownames(lncRNA))
rownames(rt)=gsub("(.*?)\\|.*","\\1",rownames(rt))

#write protein and lncRNA table
diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffSig,file="diff.xls",sep="\t",quote=F)
#write expression level of diff gene
hmExp=rt[rownames(diffSig),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="heatmap.txt",sep="\t",quote=F,col.names=F)

#write lncRNA table
lncRNASig <- lncRNA[with(lncRNA, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
write.table(lncRNASig,file="lncRNA_diff.xls",sep="\t",quote=F)
#write expression level of diff lncRNA gene
hmExp=rt[rownames(lncRNASig),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="lncRNA_heatmap.txt",sep="\t",quote=F,col.names=F)

#write protein table
proteinSig <- protein[with(protein, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
write.table(proteinSig,file="protein_diff.xls",sep="\t",quote=F)
#write expression level of diff gene
hmExp=rt[rownames(proteinSig),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="protein_heatmap.txt",sep="\t",quote=F,col.names=F)

#volcano
tiff(file="vol.tiff",
       width = 12,            
       height =12,            
       units ="cm",
       compression="lzw",
       bg="white",
       res=600)
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()