### R code from vignette source 'GeneGroupAnalysis.Rnw'
### Encoding: ISO8859-1

###################################################
### code chunk number 1: install-pkg (eval = FALSE)
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("GeneGroupAnalysis")


###################################################
### code chunk number 2: loadlib (eval = FALSE)
###################################################
## library(GeneGroupAnalysis)


###################################################
### code chunk number 3: GeneGroupAnalysishelp (eval = FALSE)
###################################################
## library(help=GeneGroupAnalysis)


###################################################
### code chunk number 4: GeneGroupAnalysis.Rnw:122-132
###################################################
library(breastCancerVDX)
library(GeneGroupAnalysis)
library(hgu133a.db)
library(annotate)

data(vdx,package="breastCancerVDX")
#-Normalized expression data set
minn.data.expr=exprs(vdx)
#--- Checking that the columns correspond to their respective phenotype data id
#all(colnames(minn.data.expr)==rownames(pData(vdx)))


###################################################
### code chunk number 5: GeneGroupAnalysis.Rnw:151-154
###################################################
#---Checking the annotation of the data
#annotation(vdx)
array.info=ArrayInfoFun(minn.data.expr,hgu133aSYMBOL)


###################################################
### code chunk number 6: GeneGroupAnalysis.Rnw:158-163
###################################################
#---ER status of patients
er.status=pData(vdx)$er
er.pos=which(er.status==1)
er.neg=which(er.status==0)
genes.max.var=GeneMaxVarFun(array.info,minn.data.expr,er.pos,er.neg)


###################################################
### code chunk number 7: GeneGroupAnalysis.Rnw:168-169
###################################################
GO2Gene.grps=GeneGrps2AffyGrpsFun("CC",3,array.info$genes.name.unique,genes.max.var$IndexMaxVar)


###################################################
### code chunk number 8: GeneGroupAnalysis.Rnw:174-175
###################################################
Minn07WrkngGrps=SizeGOAffyGrps(GO2Gene.grps$index.GO.grps,100)


###################################################
### code chunk number 9: GeneGroupAnalysis.Rnw:179-181
###################################################
Minn07MCMCData=MCMCData.cs(Minn07WrkngGrps$groups,GO2Gene.grps$index.GO.grps,
GO2Gene.grps$GO.grps,minn.data.expr,er.pos,er.neg)


###################################################
### code chunk number 10: GeneGroupAnalysis.Rnw:190-218
###################################################
nsim=40       
burn.in=10   

v.SS.A.i=matrix(0,length(Minn07MCMCData$proc.GO),nsim)
v.SS.B.i=matrix(0,length(Minn07MCMCData$proc.GO),nsim)
v.SS.Al=rep(0,nsim)
v.beta.i=matrix(0,length(Minn07MCMCData$proc.GO),nsim)
v.alfa.i=matrix(0,length(Minn07MCMCData$proc.GO),nsim)
v.pi.a.i=matrix(0,length(Minn07MCMCData$proc.GO),nsim)
v.rho.a=rep(0,nsim)
#----- Initialization of the parameters
v.beta.i[,1]=rnorm(length(Minn07MCMCData$proc.GO),0,1)
v.alfa.i[,1]=rnorm(length(Minn07MCMCData$proc.GO),0,1)
v.SS.Al[1]=0.1   
L0.alfa=0.001
v.pi.a.i[,1]=runif(length(Minn07MCMCData$proc.GO))
v.rho.a[1]=0.1
Grps.apriori.diff.exp=23 
shape=3
scale=0.3
mm.pi=0.75
aa.pi=10
often=20
cut.off=0.7

result.MCMC=GeneGroupAnalysis::GibbsFun.cs(Minn07MCMCData$y.mu.a,Minn07MCMCData$y.mu.b,
Minn07WrkngGrps$group.size,v.beta.i,v.alfa.i,v.SS.A.i,v.SS.B.i,v.SS.Al,v.rho.a,v.pi.a.i,mm.pi,shape,
scale,aa.pi,Grps.apriori.diff.exp,nsim,burn.in,often,cut.off)


###################################################
### code chunk number 11: GeneGroupAnalysis.Rnw:226-235
###################################################
#-------------------------------------------------
#- Executable example after the application of GibbsFun.cs.
#-------------------------------------------------
total=length(burn.in:nsim)
prob.a=apply(result.MCMC$Alfa[,burn.in:nsim],1,function(x) length(which(x!=0))/total)
no.groups=length(Minn07MCMCData$proc.GO)
plot(1:no.groups,prob.a,type="h",main="Probabilities of Alpha different to 0",
xlab="CC processes considered",ylab="probability of being differentially expressed")
abline(h=cut.off,lwd=2,col="red")


###################################################
### code chunk number 12: GeneGroupAnalysis.Rnw:249-257
###################################################
library(GeneGroupAnalysis)
library(annotate)
library(rheumaticConditionWOLLBOLD)
library(hgu133plus2.db)

data(wollbold,package="rheumaticConditionWOLLBOLD")
#----Normalized expression data set
woll.data.exp=exprs(wollbold)


###################################################
### code chunk number 13: GeneGroupAnalysis.Rnw:282-285
###################################################
#---Checking the annotation of the data
#annotation(wollbold)
array.info=ArrayInfoFun(woll.data.exp,hgu133plus2SYMBOL)


###################################################
### code chunk number 14: GeneGroupAnalysis.Rnw:290-294
###################################################
#---ER status of patients
TGF.exp=1:30
TNF.exp=31:60
genes.max.var=GeneMaxVarFun(array.info,woll.data.exp,TGF.exp,TNF.exp)


###################################################
### code chunk number 15: GeneGroupAnalysis.Rnw:299-300
###################################################
GO2Gene.grps=GeneGrps2AffyGrpsFun("MF",3,array.info$genes.name.unique,genes.max.var$IndexMaxVar)


###################################################
### code chunk number 16: GeneGroupAnalysis.Rnw:305-306
###################################################
Wollbold09WrkngGrps=SizeGOAffyGrps(GO2Gene.grps$index.GO.grps,100)


###################################################
### code chunk number 17: GeneGroupAnalysis.Rnw:311-314
###################################################
indexes.1=c(1,2,3,4,5,6)
Wollbold09MCMCData=MCMCData.ts(Wollbold09WrkngGrps$groups,GO2Gene.grps$index.GO.grps,
GO2Gene.grps$GO.grps,woll.data.exp,TGF.exp,TNF.exp,indexes.1,5)


###################################################
### code chunk number 18: GeneGroupAnalysis.Rnw:326-361
###################################################
nsim=40
burn.in=10
Grps.apriori.diff.exp=23
shape=3
scale=0.1
mm.pi=0.75
aa.pi=10
often=20
cut.off=0.72
df.lambda=10
no.time.pnts=5
no.pcnts=6

v.SS.i=array(0,c(length(Wollbold09MCMCData$proc.GO),no.time.pnts,no.time.pnts,2))
v.SS.Al=array(0,c((no.time.pnts-1),(no.time.pnts-1),2))
v.beta.i=array(0,c(length(Wollbold09MCMCData$proc.GO),no.time.pnts,nsim))
v.alfa.0.i=matrix(0,length(Wollbold09MCMCData$proc.GO),nsim)
v.alfa.i=array(0,c(length(Wollbold09MCMCData$proc.GO),(no.time.pnts-1),nsim))
v.pi.a.i=matrix(0,length(Wollbold09MCMCData$proc.GO),nsim)
v.rho.a=rep(0,nsim)
v.beta.i[,,1]=rnorm(length(Wollbold09MCMCData$proc.GO),0,1)
v.alfa.i[,,1]=rnorm(length(Wollbold09MCMCData$proc.GO),0,1)
v.alfa.0.i[,1]=rnorm(length(Wollbold09MCMCData$proc.GO),0,1)
v.SS.Al[,,1]=diag(0.1,(no.time.pnts-1))
v.pi.a.i[,1]=runif(length(Wollbold09MCMCData$proc.GO))
v.rho.a[1]=0.1                
L0.alfa=diag(0.005,(no.time.pnts-1))

indexes.1=c(0,6,12,18,24) 
indexes.2=c(0,6,12,18)  

results.MCMC=GeneGroupAnalysis::GibbsNPFun.ts(Wollbold09MCMCData$y.mu.a,
Wollbold09MCMCData$y.mu.b,Wollbold09WrkngGrps$group.size,v.beta.i,v.alfa.0.i,
v.alfa.i,v.SS.Al,v.rho.a,v.pi.a.i,mm.pi,aa.pi,Wollbold09MCMCData$lambda,df.lambda,L0.alfa,
indexes.1,indexes.2,no.time.pnts,no.pcnts,Grps.apriori.diff.exp,nsim,burn.in,often,cut.off)


###################################################
### code chunk number 19: GeneGroupAnalysis.Rnw:372-381
###################################################
#---------------------------------------------------
#- Executable example after the application of GibbsNPFun.ts
#---------------------------------------------------
total=length(burn.in:nsim)
prob.a=apply(results.MCMC$Alfa[,1,burn.in:nsim],1,function(x) length(which(x!=0))/total)
no.groups=length(Wollbold09MCMCData$proc.GO)
plot(1:no.groups,prob.a,type="h",main="Probabilities of Alpha different to 0",
xlab="CC processes considered",ylab="probability of being differentially expressed")
abline(h=cut.off,lwd=2,col="red")


###################################################
### code chunk number 20: GeneGroupAnalysis.Rnw:387-397
###################################################
#---------------------------------------------------
#- Executable example after the application of GibbsNPFun.ts
#---------------------------------------------------
library(annotate)
library(GO.db)
selection=which(prob.a>=cut.off)
selected.GO.groups=GO2Gene.grps$GO.grps[Wollbold09WrkngGrps$groups][selection]
print(selected.GO.groups)
funct.classes=getGOTerm(selected.GO.groups)[[1]]
print(funct.classes)


###################################################
### code chunk number 21: sessionInfo
###################################################
toLatex(sessionInfo())


