require(edgeR)
rnaseqMatrix0= read.table("/data/images/proton/DKlab/mr/Polysomes/reads/stringtie/gene_count_matrix.csv",header=T,sep=",")
r=rownames(rnaseqMatrix0)
flt=grep("ENSMUSG",r)
rnaseqMatrix0 =rnaseqMatrix0[flt,];
#data = read.table("/data/images/proton/run125/www/genes.counts.matrix", header=T, row.names=1, com='')
#col_ordering = c(1,2,6,7)

rnaseqMatrix =rnaseqMatrix0;
#rnaseqMatrix = round(rnaseqMatrix)
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
#conditions = factor(c(rep("FAGsA", 2), rep("FAGsB", 2)))
#                       1    2     3       4       5     6     7       8       9     10    11    12       13    14       15    16    17    18    19      20 
#conditions = factor(c("KO2","KO2","KOIFN","KOIFN","Ko0","Ko0","KoIL4","KoIL4","Ko6","Ko6","wt0","wt0IL4","wt0","wt0IL4","wt2","wt2","wt6","wt6","wtIFN","wtIFN"))
#                 1                  2                 3                  4                 5                  6                         7                  8                  9                  10                11                 12                13                 14                15                 16                       17                        18                 19                  20 
#,KO_0h_replicate_I,KO_0h_replicate_II,KO_2h_replicate_I,KO_2h_replicate_II,KO_6h_replicate_I,KO_6h_replicate_II,KO_IFN-gamma_replicate_II,KO_IFN_replicate_I,KO_IL4_replicate_I,KO_IL4_replicate_II,WT_0h_replicate_I,WT_0h_replicate_II,WT_2h_replicate_I,WT_2h_replicate_II,WT_6h_replicate_I,WT_6h_replicate_II,WT_IFN-gamma_replicate_I,WT_IFN-gamma_replicate_II,WT_IL4_replicate_I,WT_IL4_replicate_II
#conditions = factor(c("KO","KO","KO","KO","KO","KO","KO","KO","KO","KO","wt","wt","wt","wt","wt","wt","wt","wt","wt","wt"))

#exp_study = DGEList(counts=rnaseqMatrix, group=conditions)
#exp_study = calcNormFactors(exp_study)
#exp_study = estimateCommonDisp(exp_study)
#exp_study = estimateTagwiseDisp(exp_study)
#et = exactTest(exp_study)
#tTags = topTags(et,n=NULL,adjust.method = "fdr")
#r = rownames(tTags)
#c=data[r,col_ordering]
#tab <- topTags(lrt)
#write.table(tTags, file="edgeRgenes.csv")

col_ordering = c(1,2,11,12)#0h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
conditions = factor(c("KO","KO","WT","WT"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
detags <- rownames(exp_study)[as.logical(de)]
png("edgeRgenes-0h-MA.png",width=800,height=800)
plotSmear(et, de.tags=detags)
dev.off();
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-0h.csv")
write.table(x, file="edgeRgenesExpr-0h.csv")


col_ordering = c(3,4,13,14)#2h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
detags <- rownames(exp_study)[as.logical(de)]
png("edgeRgenes-2h-MA.png",width=800,height=800)
plotSmear(et, de.tags=detags)
dev.off();
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-2h.csv")
write.table(x, file="edgeRgenesExpr-2h.csv")

col_ordering = c(5,6,15,16)#6h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);
exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
detags <- rownames(exp_study)[as.logical(de)]
png("edgeRgenes-6h-MA.png",width=800,height=800)
plotSmear(et, de.tags=detags)
dev.off();
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-6h.csv")
write.table(x, file="edgeRgenesExpr-6h.csv")

col_ordering = c(7,8,17,18)#IFN
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);
exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
detags <- rownames(exp_study)[as.logical(de)]
png("edgeRgenes-IFN-MA.png",width=800,height=800)
plotSmear(et, de.tags=detags)
dev.off();
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-IFN.csv")
write.table(x, file="edgeRgenesExpr-IFN.csv")

col_ordering = c(9,10,19,20)#IL4
rnaseqMatrix =rnaseqMatrix0[,col_ordering];
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
#rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);
exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(tTags$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
detags <- rownames(exp_study)[as.logical(de)]
png("edgeRgenes-IL4-MA.png",width=800,height=800)
plotSmear(et, de.tags=detags)
dev.off();
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRgenes-IL4.csv")
write.table(x, file="edgeRgenesExpr-IL4.csv")