This repository has been archived by the owner on Apr 11, 2023. It is now read-only.
/
bfrmTesting2.R
82 lines (63 loc) · 2.51 KB
/
bfrmTesting2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# bfrmTesting2.R
require(synapseClient)
require(bfrm)
require(affy)
require(parallel)
require(hgu133a.db)
require(corpcor)
## READ IN BETA CATENIN DATA
egfrEnt <- downloadEntity(138511)
egfrExpr <- ReadAffy(filenames = file.path(egfrEnt$cacheDir, egfrEnt$files))
## NORMALIZE AND SUMMARIZE USING RMA
egfrExprNorm <- rma(egfrExpr, normalize = T, background = F)
egfrMat <- exprs(egfrExprNorm)
egfrSamp <- sapply(strsplit(egfrEnt$files, "_", fixed = T), "[[", 3)
egfrVec <- as.numeric(grepl("wtegf", egfrEnt$files, ignore.case = T))
## BUILD THE H MATRIX
# Find the AFFX control probesets
affxInd <- grep('AFFX', rownames(egfrMat), ignore.case = F)
controlMat <- egfrMat[affxInd, ]
controlSVD <- fast.svd(controlMat)
hRightNine <- controlSVD$v[ , 1:9]
# hMatrix <- cbind(rep(1, 19), egfrVec, hRightFive)
## RUN BFRM IN THE NONEVOLUTIONARY SPARSE ANOVA MODE
foo <- bfrm(egfrMat, design = egfrVec, control = t(hRightNine))
mPPib <- foo@results$mPostPib
topProbeLogical <- mPPib[ , 2] >= 0.99
topProbeInd <- grep("TRUE", topProbeLogical)
topProbeIDs <- rownames(egfrMat)[topProbeInd]
topEGFRSyms <- as.character(mget(topProbeIDs, hgu133aSYMBOL,
ifnotfound = NA))
## BRING IN THE GGHEAT CODE ENTITY
ggheatEnt <- loadEntity(274063)
attach(ggheatEnt)
## VISUALIZING THE EGFR DATA SUBSETTED TO THE ANOVA SELECTED INDICES
anovaHeatPlot <- ggheat2(egfrMat[topProbeInd, ], clustering = 'both')
## RUN BFRM IN EVOLUTIONARY MODE
# In Asynchronous mode using {parallel}
# asyncJob <- parallel(fooFactor <- evolve(egfrMat,
# init = as.numeric(topProbeInd),
# maxVarIter = 30,
# maxVars = length(topProbeInd))
# )
# In conventional mode
fooFactor <- evolve(egfrMat,
init = as.numeric(topProbeInd),
control = t(hRightNine),
priorpsia = 2,
priorpsib = 0.005,
varThreshold = 0.85,
facThreshold = 0.95,
maxVarIter = 30,
minFacVars = 10,
maxFacVars = length(topProbeInd),
maxFacs = 50,
maxVars = length(topProbeInd)
)
# OK, this seems to work. We therefore will need to change the default values
# for priorpsia and priorpsib
# Now to evaluate the 'projection()' function:
projFoo <- projection(fooFactor, egfrMat)
par(mfrow = c(2,1))
image(fooFactor@results$mF)
image(t(projFoo)) # The projection result is transposed