#!/bin/Rscript
# ***************************************************************
# Name:      taxo_NMDS_plot.R
# Purpose:   Generates the NMDS plot for microbial data
# Version:   0.2
# Authors:   Umer Zeeshan Ijaz (Umer.Ijaz@glasgow.ac.uk)
#                 http://userweb.eng.gla.ac.uk/umer.ijaz
#            Christopher Quince (Christopher.Quince@glasgow.ac.uk)
#                 http://userweb.eng.gla.ac.uk/christopher.quince
# Created:   2012-09-10
# License:   Copyright (c) 2012 Computational Microbial Genomics Group, University of Glasgow, UK
# 
#            This program is free software: you can redistribute it and/or modify
#            it under the terms of the GNU General Public License as published by
#            the Free Software Foundation, either version 3 of the License, or
#            (at your option) any later version.
#
#            This program is distributed in the hope that it will be useful,
#            but WITHOUT ANY WARRANTY; without even the implied warranty of
#            MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#            GNU General Public License for more details.
#
#            You should have received a copy of the GNU General Public License
#            along with this program.  If not, see <http://www.gnu.org/licenses/>.
# **************************************************************/
suppressPackageStartupMessages(library("optparse"))
suppressPackageStartupMessages(library("calibrate"))
suppressPackageStartupMessages(library("vegan"))
suppressPackageStartupMessages(library("RColorBrewer"))
suppressPackageStartupMessages(library("gplots"))

#specify desired options in a list
option_list <- list(
make_option("--ifile", action="store",default=NULL, help="CSV file"),
make_option("--opath", action="store",default=NULL, help="Output path"),
make_option("--fsize", action="store", default="1.2", help="Font size [default %default]"),
make_option("--width", type="integer",default=800, help="Width of jpeg files [default %default]"),
make_option("--height", type="integer",default=800, help="Height of jpeg files [default %default]"),
make_option("--mreads", type="integer",default=200, help="Minimum reads [default %default]"),
make_option("--ccol", type="integer",default=1, help="Minimum reads [default %default]"),
make_option("--method", type="integer",default=6, help="Dissimilarity matrix: 1=euclidean, 2=manhattan, 3=gower, 4=altGower, 5=canberra, 6=bray, 7=kulczynski, 8=morisita,9=horn, 10=binomial, 11=cao [default %default]"),
make_option("--rmode", action="store_true",default=FALSE, help="Mode: TRUE=R mode, FALSE=Q mode [default %default]"),
make_option("--dmode", action="store_true",default=FALSE, help="Site names display modes [default %default]")
)

#get command line options
opt<-parse_args(OptionParser(usage="%prog [options] file", option_list=option_list))

if(is.null(opt$ifile))
quit()

#generate colors function
generateColors <- function(siteNames,column) {
    
    #generate the indices
    colInd<-""
    
    if(column==0){
        colInd<-rep("blue",length(siteNames))
    }
    else{
        #first split the names with respect to underscore and construct a parsedString matrix
        for(i in seq(1:length(siteNames)))
        {
            tmp<-t(data.matrix(unlist(strsplit(siteNames[i],"_"))))
            if(column>length(tmp)){
                column<-length(tmp)
            }
            if(i==1){
                parsedString<-tmp
            }
            else{
                parsedString<-rbind(parsedString,tmp)
            }
        }
        #generate template colors
        crp<- colorRampPalette(c("blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red"),space = "Lab")(length(unique(parsedString[,column])))        
        for(i in seq(1:length(siteNames)))
        {
            colInd[i]<-crp[match(parsedString[i,column],unique(parsedString[,column]))]
        }
    }
    return(colInd)
}



#Importing data now
AS_C <-read.csv(opt$ifile,header=TRUE,row.names=1)

#transpose the data
AS <- t(AS_C)

#filter the data for reads smaller than a particular number
AS<-AS[rowSums(AS)>opt$mreads,]

#Transpose the data if R Mode
if(opt$rmode){
    AS=t(AS)
}

method<-switch(opt$method,"euclidean","manhattan","gower","altGower","canberra","bray","kulczynski","morisita","horn","binomial","cao")

MDS_res=metaMDS(AS, distance = method, k = 2, trymax = 50)
MDS_res.nmds.plot <- plot(MDS_res, dis = "sites", type = "n")

# save in a file
jpeg(filename = paste(opt$opath,substr(basename(opt$ifile),1,nchar(basename(opt$ifile))-4),"_NMDS.jpg",sep=""),width = as.numeric(opt$width), height = as.numeric(opt$height), quality=100)


par(mai=c(2,2,1,0.2)
    ,bg="white")

plot(MDS_res$points, t="n",
     cex.axis=as.numeric(opt$fsize)*0.6,
     cex.lab=as.numeric(opt$fsize),
     cex=as.numeric(opt$fsize)*0.7,
     cex.main=as.numeric(opt$fsize),
     xlim=c(min(MDS_res.nmds.plot$sites[,1])-1,max(MDS_res.nmds.plot$sites[,1]))+0.5,ylim=c(min(MDS_res.nmds.plot$sites[,2]),max(MDS_res.nmds.plot$sites[,2])),
     xlab="NMDS1", ylab="NMDS2")
points(MDS_res.nmds.plot$sites[,1],MDS_res.nmds.plot$sites[,2],xlim=c(min(MDS_res.nmds.plot$sites[,1])-1,max(MDS_res.nmds.plot$sites[,1]))+0.5,ylim=c(min(MDS_res.nmds.plot$sites[,2]),max(MDS_res.nmds.plot$sites[,2])),col=generateColors(rownames(AS),opt$ccol), pch = 19,cex=as.numeric(opt$fsize)*0.7)
if(opt$dmode){
    textxy(MDS_res.nmds.plot$sites[,1],MDS_res.nmds.plot$sites[,2],rownames(MDS_res.nmds.plot$sites),cx=as.numeric(opt$fsize)*0.6)
}
dev.off()