# Rscript cummeRbund_Plots.R </path/to/cuffdiff/results> <1/0>


# Read data location from command line
args <- commandArgs(trailingOnly = TRUE)

if (length(args)!=2)
{
print("Wrong Arguements!!!")
print("Rscript cummeRbund_Plots.R </path/to/cuffdiff/results> <1 (Replicates) or 0 (Non Replicates)>")
q()
}

if (args[2]!=0)
	if (args[2]!=1)
	{
	print("Wrong Arguements!!!")
	print("Rscript cummeRbund_Plots.R </path/to/cuffdiff/results> <1 (Replicates) or 0 (Non Replicates)>")
	q()
	}



# Set the working directory accordingly
setwd(args[1])
Replicates=args[2]
# Load the cummeRbund and the Cuffdiff data
library(cummeRbund)
cuff<-readCufflinks()


# Output file is Output.pdf
Filename=paste("Output",unclass(Sys.time()),".pdf",sep="")
pdf(Filename,onefile=TRUE)




#how many significantly up or down regulated genes are there at various thresholds?
sig0 <- getSig(cuff, alpha=0.001, level='genes')
sig1 <- getSig(cuff, alpha=0.01, level='genes')
sig2 <- getSig(cuff, alpha=0.05, level='genes')

sample.names<-samples(genes(cuff))

# Count vs dispersion plot by condition for all genes
#disp<-dispersionPlot(genes(cuff))
#print(disp)

if (Replicates==1)
{
# squared coefficient of variation at the gene level
genes.scv<-fpkmSCVPlot(genes(cuff))
print(genes.scv)

# squared coefficient of variation at the isoform level
isoforms.scv<-fpkmSCVPlot(isoforms(cuff))
print(isoforms.scv)
}
# distributions of FPKM scores across samples
dens<-csDensity(genes(cuff))
print(dens)

# distributions of FPKM scores across samples with replicates
densRep<-csDensity(genes(cuff),replicates=T)
print(densRep)

# Boxplots
b<-csBoxplot(genes(cuff))
print(b)

# Boxplots with replicates
brep<-csBoxplot(genes(cuff),replicates=T)
print(brep)

s<-csScatterMatrix(genes(cuff))
print(s)

s<-csScatterMatrix(genes(cuff),replicates=T)
print(s)

# For all combination between different samples
for (i in 1:length(sample.names))
for (j in 1:length(sample.names))
{
if (i==j) next
#  individual Pairwise comparisons 
s<-csScatter(genes(cuff),sample.names[i], sample.names[j])
print(s)

# intensity vs fold-change plots
m<-MAplot(genes(cuff),sample.names[i], sample.names[j])
plot.new()
plot(m, main="")
title(main = paste(sample.names[i], sample.names[j]))

# normalized count value
mCount<-MAplot(genes(cuff),sample.names[i], sample.names[j],useCount=T)
plot.new()
plot(mCount, main="")
title(main = paste("normalized count values",sample.names[i], sample.names[j]))

if (Replicates==1)
{
#     csVolcano(object, x, y, alpha=0.05, showSignificant=TRUE,features=FALSE, xlimits = c(-20, 20), ...)
v<-csVolcano(genes(cuff),sample.names[i], sample.names[j])
print(v)
}

}
# Overall Volcano plot
v<-csVolcanoMatrix(genes(cuff))
print(v)

# dendrogram
dend<-csDendro(genes(cuff))
print(dend)
# dendrogram with Replicates
dend.rep<-csDendro(genes(cuff),replicates=T)
print(dend.rep)

# Distance matrix
myDistHeat<-csDistHeat(genes(cuff))
print(myDistHeat)

# Distance matrix with Replicates
myRepDistHeat<-csDistHeat(genes(cuff),replicates=T)
print(myRepDistHeat)

# Dimensionality reduction (PCA)
genes.PCA<-PCAplot(genes(cuff),"PC1","PC2")
print(genes.PCA)

if (Replicates==1)
{
# Dimensionality reduction (MDS)
#genes.MDS<-MDSplot(genes(cuff))
#print(genes.MDS)
}

# Dimensionality reduction (PCA) with Replicates
genes.PCA.rep<-PCAplot(genes(cuff),"PC1","PC2",replicates=T)
print(genes.PCA.rep)


#  at the gene level
if (length(sig0)>0)
sigGenes0 <- getGenes(cuff, sig0)

if (length(sig1)>0)
sigGenes1 <- getGenes(cuff, sig1)

if (length(sig2)>0)
sigGenes2 <- getGenes(cuff, sig2)


if (Replicates==1)
{
# Dimensionality reduction (MDS) with Replicates
genes.MDS.rep<-MDSplot(genes(cuff),replicates=T)
print(genes.MDS.rep)

# Overview of significant features
if (length(sig0)>0)
if (length(sigGenes0)>2)
{
mySigMat<-sigMatrix(cuff,level='genes',alpha=0.001)
print(mySigMat)
}

# Overview of significant features
if (length(sig1)>0)
if (length(sigGenes1)>2)
{
mySigMat<-sigMatrix(cuff,level='genes',alpha=0.01)
print(mySigMat)
}

# Overview of significant features
if (length(sig2)>0)
if (length(sigGenes2)>2)
{
mySigMat<-sigMatrix(cuff,level='genes',alpha=0.05)
print(mySigMat)
}
}


if (length(sig0)>0)
print(expressionBarplot(sigGenes0))

if (length(sig1)>0)
print(expressionBarplot(sigGenes1))



if (length(sig2)>0)
print(expressionBarplot(sigGenes2))


if (length(sig0)>2)
print(csHeatmap(sigGenes0, cluster='both'))

if (length(sig1)>2)
print(csHeatmap(sigGenes1, cluster='both'))


if (length(sig2)>2)
print(csHeatmap(sigGenes2, cluster='both'))


if (length(sig0)>0)
{
for (i in 1:length(sample.names))
for (j in 1:length(sample.names))
{
if (i==j) next
print(csScatter(sigGenes0, sample.names[i], sample.names[j], smooth=T))
}
}

if (length(sig1)>0)
{
for (i in 1:length(sample.names))
for (j in 1:length(sample.names))
{
if (i==j) next
print(csScatter(sigGenes1, sample.names[i], sample.names[j], smooth=T))
}
}


if (length(sig2)>0)
{
for (i in 1:length(sample.names))
for (j in 1:length(sample.names))
{
if (i==j) next
print(csScatter(sigGenes2, sample.names[i], sample.names[j], smooth=T))
}
}



dev.off()
print(paste("Output File: ",getwd(),"/",Filename,sep=""))