#!/bin/Rscript
# ***************************************************************
# Name:      taxo_diversity_plot.R
# Purpose:   This script generates the ecological diversity indices and 
#            rarefaction species richness. The following indices are supported:
#                 -Shanon Index
#                 -Simpson Index
#                 -Inverse Simpson’s Index
#                 -Fisher’s logarithmic series’ alpha parameters
#                 -Pielou’s evenness
#            Furthermore, it also generates a csv file *_div.csv for these indices 
# 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("vegan"))
suppressPackageStartupMessages(library("mgcv"))

#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("--rmode", action="store_true",default=FALSE, help="Mode: TRUE=R mode, FALSE=Q mode [default %default]")
)

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

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

#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,]

#use complete.cases function to filter out those rows that have NA entry in AS
AS<-AS[complete.cases(AS),]

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

H <- diversity(AS)
simp <- diversity(AS, "simpson")
invsimp <- diversity(AS, "inv")
r.2 <- rarefy(AS, 2)
alpha <- fisher.alpha(AS)

# save in a file
jpeg(filename = paste(opt$opath,substr(basename(opt$ifile),1,nchar(basename(opt$ifile))-4),"_DI_BP.jpg",sep=""),width = as.numeric(opt$width), height = as.numeric(opt$height), quality=100)
par( mai=c(2,2,2,0.2)
,bg="white")
pairs(cbind(H, simp, invsimp, r.2, alpha), labels=c("Shanon","Simpson","Inverse Simpson","Rarefied",expression(paste("Fisher's ",alpha))),pch=19,
    cex.labels=as.numeric(opt$fsize),
    font.labels=as.numeric(opt$fsize),
    cex.axis=as.numeric(opt$fsize)*0.6,  
    cex=as.numeric(opt$fsize)*0.7, col="black")
dev.off()

## Species richness (S) and Pielou's evenness (J):
S <- specnumber(AS) ## rowSums(BCI > 0) does the same...
J <- H/log(S)



# save in a file
jpeg(filename = paste(opt$opath,substr(basename(opt$ifile),1,nchar(basename(opt$ifile))-4),"_DI_INDICES.jpg",sep=""),width = as.numeric(opt$width), height = as.numeric(opt$height), quality=100)
par( mfrow=c(5,1)
     ,mai=c(as.numeric(opt$fsize),2,2,0.2)
     ,bg="white")
barplot(H, main="Shanon index", plot = TRUE, beside=TRUE, width=5, 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),col="black", las=2)


barplot(simp, main="Simpson's index", plot = TRUE, beside=TRUE, width=5, 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),col="black", las=2)


barplot(invsimp,main="Inverse Simpson's index", plot = TRUE, beside=TRUE, width=5,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),col="black", las=2)


barplot(alpha, main=expression(paste("Fisher's logarithmic series ",alpha, " parameter")),plot = TRUE, beside=TRUE, width=5, cex.main=as.numeric(opt$fsize),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),col="black", las=2)


barplot(J, main="Pielou's evenness",plot = TRUE, beside=TRUE, width=5, 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),col="black", las=2)
dev.off()

dI<-data.frame(H,simp,invsimp,alpha,J)
colnames(dI)<-c("Shanon","Simpson","Inverse_Simpson","Fisher_Alpha","Pielou")
write.csv(dI,file=paste(opt$opath,substr(basename(opt$ifile),1,nchar(basename(opt$ifile))-4),"_DI_INDICES.csv",sep=""),quote=F)