
############### R package and example data loading  ###############

# source("http://bioconductor.org/biocLite.R")
# biocLite("qvalue")

# source("http://bioconductor.org/biocLite.R")
# biocLite("sva")

library(qvalue)
library(sva)

# getwd()
temp="D:/My Documents/C_language/TBI_report/workshop/"
setwd(temp)

load(file=paste(temp,"SVA.pheno.RData",sep=""))
load(file=paste(temp,"SVA.exp.preprocessed.RData",sep=""))
load(file=paste(temp,"exp.annotation.RData",sep=""))

SVA.pheno
head(SVA.exp.preprocessed)
head(exp.annotation)

############### check subgroup within BRCA2 samples  ###############


BRCA1=which(SVA.pheno[,"Mutation"]=="BRCA1")
BRCA2=which(SVA.pheno[,"Mutation"]=="BRCA2")

mod=model.matrix(~as.factor(Mutation),data=SVA.pheno)

mod0=model.matrix(~1,data=SVA.pheno)

svaObj=sva(log2(SVA.exp.preprocessed[,-1:-3]),mod,mod0)

symbol<-rep(2,15)
symbol[BRCA2]<-3
svaObj$sv

plot(1:15,svaObj$sv[,1],xlab="Array",ylab="Covariate Value",pch=symbol)
legend(4,-0.4,c("BRCA1","BRCA2"),pch=c(2,3),bty = "n")


############### compare q-value before/after batch effect  ###############


pValues.before=f.pvalue(log2(SVA.exp.preprocessed[,-1:-3]),mod,mod0)
qValuesObj.before<-qvalue(pValues.before)

qsummary(qValuesObj.before)

modSv=cbind(mod,svaObj$sv[,1])
mod0Sv=cbind(mod0,svaObj$sv[,1])

pValues.after=f.pvalue(log2(SVA.exp.preprocessed[,-1:-3]),modSv,mod0Sv)
qValuesObj.after<-qvalue(pValues.after)

qsummary(qValuesObj.after)


############### a example (eukaryotic translation factor 2 (EIFS2))  ###############


gene.no<-43

exp.annotation[gene.no,]

plot(1:15,log2(SVA.exp.preprocessed[gene.no,-1:-3]),xlab="Array",ylab="log2(Expression)",pch=symbol)
legend(2,-1.5,c("BRCA1","BRCA2"),pch=c(2,3),bty = "n")


pValues.before[gene.no]
qValuesObj.before$qvalue[gene.no]


pValues.after[gene.no]
qValuesObj.after$qvalue[gene.no]





############### temp  ###############




# pValues=f.pvalue(log2(SVA.exp.preprocessed[,-1:-3]),mod,mod0)
pValues=f.pvalue(as.matrix(SVA.exp.preprocessed),mod,mod0)



qValuesObj<-qvalue(pValues)

qValuesObj$qvalue[1:10]



qsummary(qValuesObj)


SVA.exp.preprocessed[1:5,]



modSv=cbind(mod,sva.result$sv[,1])

mod0Sv=cbind(mod0,sva.result$sv[,1])



# lm(log2(SVA.exp.preprocessed[1,-1:-3])~sva.result$sv[,1])
# SVA.exp.preprocessed.adjusted<-t(sapply(1:3170,function(gene) residuals(lm(log2(SVA.exp.preprocessed[gene,-1:-3])~sva.result$sv[,1]))))






# pValues=f.pvalue(log(SVA.exp.preprocessed[,-1:-3]),modSv,mod0Sv)
pValues=f.pvalue(as.matrix(SVA.exp.preprocessed),modSv,mod0Sv)


qValues<-qvalue(pValues)


qsummary(qValues)




qValues$qvalues[43]


SVA.data[43,]


plot(1:15,SVA.exp.preprocessed[46,-1:-3],xlab="Array",ylab="log2(Expression)",pch=symbol)
legend(locator(1),c("BRCA1","BRCA2"),pch=c(2,3),bty = "n")









# Calibrated.Ratios<-c("BRCA1@1996","BRCA1@1822","BRCA1@1714","BRCA1@1224","BRCA1@1252","BRCA1@1510","BRCA1@1905",
# "BRCA2@1900","BRCA2@1787","BRCA2@1721","BRCA2@1486","BRCA2@1816","BRCA2@1616","BRCA2@1063","BRCA2@1936",
# "Sporadic@1572","Sporadic@1324","Sporadic@1649","Sporadic@1320","Sporadic@1542","Sporadic@1281","Sporadic@1321")

SVA.data=read.table(file=paste(temp,"SVA_example2.data",sep=""),sep="\t")


SVA.data.temp=read.table(file=paste(temp,"SVA_example3.data",sep=""),sep="\t",head=T)[,-12:-18]

dim(SVA.data.temp)












dim(SVA.data)

SVA.exp=matrix(as.numeric(as.matrix(SVA.data[4:3229,4:21])),3226,18)

colnames(SVA.exp)=sapply(SVA.data[3,4:21],as.character)

head(SVA.exp)

SVA.pheno<-data.frame(t(SVA.data[1:2,7:21]))

colnames(SVA.pheno)<-SVA.data[1:2,1]

BRCA1=which(SVA.pheno[,"Mutation"]=="BRCA1")
BRCA2=which(SVA.pheno[,"Mutation"]=="BRCA2")



eliminated<-which(rowSums(SVA.exp[,-1:-3]>20)>=1)


exp.annotation<-SVA.data[c(-1:-3,-eliminated-3),1:6]

exp.annotation<-SVA.data[c(-1:-3,-eliminated-3),1:6]


colnames(exp.annotation)<-sapply(SVA.data[3,1:6],as.character)
save()

SVA.exp.preprocessed<-SVA.exp[c,]

SVA.exp.preprocessed<-SVA.data.temp[-eliminated,-1]

dim(SVA.exp.preprocessed)

exp.annotation<-exp.annotation[,1:3]

save(exp.annotation, file = "exp.annotation.RData")
save(SVA.exp.preprocessed, file = "SVA.exp.preprocessed.RData")
save(SVA.pheno, file = "SVA.pheno.RData")

load(file=paste(temp,"SVA.exp.preprocessed.RData",sep=""))


dim(SVA.exp.preprocessed)


# SVA.example.pheno<-data.frame(t(matrix(unlist(strsplit(Calibrated.Ratios,"@")),2,22))[1:15,])
# rownames(SVA.example.pheno)<-Calibrated.Ratios[1:15]
# colnames(SVA.example.pheno)<-c("sample","mutation")

# save(SVA.example.pheno, file = "SVA.example.pheno.RData")
load(file=paste(temp,"SVA.pheno.RData",sep=""))


BRCA1=which(SVA.pheno[,"Mutation"]=="BRCA1")
BRCA2=which(SVA.pheno[,"Mutation"]=="BRCA2")


# p.value=apply(SVA.example.preprocessed,1,function(data) t.test(as.numeric(data[BRCA1]),as.numeric(data[BRCA2]))[3]$p.value)

# q.value=qvalue(p.value)


# qsummary(q.value)




# SVA.exp.preprocessed.norm<-t(apply(log2(SVA.exp.preprocessed[,-1:-3]),1,function(data) (data-mean(data))/sd(data)))
SVA.exp.preprocessed.norm<-t(apply(SVA.exp.preprocessed,1,function(data) (data-mean(data))/sd(data)))



mod=model.matrix(~as.factor(Mutation),data=SVA.pheno)


mod0=model.matrix(~1,data=SVA.pheno)


# sva.result=sva(as.matrix(SVA.data.temp[,-1]),mod,mod0,method="irw",B=100)


sva.result=sva(as.matrix(SVA.exp.preprocessed),mod,mod0)
# sva.result=irwsva.build(SVA.exp.preprocessed.norm,mod,mod0,n.sv=4,B=100)

symbol<-rep(2,15)
symbol[BRCA2]<-3
sva.result$sv

plot(1:15,sva.result$sv[,1],xlab="Array",ylab="Covariate Value",pch=symbol)
# legend(locator(1),c("BRCA1","BRCA2"),pch=c(2,3),bty = "n")





# pValues=f.pvalue(log2(SVA.exp.preprocessed[,-1:-3]),mod,mod0)
pValues=f.pvalue(as.matrix(SVA.exp.preprocessed),mod,mod0)



qValuesObj<-qvalue(pValues)

qValuesObj$qvalue[1:10]



qsummary(qValuesObj)


SVA.exp.preprocessed[1:5,]



modSv=cbind(mod,sva.result$sv[,1])

mod0Sv=cbind(mod0,sva.result$sv[,1])



# lm(log2(SVA.exp.preprocessed[1,-1:-3])~sva.result$sv[,1])
# SVA.exp.preprocessed.adjusted<-t(sapply(1:3170,function(gene) residuals(lm(log2(SVA.exp.preprocessed[gene,-1:-3])~sva.result$sv[,1]))))






# pValues=f.pvalue(log(SVA.exp.preprocessed[,-1:-3]),modSv,mod0Sv)
pValues=f.pvalue(as.matrix(SVA.exp.preprocessed),modSv,mod0Sv)


qValues<-qvalue(pValues)


qsummary(qValues)

exp.annotation[43,]


qValues$qvalues[43]


SVA.data[43,]


plot(1:15,SVA.exp.preprocessed[46,-1:-3],xlab="Array",ylab="log2(Expression)",pch=symbol)
legend(locator(1),c("BRCA1","BRCA2"),pch=c(2,3),bty = "n")






