library(edgeR)

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 = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
conditions = factor(c(rep("FAGsA", 2), rep("FAGsB", 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)
c=data[r,col_ordering]

# Do a volcano plot
pdf("genes.counts.matrix.FAGsA_vs_FAGsB.edgeR.DE_results.MA_n_Volcano2.pdf", w=11, h=8.5)
plot(tTags$table$logFC, -log(tTags$table$PValue,10), main="Volcano plot of raw p-values", 
  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (raw p-value)")
#plot(tTags$table$logFC, -log(tTags$table$FDR,10), main="Volcano plot of FDR-adjusted p-values", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (FDR p-value)")
#plot(tTags$table$logFC, -log(tTags$table$PValue,10), main="Volcano plot showing genes with a raw p-value < 1e-05", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (raw p-value)", ylim=c(0,5))
#plot(tTags$table$logFC, -log(tTags$table$FDR,10), main="Volcano plot showing genes with a FDR p-value < 1e-05", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (FDR p-value)", ylim=c(0,5))
dev.off()


write.table(cbind(tTags,c), file='genes.counts.matrix.FAGsA_vs_FAGsB.edgeR.DE_results', sep='	', quote=F, row.names=T)
source("/data/results/tools/align/trinity/trinityrnaseq_r20140413p1/Analysis/DifferentialExpression/R/rnaseq_plot_funcs.R")
pdf("genes.counts.matrix.FAGsA_vs_FAGsB.edgeR.DE_results.MDS.pdf")
plotMDS(exp_study)
dev.off()
pdf("genes.counts.matrix.FAGsA_vs_FAGsB.edgeR.DE_results.MA_n_Volcano.pdf")
result_table = tTags$table
plot_MA_and_Volcano(result_table$logCPM, result_table$logFC, result_table$FDR)
dev.off()