1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# compact_data = readRDS("/storage11_7T/fuy/TADA-A/annotation/spliceai/2021-03-15_spliceai_DS_AG_6788SNV_sf_uni_prior_compact.rds")
compact_data = readRDS("/storage11_7T/fuy/TADA-A/cell_WES/DNM/2021-01-19cp6_all_gene_PO_6788SNV_sf_uni_prior_compact.rds")
data=compact_data$base_info

gene_prior_file = "/storage11_7T/fuy/TADA-A/db/MS_data/prior/new_uniform_gene_prior.txt"
optimization_iteration = 2000

library(Rcpp)

sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/ll_sum.cpp")

sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/multi_annota.cpp")

library(data.table)
gene_prior = fread("/storage11_7T/fuy/TADA-A/db/MS_data/prior/new_uniform_gene_prior.txt")
gene_prior = gene_prior[order(gene_prior$genename),]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
fr <- function(x){
all_rr = x

logP_Zg0 = sumall0(data)
logP_Zg1 = sumall1(data,selected_annotations,all_rr)

logP_table<-data.table(logP_Zg1 = logP_Zg1, logP_Zg0 = logP_Zg0, genename = names(data))
logP_table <- logP_table[gene_prior, on = "genename"]
logP_table <- logP_table[complete.cases(logP_table)]
idx = match(unique(logP_table$genename),logP_table$genename)
idx2 = c(idx,nrow(logP_table) + 1)
pr = logP_table[idx,]$prior
ll_sum1 <- ll_sum(idx2,logP_table$logP_Zg1,logP_table$logP_Zg0,pr)
ll_sum1 
    }
1
2
selected_annotations = c(1)
fr(.1)

-36676.782281677

1
2
selected_annotations = seq(1,2)
fr(rep(.1,2))

-36675.829491401

1
2
selected_annotations = seq(1,5)
fr(rep(.1,5))

-36673.9502840269

1
2
3
4
5
6
7
8
selected_annotations = seq(1,2)
tmm = proc.time()

gama = optim(rep(0.1, 2), fr ,control=list("fnscale"=-1, "maxit" = optimization_iteration), hessian = TRUE)$par

proc.time() - tmm

gama
   user  system elapsed 
882.703   0.333 810.098
  1. 1.19776143563448
  2. 1.67173709266479
1
810/60

13.5

1
2
3
4
5
6
7
8
selected_annotations = seq(1,13)
tmm = proc.time()

gama = optim(rep(0.1,13), fr ,control=list("fnscale"=-1, "maxit" = optimization_iteration), hessian = TRUE)$par

proc.time() - tmm

gama
     user    system   elapsed 
16817.386     7.781 15419.179
  1. 0.859144612544029
  2. 0.610726150706316
  3. 0.110649587883844
  4. 2.06212251131498
  5. 0.221748031808197
  6. 1.3862515026335
  7. -2.01959850632312
  8. -0.314347856181316
  9. -1.0356541085996
  10. -1.47924971831386
  11. 1.13112632604613
  12. 0.362758870702958
  13. 0.899116918698962
1
15419/3600

4.28305555555556