### R code from vignette source 'vignettes/GeneGroupAnalysis/inst/doc/GeneGroupAnalysis.Rnw' ### Encoding: UTF-8 ################################################### ### 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())