foldChange=1
pval=0.05

setwd("C:\\Users\\ZXJ\\Desktop\\miRNA exosome\\peerj\\1差异分析") #Set up the working directory
library("edgeR")
rt=read.table("miRNA.7ren.txt",sep="\t",header=T,check.names=F)  #Load data set
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
dim(data)

# Filtering criteria: at least 75% of the samples had gene expression greater than 3
quant <- apply(data,1,quantile,0.25)
keep <- which((quant >= 3) == 1)
data <- data[keep,]
dim(data)
write.table(data,file="data.txt",sep="\t",quote=F,col.names=T)
group=c(rep("group1",2),rep("group2",5))                         #grouping
design <- model.matrix(~group)
y <- DGEList(counts=data,group=group)
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)
et <- exactTest(y,pair = c("group1","group2"))
topTags(et)
ordered_tags <- topTags(et, n=100000)

allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$PValue)==FALSE,]
diff=allDiff
newData=y$pseudo.counts

write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
diffSig = diff[(diff$PValue < pval & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
diffUp = diff[(diff$PValue < pval & (diff$logFC>foldChange)),]
write.table(diffUp, file="up.xls",sep="\t",quote=F)
diffDown = diff[(diff$PValue < pval & (diff$logFC<(-foldChange))),]
write.table(diffDown, file="down.xls",sep="\t",quote=F)

normalizeExp=rbind(id=colnames(newData),newData)
write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F)   
diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F)         

heatmapData <- newData[rownames(diffSig),]

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

#heatmap
hmExp=log10(heatmapData+0.001)
library('gplots')
hmMat=as.matrix(hmExp)
pdf(file="heatmap.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
heatmap.2(hmMat,col='bluered',trace="none")
dev.off()