Pigengene_inference.R
In Pigengene: Infers biological signatures from gene expression data

## ----style, eval=TRUE, echo=FALSE, results="asis"--------------------------
BiocStyle::latex()

## ----loading, fig.width=6, fig.height=6, echo=TRUE-------------------------
library(Pigengene)
data(aml)
data(mds)
d1 <- rbind(aml,mds)
Labels <- c(rep("AML",nrow(aml)),rep("MDS",nrow(mds)))
names(Labels) <- rownames(d1)
Disease <- as.data.frame(Labels)
h1 <- pheatmap.type(d1[,1:20],annRow=Disease,show_rownames=FALSE)

## ----oneStep, echo=TRUE----------------------------------------------------
p1 <- one.step.pigengene(Data=d1,saveDir='pigengene', bnNum=0, verbose=1,
      seed=1, Labels=Labels, toCompact=FALSE, doHeat=FALSE)

## ----tree, fig.width=5, fig.height=5, echo=TRUE----------------------------
plot(p1$c5treeRes$c5Trees[["34"]])

## ----pigengene, fig.width=5, fig.height=5, echo=TRUE-----------------------
dim(p1$pigengene$eigengenes)
p1 <- pheatmap.type(p1$pigengene$eigengenes,annRow=Disease,show_rownames=FALSE)

## ----qc, echo=TRUE---------------------------------------------------------
## QC:
checked <- check.pigengene.input(Data=d1, Labels=Labels)
DataI <- checked$Data
LabelsI <- checked$Labels

## ----balance, echo=TRUE----------------------------------------------------
wData <- balance(Data=DataI, Labels=LabelsI, verbose=1)$balanced

## ----beta, echo=TRUE-------------------------------------------------------
betaI <- calculate.beta(RsquaredCut=0.8, Data=wData, verbose=1)$power
saveDir <- "steps" ## Results will be saved in this folder.
dir.create(saveDir)

## ----wgcna, echo=TRUE------------------------------------------------------
## WGCNA
wgRes <- wgcna.one.step(Data=wData, seed=1, power=betaI,
      saveDir=saveDir, verbose=1)

## ----eigengenes, echo=TRUE-------------------------------------------------
## Eigengenes:
pigengene <- compute.pigengene(Data=DataI, Labels=LabelsI,
	     	saveFile=file.path(saveDir, 'pigengene.RData'),
                modules=wgRes$modules, verbose=1)
class(pigengene)
print(dim(pigengene$eigengenes)) ##This is the eigengenes matrix.
print(pigengene$eigengenes[1:3,1:4])

## ----bn, echo=TRUE---------------------------------------------------------
## Learning the BN structure:
library(bnlearn)
learnt <- learn.bn(pigengene=pigengene, bnPath=file.path(saveDir, "bn"),
       	  	bnNum=10, ## In real applications, at least 100-1000.
		seed=1, verbose=1)
BN <- learnt$consensus1$BN

## ----draw, echo=TRUE, eval=FALSE-------------------------------------------
#  drawn <- draw.bn(BN)
#  ## By construction, the Disease node can have no parents.

## ----bnParam, echo=TRUE----------------------------------------------------
## Fit the parameters of the Bayesian network:
fit <- bn.fit(x=BN, data=learnt$consensus1$Data, method="bayes", iss=10)
##where learnt$consensus1$Data is the discretized data matrix.

## The conditional probability table for one of the children of the Disease node:
selectedFeatures <- children("Disease", x=BN)
print(fit[[selectedFeatures[1]]])

## ----bnPrediction, echo=TRUE-----------------------------------------------
l2 <- predict(object=fit, node="Disease", data=learnt$consensus1$Data, method="bayes-lw")
table(LabelsI, l2)

## ----trees, echo=TRUE------------------------------------------------------
## Decision trees:
treePath <- file.path(saveDir, 'C5Trees')
dir.create(path=treePath)
treeRes <- make.decision.tree(pigengene=pigengene, Data=DataI,
	       selectedFeatures=selectedFeatures, saveDir=treePath,
       	       verbose=1, toCompact=FALSE)

## ----citation, results='asis', eval=TRUE-----------------------------------
citation("Pigengene")

## ----sessionInfo, results='asis', eval=TRUE--------------------------------
toLatex(sessionInfo())