1
2
3
4
5
6
# 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
1
2
3
4
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")
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/post.cpp")
1
2
3
4
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),]
selected_annotations = c(1)
1
2
3
4
5
a = c(0.1,.05)
tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gama
proc.time() - tmm
  1. 0.255723407377646
  2. 0.999829134189715
   user  system elapsed 
338.867   0.072 297.475 
1
297/60

4.95

1
2
3
4
5
a = c(exp(1),.05)
tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gama
proc.time() - tmm
  1. 1.89123499319297
  2. 0.0504755152875449
   user  system elapsed 
283.513   0.127 248.781 
1
2
3
4
5
a = c(exp(1),.9)
tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gama
proc.time() - tmm
  1. 1.89459951460619
  2. 0.050258497374786
   user  system elapsed 
291.112   0.080 257.158 
1
2
3
4
5
6
selected_annotations = c(1)
a = c(.1,.05)
tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gama
proc.time() - tmm
  1. 0.255723407377646
  2. 0.999829134189715
   user  system elapsed 
354.290   0.116 313.294 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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),]

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
3
4
5
6
7
8
selected_annotations = seq(1,2)
tmm = proc.time()

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

proc.time() - tmm

gama
   user  system elapsed 
449.122   0.208 412.240
  1. 1.19516468584344
  2. 1.6678512493307
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
logP_Zg0 = sumall0(data)

fr_pi <- function(par){

    all_rr = par[1:length(selected_annotations)] 
    gene_prior$prior <- rep(par[length(par)],nrow(gene_prior))
    
    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
3
4
5
6
7
8
selected_annotations = seq(1,2)
tmm = proc.time()

gama = optim(c(exp(1),exp(1),.05), fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par

proc.time() - tmm

gama
   user  system elapsed 
540.017   0.233 478.671
  1. 1.09919237592268
  2. 1.49613778568841
  3. 0.076776311502852
1
2
3
4
5
6
7
8
selected_annotations = seq(1,2)
tmm = proc.time()

gama = optim(c(.1,.1,.05), fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par

proc.time() - tmm

gama
   user  system elapsed 
940.195   0.363 832.371
  1. 1.10186122803905
  2. 1.5013327085674
  3. 0.0764190730631155
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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),]
gene_prior$prior = rep(0.0764190730631155,nrow(gene_prior))

post_fr <- function(x){
    all_rr = x
    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)]
    u = unique(logP_table$genename)
    idx = match(u,logP_table$genename)
    idx2 = c(idx,nrow(logP_table) + 1)
    pr = logP_table[idx,]$prior 
    post = post(idx2,logP_table$logP_Zg1,logP_table$logP_Zg0,pr)
    post.dt = data.table(genename =u,prior = post )   
    post.dt
}
1
max(mm$prior)

0.36704126967405

1
2
all_rr = c(1.19516468584344,1.6678512493307)
selected_annotations = c(1,2)
1
logP_Zg1 = sumall1(data,selected_annotations,all_rr)
1
2
3
4
5
6
7
8
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)]
    u = unique(logP_table$genename)
    idx = match(u,logP_table$genename)
    idx2 = c(idx,nrow(logP_table) + 1)
    pr = rep(0.05,length(idx))
    post = post(idx2,logP_table$logP_Zg1,logP_table$logP_Zg0,pr)