#########################################
## 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
#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("Rat_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)
unique(raw$PeptideSequence)


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

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

quantData<-dataProcess(raw)

quantData<-dataProcess(raw, normalization=FALSE)

quantData[1:5,]


#=====================
# Function: dataProcessPlots
# visualization 
? dataProcessPlots
dataProcessPlots(data=quantData.q,type="ProfilePlot",text.size=1.2,width=10,height=8,address="fig_rat/")
dataProcessPlots(data=testquantile,type="QCPlot",text.size=1.2,width=10, height=8,address="fig_rat/")


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

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

## fixed Subject, without interaction
result.fS<-groupComparison(contrast.matrix=comparison,data=quantData,labeled=FALSE,scopeOfTechReplication="restricted",interference=FALSE)

## random Subject, without interaction
result.rS<-groupComparison(contrast.matrix=comparison,data=quantData,scopeOfBioReplication="expanded",labeled=FALSE, interference=FALSE)

## save the result
write.csv(result.fS, file="fig_rat/result_fixedSub.csv")



#=====================
# 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.rS,type="VolcanoPlot",ProteinName=FALSE,address="fig_rat/randomSub_")

# (2) Both FDR cutoff = 0.05; fold change cutoff = 1.2
groupComparisonPlots(data=result.fS.noInfer,type="VolcanoPlot",ProteinName=FALSE,FCcutoff=1.2,address="fig_rat/randomSub_FC_")


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

# with fixed Subject
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


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

designSampleSizePlots(data=result.sample,address="fig_rat/Rat_Sample_")


#########################################
## Part 4: poor quality peaks or inconsistent peptides
#########################################

#### 1. remove the feature entirely,
raw1<-raw[raw$PeptideSequence!="CSSLLWAGAAWLR",] ## NP_001007697
raw1<-raw1[raw1$PeptideSequence!="TGTNLMDFLSR",]  ## NP_037244

quantData_exclude<-dataProcess(raw1)

dataProcessPlots(data=quantData_exclude,type="ProfilePlot",text.size=1.2,width=10,height=8,address="fig_rat/exclude_")


## fixed Subject
result.fS.exclude<-groupComparison(contrast.matrix=comparison,data=quantData_exclude,labeled=FALSE,scopeOfTechReplication="restricted",interference=FALSE)

## random Subject
result.rS.exclude<-groupComparison(contrast.matrix=comparison,data=quantData_exclude,scopeOfBioReplication="expanded",scopeOfTechReplication="restricted",labeled=FALSE, interference=FALSE)

write.csv(result.fS.exclude, file="fig_rat/result_fixedSub_exclude.csv")

groupComparisonPlots(data=result.fS.exclude,type="VolcanoPlot",ylimUp=50,ProteinName=FALSE,address="fig_rat/fixedSub_exclude_")

