#########################################
## Install MSstats.daily (only need to install once)
#########################################
# 1. Please install all dependancy packages first.
install.packages(c("gplots","lme4","ggplot2","ggrepel","reshape","reshape2","data.table","Rcpp","survival"))
source("http://bioconductor.org/biocLite.R")
biocLite(c("limma","marray","preprocessCore","MSnbase"))

# 2. Install MSstats
install.packages(pkgs = "MSstats.daily_3.3.3.tar.gz", repos = NULL, type ="source") 

# MSstats installation is complete ~~~

#########################################
## Run MSstats (start here if you already install MSstats)
#########################################
# load the library
library(MSstats)

# Help file
?MSstats

#=====================
# Function: dataProcess
# pre-processing data: quality control of MS runs
head(SRMRawData)

QuantData<-dataProcess(SRMRawData)
head(QuantData$ProcessedData)


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

# 1. Profile plot
dataProcessPlots(data=QuantData,type="ProfilePlot",width=20, height=10)

# 2. Quality control plot 
dataProcessPlots(data=QuantData,type="QCPlot")	

# 3. Quantification plot for conditions
dataProcessPlots(data=QuantData,type="ConditionPlot")


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

levels(QuantData$ProcessedData$GROUP_ORIGINAL)

## based on multiple comparisons  (T1 vs T3; T1 vs T7; T1 vs T9)
comparison1<-matrix(c(-1,0,1,0,0,0,0,0,0,0),nrow=1)
comparison2<-matrix(c(-1,0,0,0,0,0,1,0,0,0),nrow=1)
comparison3<-matrix(c(-1,0,0,0,0,0,0,0,1,0),nrow=1)
comparison<-rbind(comparison1,comparison2, comparison3)
row.names(comparison)<-c("T3-T1","T7-T1","T9-T1")

testResultMultiComparisons<-groupComparison(contrast.matrix=comparison,data=QuantData)

testResultMultiComparisons$ComparisonResult


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

# Volcano plot with FDR cutoff = 0.05 and no FC cutoff
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="VolcanoPlot",logBase.pvalue=2,address="Ex1_")

# Volcano plot with FDR cutoff = 0.05, FC cutoff = 70, upper y-axis limit = 100, and no protein name displayed
# FCcutoff=70 is for demonstration purpose
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="VolcanoPlot",FCcutoff=70, logBase.pvalue=2, ylimUp=100, ProteinName=FALSE,address="Ex2_")

# Heatmap with FDR cutoff = 0.05
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="Heatmap", logBase.pvalue=2, address="Ex1_")

# Heatmap with FDR cutoff = 0.05 and FC cutoff = 70
# FCcutoff=70 is for demonstration purpose
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="Heatmap",FCcutoff=70, logBase.pvalue=2, address="Ex2_")

# Comparison Plot
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="ComparisonPlot",address="Ex1_")

# Comparison Plot
groupComparisonPlots(data=testResultMultiComparisons$ComparisonResult,type="ComparisonPlot",ylimUp=8,ylimDown=-1,address="Ex2_")


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

#(1) Sample size calculation
result.sample<-designSampleSize(data=testResultMultiComparisons$fittedmodel,numSample=TRUE,desiredFC=c(1.25,1.75),FDR=0.05,power=0.9)

result.sample

#(2) Power calculation
result.power<-designSampleSize(data=testResultMultiComparisons$fittedmodel,numSample=20,desiredFC=c(1.25,1.75),FDR=0.05,power=TRUE)

result.power


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

pdf("SampleSizePlot.pdf")
designSampleSizePlots(data=result.sample)
dev.off()


#=====================
# Function: quantification
# model-based quantification

# (1) subject quantification

sampleQuant<-quantification(QuantData)
head(sampleQuant)

# (2) Group quantification

groupQuant<-quantification(QuantData, type="Group")
head(groupQuant)