
runSRMstats<-function(){
#=====================
# install package
# pre-processing data: quality control of MS runs



cat("\n\n =======================================")
cat("\n ** Loading the required library..... \n")


# load the library
library(SRMstats)

# Input data
filename<-commandArgs(trailingOnly=TRUE)[1];

raw<-read.csv(filename)
colnames(raw)[10]<-"Intensity"
raw$Intensity<-as.character(raw$Intensity)
raw$Intensity<-as.numeric(raw$Intensity)

## impute zero to NA
raw[!is.na(raw$Intensity)&raw$Intensity==0,"Intensity"]<-NA

#=====================
# Function: dataProcess
# pre-processing data: quality control of MS runs

cat("\n\n =======================================")
cat("\n ** Data Processing for analysis..... \n")

quantData<-dataProcess(raw)

write.csv(quantData,file="dataProcessedData.csv")
cat("\n Saved dataProcessedData.csv \n")



#=====================
# Function: dataProcessPlots
# visualization 

cat("\n\n =======================================")
cat("\n ** Generating dataProcess Plots..... \n")


dataProcessPlots(data=quantData,type="ProfilePlot",text.size=1.5,address="")
cat("\n Saved ProfilePlot.pdf \n")

dataProcessPlots(data=quantData,type="QCPlot",text.size=1.5,address="")
cat("\n Saved QCPlot.pdf \n")


#=====================
# Function: groupComparison
# generate testing results of protein inferences across concentrations

cat("\n\n =======================================")
cat("\n ** Let's compare Groups..... \n")


## input from groupComparison
# comparison matrix
namelevel<-levels(quantData$GROUP_ORIGINAL)[1]

for(i in 2:length(levels(quantData$GROUP_ORIGINAL))){
	namelevel<-paste(namelevel, levels(quantData$GROUP_ORIGINAL)[i], sep=", ")
}

cat(paste("\n\n Groups are (",namelevel, ")\n",sep=""))

## first, number of comparisons
cat("\n How many comparisons do you want? : ")
numcomparison<-readLines(file("stdin"),1)
numcomparison<-as.integer(numcomparison)

## then, what comparisons?

alllevel<-levels(quantData$GROUP_ORIGINAL)

comparison<-NULL

for(k in 1:numcomparison){

cat(paste("\n Comparison ",k," : Enter constants for each group \n",sep=""))

temp<-NULL

for(i in 1:length(alllevel)){
cat(paste(" ", alllevel[i], " : ",sep=""))
inputtemp<-readLines(file("stdin"),1)
#inputtemp<-readline()
temp<-c(temp,as.numeric(inputtemp))
}

comparison<-rbind(comparison, matrix(temp,nrow=1))

cat(" Name of this comparison : ")
inputrowname<-readLines(file("stdin"),1)
#inputrowname<-readline()
row.names(comparison)[k]<-inputrowname

}

## input for options
cat("\n\n ** Enter your choices for model..... \n")

# labeled or label-free?
cat("\n Label-free(FALSE-F)? or Labeled(TRUE-T)?. Default is TRUE : ")
inputlabel<-readLines(file("stdin"),1)
#inputlabel<-readline()


if(inputlabel=="" | inputlabel=="T") inputlabel<-"TRUE"
if(inputlabel=="F") inputlabel<-"FALSE"
if(inputlabel!="TRUE" & inputlabel!="FALSE"){
	cat("\n Wrong input. will use the default=TRUE.\n")
	inputlabel<-"TRUE"
}


# scope of biological replicate?
cat("\n Scope of Biological Replicate (expanded-E/restricted-R)? Default is R : ")
inputbio<-readLines(file("stdin"),1)
#inputbio<-readline()

if(inputbio=="" | inputbio=="R") inputbio<-"restricted"
if(inputbio=="E") inputbio<-"expanded"
if(inputbio!="expanded" & inputbio!="restricted"){
	cat("\n Wrong input. will use the default=restricted. \n")
	inputbio<-"restricted" 
}


# scope of technical replicate?
cat("\n Scope of Technical Replicate (expanded-E/restricted-R)? Default is E : ")
inputtech<-readLines(file("stdin"),1)
#inputtech<-readline()

if(inputtech=="" | inputtech=="E") inputtech<-"expanded"
if(inputtech=="R") inputtech<-"restricted"
if(inputtech!="expanded" & inputtech!="expanded"){
	cat("\n Wrong input. will use the default=expanded. \n")
	inputtech<-"expanded" 
}



# interference or not?
cat("\n Contain inference transitions (TRUE) or not (FALSE)? : ")
inputinfer<-readLines(file("stdin"),1)
#inputinfer<-readline()

if(inputinfer=="" | inputinfer=="T") inputinfer<-"TRUE"
if(inputinfer=="F") inputinfer<-"FALSE"
if(inputinfer!="TRUE" & inputinfer!="FALSE"){
	cat("\n Wrong input. will use the default=TRUE.\n")
	inputinfer<-"TRUE"
}



## then testing with inputs from users
cat("\n Starting comparison... \n")

resultComparison<-groupComparison(contrast.matrix=comparison,data=quantData,label=inputlabel,scopeOfBioReplication=inputbio, scopeOfTechReplication=inputtech, interference=inputinfer)

if(class(resultComparison)!="try-error"){
	write.csv(resultComparison,"TestingResult.csv")
	cat("\n Saved the testing result. \n")
}else{
	cat("\n Error : Can't analyze. \n")
}

#=====================
# Function: groupComparisonPlots
# visualization for testing results

cat("\n\n =======================================")
cat("\n ** Visualizing testing results..... \n")


# Visualization 1: Volcano plot
# default setup: FDR cutoff = 0.05; fold change cutoff = NA
groupComparisonPlots(data=resultComparison,type="VolcanoPlot",address="")
cat("\n Saved VolcanoPlot.pdf \n")

# Visualization 2: Heatmap (required more than one comparisons)
if(numcomparison>1){
  groupComparisonPlots(data=resultComparison,type="Heatmap",address="")
  cat("\n Saved Heatmap.pdf \n")
} else{
  cat("\n No Heatmap. Need more than 1 comparison for Heatmap. \n")
}


# Visualization 3: Comparison plot
groupComparisonPlots(data=resultComparison,type="ComparisonPlot",address="")
cat("\n Saved ComparisonPlot.pdf \n")


#=====================
# Function: designSampleSize
# calulate number of biological replicates per group for your next experiment

cat("\n\n =======================================")
cat("\n ** Calculating sample size..... \n")

## input for options
# num of sample?
cat("\n Number of Biological Replicates :  ")
inputnsample<-readLines(file("stdin"),1)
#inputnsample<-readline()
if(inputnsample!="") inputnsample<-as.numeric(inputnsample)
if(inputnsample=="") inputnsample<-"TRUE"

# num of peptide?
cat("\n Number of Peptides per protein :  ")
inputnpep<-readLines(file("stdin"),1)
#inputnpep<-readline()
if(inputnpep!="") inputnpep<-as.numeric(inputnpep)
if(inputnpep=="") inputnpep<-"TRUE"

# num of peptide?
cat("\n Number of Transitions per peptides :  ")
inputntran<-readLines(file("stdin"),1)
#inputntran<-readline()
if(inputntran!="") inputntran<-as.numeric(inputntran)
if(inputntran=="") inputntran<-"TRUE"

# power?
cat("\n power : ")
inputpower<-readLines(file("stdin"),1)
#inputpower<-readline()
if(inputpower!="") inputpower<-as.numeric(inputpower)
if(inputpower=="") inputpower<-"TRUE"

# FDR?
cat("\n FDR : (Default is 0.05)")
inputfdr<-readLines(file("stdin"),1)
#inputfdr<-readline()
if(inputfdr=="") inputfdr<-0.05
inputfdr<-as.numeric(inputfdr)

# desiredFC?
cat("\n lower desired Fold Change : (Default is 1.25) ")
inputlower<-readLines(file("stdin"),1)
#inputlower<-readline()
if(inputlower=="") inputlower<-1.25
inputlower<-as.numeric(inputlower)

cat("\n upper desired Fold Change : (Default is 1.75)")
inputupper<-readLines(file("stdin"),1)
#inputupper<-readline()
if(inputupper=="") inputupper<-1.75
inputupper<-as.numeric(inputupper)


result.sample<-designSampleSize(data=quantData,numSample=inputnsample,numPep=inputnpep,numTran=inputntran,desiredFC=c(inputlower,inputupper),FDR=inputfdr,power=inputpower)
#result.sample<-designSampleSize(data=quantData,numSample=TRUE,numPep=2,numTran=3,desiredFC=c(1.25,1.75),FDR=0.05,power=0.8)

if(class(result.sample)!="try-error"){
	write.csv(result.sample,"SampleSizeCalculation.csv")
	cat("\n Saved the Sample Size Calculation. \n")
}else{
	cat("\n Error : Can't analyze. \n")
}



#=====================
# Function: designSampleSizePlots
# visualization for sample size calculation

designSampleSizePlots(data=result.sample, address="")
cat("\n Saved SampleSizePlot.pdf \n")

}


#=====================
tryCatch({runSRMstats()}, 
finally = {
cat("Press enter to finish")
finish<-readLines(file("stdin"),n=1)
})


#cat("\n\n Done. Press enter to finish.")
#finish<-readLines(file("stdin"),1)



