inst/doc/BAGS.R

In BAGS: A Bayesian Approach for Geneset Selection

### R code from vignette source 'BAGS.Rnw'

###################################################
### code chunk number 1: install-pkg (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install("BAGS")


###################################################
### code chunk number 2: BAGS.Rnw:89-90 (eval = FALSE)
###################################################
## library(BAGS)


###################################################
### code chunk number 3: BAGS.Rnw:97-98 (eval = FALSE)
###################################################
## library(help=BAGS)


###################################################
### code chunk number 4: BAGS.Rnw:125-166
###################################################
library(BAGS)
library(breastCancerVDX)
library(Biobase)
 
data(vdx)
gene.expr=exprs(vdx)   # Gene expression of the package
vdx.annot=fData(vdx)   # Annotation associated to the dataset
vdx.clinc=pData(vdx)   # Clinical information associated to the dataset 
 
# Identifying the sample identifiers associated to ER+ and ER- breast cancer
er.pos=which(vdx.clinc$er==1)
er.neg=which(vdx.clinc$er==0)
 
# Only keep columns 1 and 3, probeset identifiers and Gene symbols respectively
vdx.annot=vdx.annot[,c(1,3)]

# Checking if the probeset are ordered with respect to the dataset
all(rownames(gene.expr)==as.character(vdx.annot[,1]))  
# Checking if the sample identifiers are order with respect to the dataset
all(colnames(gene.expr)==as.character(vdx.clinc[,1]))  
# Changing the row identifiers to the gene identifiers of interest
rownames(gene.expr)=as.character(vdx.annot[,2])        
 
#= Because we have several measurements for a gene, we filter the genes
# Function to obtain the genes with highest variabilty among phenotypes
gene.nms.u=unique(rownames(gene.expr))
gene.nms=rownames(gene.expr)
indices=NULL
for(i in 1:length(gene.nms.u))
{
	aux=which(gene.nms==gene.nms.u[i])
 	if(length(aux)>1){
 		var.r = apply(cbind(apply(gene.expr[aux,er.pos],1,mean),
                         apply(gene.expr[aux,er.neg],1,mean)),1,var)
		aux=aux[which.max(var.r)]
  	}
  indices=c(indices,aux)
}
#===== Only keep the genes with most variability among the phenotypes of interest
gene.expr=gene.expr[indices,]
gene.nams=rownames(gene.expr)     # The gene symbols of interest are stored here


###################################################
### code chunk number 5: BAGS.Rnw:171-175
###################################################
data(AnnotationMFGO,package="BAGS")
data.gene.grps=DataGeneSets(AnnotationMFGO,gene.nams,5)
phntp.list=list(er.pos,er.neg)
data.mcmc=MCMCDataSet(gene.expr,data.gene.grps$DataGeneSetsIds,phntp.list)


###################################################
### code chunk number 6: BAGS.Rnw:184-200
###################################################
noRow=dim(data.mcmc$y.mu)[1]
noCol=unlist(lapply(phntp.list,length))
iter=10000
GrpSzs=data.gene.grps$Size
YMu=data.mcmc$y.mu
L0=rep(2,2)
V0=rep(4,2)
L0A=rep(3,1)
V0A=rep(3,1)
MM=0.55
AAPi=10
ApriDiffExp=floor(dim(data.mcmc$y.mu)[1]*0.03)
results=matrix(0,noRow,iter)
 		
mcmc.chains=Gibbs2(noRow,noCol,iter,GrpSzs,YMu,L0,V0,L0A,V0A,MM,AAPi,ApriDiffExp,
                   results)


###################################################
### code chunk number 7: BAGS.Rnw:202-211
###################################################
burn.in=2000
alfa.pi=apply(mcmc.chains[[1]][,burn.in:iter],1,function(x){
 	                        y=length(which(x!=0))/length(burn.in:iter);return(y)})

plot(alfa.pi,type="h",main="Probabilities of MF differentially expressed
     between ER status")                                   
cut.off=0.9
abline(h=cut.off,col="red")
differential.processes=names(data.gene.grps$Size)[which(alfa.pi>cut.off)]


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