mnp_training/training.R at master · mwsill/mnp_training · GitHub

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
80
81
82
83
84
85
86
87
88
89
90
91
#-----------------------------------------------------------------------------------
# training Random Forest classifier
#
# Martin Sill
# m.sill@dkfz.de
#
# Note, in this example we reduce the number of features to 20k probes by sd filtering before applying the random forest for
# feature selection. To perform feature selection as described in the paper remove line 34,
# this will increase the computation time significantly.
#
# 2018-03-14 UTC
#------------------------------------------------------------------------------------
#options(max.print = 1000)
#options(stringsAsFactors = FALSE)
#options(scipen = 999)
rm(list=ls())

library(randomForest)
library(parallel)

ntrees <- 500  # 10000 in the paper, here 500 to speed up the example
cores <- 4
seed <- 180314
p <- 10000

message("loading preprocessed data ...",Sys.time())
load(file.path("results","betas.ba.RData"))

message("performing variable selection ...",Sys.time())
source(file.path("R","train.R"))
y <- as.factor(anno$`methylation class:ch1`)

# sd pre filtering to 20k probes, to speed up the example
betas <- betas[,order(-apply(betas,2,sd))[1:20000]]

set.seed(seed,kind ="L'Ecuyer-CMRG")
message("seed: ",seed)
message("cores: ",cores)
message("ntrees: ",ntrees)
message("n: ",nrow(betas))
message("p: ",ncol(betas))

rf.varsel <- rfp(betas,
                 y,
                 mc=cores,
                 ntree=ntrees,
                 sampsize=rep(min(table(y)),length(table(y))),
                 importance=TRUE)

# get permutation variable importance
imp.meandecrease <- rf.varsel$importance[,dim(rf.varsel$importance)[2]-1]

# save selection forest
save(rf.varsel,file=file.patht("results","varsel.RData"))
rm(rf.varsel)

# reduce data matrix
or <- order(imp.meandecrease,decreasing=T)
betasy <- betas[,or[1:p]]

gc()

message("finished ...",Sys.time())

message("training classifier ...",Sys.time())

message("single core")
message("ntrees: ",ntrees)
message("n: ",nrow(betasy))
message("p: ",ncol(betasy))

rf.pred <- randomForest(betasy,
                        y,
                        #mc=cores,
                        ntree=ntrees,
                        #strata=y,
                        #mtry=sqrt(ncol(betas)),
                        sampsize=rep(min(table(y)),length(table(y))),
                        proximity=TRUE,
                        oob.prox=TRUE,
                        importance=TRUE,
                        keep.inbag=TRUE,
                        do.trace=FALSE,
                        seed=seed
)

message("finished ...",Sys.time())

save(rf.pred,file=file.path("results","rf.pred.RData"))

message("finished ...",Sys.time())