/BTR_resultscripts

Permalink
BTR_resultscripts/krum_wilson_grow_hsc.R
Newer
Older
Raw Normal view History
100644 42 lines (31 sloc) 1.39 KB
Dec 4, 2015 @cheeyeelim Initial commit
1 #(1) Setup paths and environment.
Dec 11, 2015 @cheeyeelim Update name
2 library(BTR)
Dec 4, 2015 @cheeyeelim Initial commit
3 library(doParallel)
4
Mar 16, 2016 @cheeyeelim Changes for revision.
5 path='~/btr_output/'
Dec 4, 2015 @cheeyeelim Initial commit
6 setwd(path)
7
8 #Setting up for parallel processing.
9 registerDoParallel(cores=4)
10
11 #(2)Load data.
12 data(krum_bmodel)
13 data(krum_istate)
14 data(wilson_raw_data)
15
16 bmodel = initialise_model(krum_bmodel)
17 istate = initialise_data(krum_istate)
18
19 #(3)Filtering expression data.
Mar 16, 2016 @cheeyeelim Changes for revision.
20 cdata = initialise_raw_data(wilson_raw_data, max_expr = 'low') #do not filter data at this stage. keep the whole matrix.
Dec 4, 2015 @cheeyeelim Initial commit
21 cell_ind = grepl('cmp', rownames(cdata)) | grepl('gmp', rownames(cdata)) | grepl('mep', rownames(cdata))
22 fcdata = cdata[cell_ind,] #select only relevant cells.
23
24 #(4)Adding extra genes into the model.
25 extra_genes = setdiff(colnames(wilson_raw_data), bmodel@target)
26 grown_bmodel = grow_bmodel(extra_genes[c(19, 20)], bmodel) #genes to be added: ldb1, lmo2
27
28 #(5)Re-estimate initial state.
29 #Since CMPs are upsteam of GMPs, and MEPs, the initial state of extra genes can be estimated from CMPs.
30 tmp_istate = colMeans(cdata[grepl('cmp', rownames(cdata)), extra_genes[c(19, 20)]])
31 tmp_istate = matrix(round(tmp_istate), nrow=1, dimnames=list(1, names(tmp_istate)))
32 grown_istate = cbind(istate, tmp_istate)
33 grown_istate = initialise_data(grown_istate)
34
Mar 16, 2016 @cheeyeelim Changes for revision.
35 result = model_train(cdata=fcdata, bmodel=grown_bmodel, istate=grown_istate, verbose=T)
Dec 4, 2015 @cheeyeelim Initial commit
36
37 filename = 'krum_wilson_myeloid_grow_trained_bmodel.rda'
38 save(result, file=filename)
39
40 #Cleaning up.
41 stopImplicitCluster()