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

#biocLite("edgeR")

#install.packages("gplots")

foldChange=2

padj=0.01

setwd("C:\\Users\\huagb\\Desktop\\TCGA\\03miRNA")                    #???ù???¼

library("edgeR")

rt=read.table("miRNAmatrix.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))

data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)

data=avereps(data)

data=data[rowMeans(data)>1,]

group=c(rep("normal",3),rep("tumor",309))                         #????????????????

design <- model.matrix(~group)

y <- DGEList(counts=data,group=group)

y <- calcNormFactors(y)

y <- estimateCommonDisp(y)

y <- estimateTagwiseDisp(y)

et <- exactTest(y,pair = c("normal","tumor"))

topTags(et)

ordered_tags <- topTags(et, n=100000)

allDiff=ordered_tags$table

allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]

diff=allDiff

newData=y$pseudo.counts

write.table(diff,file="edgerOut.xls",sep="\t",quote=F)

diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]

write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)

diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]

write.table(diffUp, file="up.xls",sep="\t",quote=F)

diffDown = diff[(diff$FDR < padj & (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)   #???????л???У??????????normalizeExp.txt??

diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])

write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F)         #????????????У??????????diffmRNAExp.txt??

heatmapData <- newData[rownames(diffSig),]

#volcano

#pdf(file="vol.pdf")

tiff(file="vol.tiff",width =12,height =12,units ="cm",compression="lzw",bg="white",res=400)

xMax=max(-log10(allDiff$FDR))+1

yMax=12

plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",

     main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)

diffSub=allDiff[allDiff$FDRfoldChange,]

points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)

diffSub=allDiff[allDiff$FDR<(-foldChange),]

points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)

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)

tiff(file="heatmap.tiff",width =60,height =90,units ="cm",compression="lzw",bg="white",res=400)

par(oma=c(10,3,3,7))

heatmap.2(hmMat,col='greenred',trace="none")

dev.off()

#install.packages("VennDiagram")

setwd("C:\\Users\\huagb\\Desktop\\TCGA\\05miRNA_Target\\03mRNA_Venny")

files=dir()

files=grep("txt",files,value=T)

targetList=list()

for(i in 1:length(files)){

    inputFile=files[i]

    rt=read.table(inputFile,header=F)

    header=unlist(strsplit(inputFile,"\\.|\\-"))

    targetList[[header[1]]]=as.vector(rt[,1])

    uniqLength=length(unique(as.vector(rt[,1])))

    print(paste(header[1],uniqLength,sep=" "))

}

library(VennDiagram)

venn.diagram(targetList,filename="venny.tiff",imagetype = "tiff",

             fill=rainbow(length(targetList)),cat.cex=0.6)

intersectGenes=Reduce(intersect,targetList)

write.table(file="target.xls",intersectGenes,sep="\t",quote=F,col.names=F,row.names=F)