#########################################
## Part 1: Install SRMstats (only need to install once)
#########################################
# 1. Please install all dependancy packages first.
install.packages(c("gplots","lme4","lattice"))
source("http://bioconductor.org/biocLite.R")
biocLite(c("limma","marray"))

# 2. Install SRMstats or MSstats2
#install.packages(pkgs = "location/SRMstats.tar.gz", repos = NULL, type ="source")
#install.packages(pkgs = "location/SRMstats.tar.gz", repos = NULL, type ="source") 

# SRMstats installation is complete ~~~

# Please exit R and re-enter R again. 
# Once re-enter R, please start with the below commands to exercute SRMstats.


#########################################
## Part 3: Run SRMstats (start here if you already install SRMstats)
#########################################
# load the library
library(SRMstats)
#library(MSstats2)

# Help file
?SRMstats

# Input data
setwd("/Users/Meena/Desktop/UShupo2013/instruction")
raw<-read.csv("Human_plasma.csv")
head(raw)
dim(raw)
colnames(raw)[10]<-"Intensity"
raw$Intensity<-as.character(raw$Intensity)
raw$Intensity<-as.numeric(raw$Intensity)

## remove reference peptides
unique(raw$ProteinName)
raw<-raw[raw$ProteinName!="Tr_peps_set2",]

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

? dataProcess
quantData<-dataProcess(raw)
quantData[1:5,]


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

? dataProcessPlots
dataProcessPlots(data=quantData,type="ProfilePlot",text.size=1.2,width=10,height=8,address="fig_OV/")
dataProcessPlots(data=quantData,type="QCPlot",text.size=1.2,width=10, height=8,address="fig_OV/")


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

?groupComparison
levels(quantData$GROUP_ORIGINAL)
comparison<-matrix(c(-1,1),nrow=1)
row.names(comparison)<-"Disease-Healthy"

# random Run, fixed Subject
result.rRfS<-groupComparison(contrast.matrix=comparison,data=quantData)

# save the result
write.csv(result.rRfS, file="testing_result_rRfS.csv")

# with other options
# fixed Run, fixed Subject
result.fRfS<-groupComparison(contrast.matrix=comparison,scopeOfTechReplication="restricted",data=quantData)

# random Run, random Subject
result.rRrS<-groupComparison(contrast.matrix=comparison,scopeOfBioReplication="expanded",data=quantData)

# fixed Run, random Subject
result.fRrS<-groupComparison(contrast.matrix=comparison,scopeOfTechReplication="restricted",scopeOfBioReplication="expanded",data=quantData)

# fixed Subject, label-free
result.rRfS.LF<-groupComparison(contrast.matrix=comparison,data=quantData,label=FALSE)


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

?groupComparisonPlots

# Visualization 1: Volcano plot
# (1) default setup: FDR cutoff = 0.05; fold change cutoff = NA
groupComparisonPlots(data=result.rRfS,type="VolcanoPlot",address="fig_OV/randomRfixedS_")

# (2) Both FDR cutoff = 0.05; fold change cutoff = 1.2
groupComparisonPlots(data=result.rRfS,type="VolcanoPlot",FCcutoff=1.2,ylimUp=50,xlimUp=2,address="fig_OV/randomRfixedS_wFC_")


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

?designSampleSize

result.sample<-designSampleSize(data=quantData,numSample=TRUE,numPep=2,numTran=3,desiredFC=c(1.05,1.3),FDR=0.05,power=0.8)
result.sample

result.sample.randomS<-designSampleSize(data=quantData,numSample=TRUE,numPep=2,numTran=3,desiredFC=c(1.05,1.3),FDR=0.05,power=0.8,scopeOfBioReplication="expanded")
result.sample.randomS


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

designSampleSizePlots(data=result.sample,address="fig_OV/OV_Sample_")



#########################################
## Part 4: poor quality peaks
#########################################

#### 1. remove the feature entirely,
raw1<-raw[raw$PeptideSequence!="IHWESASLLR",] ## CO3

## data Processing
quantData_exclude<-dataProcess(raw1)

## Profile plot
dataProcessPlots(data=quantData_exclude,type="ProfilePlot",text.size=1.2,width=10,height=8,address="fig_OV/excludeCO3_")

## Comparison with random Run and fixed Subject
result.rRfS.exclude<-groupComparison(contrast.matrix=comparison,data=quantData_exclude)



#### 2. replace the intensities with missing values,

raw2<-raw
raw2[raw2$PeptideSequence=="IHWESASLLR" & !is.na(raw2$Intensity) & raw2$Intensity<2^16,"Intensity"]<-NA

## data Processing
quantData_missing<-dataProcess(raw2)

## Profile plot
dataProcessPlots(data=quantData_missing,type="ProfilePlot",text.size=1.2,width=10,height=8,address="fig_OV/missingCO3_")

## Comparison with random Run and fixed Subject
result.rRfS.missing<-groupComparison(contrast.matrix=comparison,data=quantData_missing)


#### 3. with missing value, remove the interaction

result.rRfS.missing.noInter<-groupComparison(contrast.matrix=comparison,data=quantData_missing,interference=FALSE)


