Ejercicio Microarreglos

intro.R

#Analisis de microarreglos

#Instalar programas requeridos
source("http://bioconductor.org/biocLite.R")
biocLite(c("affy","ALL","annotate","hgu95av2.db","multtest","topGO","genefilter"))

#Cargar la libreria ALL
library(ALL)
data(ALL)
#Dimension del objeto
dim(exprs(ALL))

#Definir un objeto ExpressionSet con todas las muestras
table(ALL$BT)
table(ALL$mol.biol)
subset <- intersect(grep("^B", as.character(ALL$BT)),
                    which(as.character(ALL$mol.biol) %in% c("BCR/ABL", "NEG")))
eset <- ALL[, subset]
eset$mol.biol <- factor(eset$mol.biol)
table(eset$mol.biol)

#Filtro 

library(genefilter)
f1 <- pOverA(0.25, log2(100))
f2 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1, f2)
selected <- genefilter(eset, ff)
sum(selected)
esetSub <- eset[selected, ]

library(multtest)
cl <- as.numeric(esetSub$mol.biol == "BCR/ABL")
t <- mt.teststat(exprs(esetSub), classlabel = cl, test = "t.equalvar")
pt <- 2 * pt(-abs(t), df = ncol(exprs(esetSub)) - 2)
hist(pt, 50)

pAdjusted <- mt.rawp2adjp(pt, proc = c("BH"))
sum(pAdjusted$adjp[, "BH"] < 0.1)

pBH <- pAdjusted$adjp[order(pAdjusted$index), "BH"]
names(pBH)<-featureNames(esetSub)

library(annotate)
library(hgu95av2.db)
diff <- pAdjusted$index[1:10]
genesymbolsDiff <- unlist(mget(featureNames(esetSub)[diff], hgu95av2SYMBOL))
genesymbolsDiff

geneSymbols = unlist(mget(featureNames(ALL), hgu95av2SYMBOL))
ABL1probes <- which(geneSymbols == "ABL1")
selected[ABL1probes]

tABL1 <- mt.teststat(exprs(eset)[ABL1probes, ], classlabel = cl,
                     test = "t.equalvar")
ptABL1 <- 2 * pt(-abs(tABL1), df = ncol(exprs(esetSub)) - 2)
sort(ptABL1)

gN <- featureNames(esetSub)
tykin <- unique(unlist(lookUp("GO:0004713", "hgu95av2", "GO2ALLPROBES")))
str(tykin)
sel <- (gN %in% tykin)

tab <- table(pt < 0.05, sel, dnn = c("p < 0.05", "tykin"))
print(tab)
fisher.test(tab)

library(topGO) #load the package

topDiffGenes <- function(allScore) return(allScore < 0.05)

GOdataMF <- new("topGOdata", ontology = "MF",
                allGenes = pBH, geneSel = topDiffGenes, annot = annFUN.db,
                affyLib = "hgu95av2.db")

resultFisher <- runTest(GOdataMF, algorithm = "classic", statistic = "fisher")

resultKS<- runTest(GOdataMF, algorithm = "classic", statistic = "ks")
resultFisher.elim<- runTest(GOdataMF, algorithm = "elim", statistic = "fisher")
resultKS.elim<- runTest(GOdataMF, algorithm = "elim", statistic = "ks")

allRes <- GenTable(GOdataMF, classicFisher = resultFisher,
                   classicKS = resultKS,elimFisher = resultFisher.elim, elimKS = resultKS.elim,
                   orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 20)
allRes

gt <- printGenes(GOdataMF, whichTerms = "GO:0004713", chip = "hgu95av2.db",
                 numChar = 40)
gt