In [2]:
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")
sourceCpp(file="/storage11_7T/fuy/TADA-A/cell_WES/DNM/ll_sum2.cpp")
In [3]:
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
In [ ]:
library(data.table)
gene_prior = fread(gene_prior_file)
gene_prior = gene_prior[order(gene_prior$genename),]
selected_annotations = c(1,2,3)
logP_Zg0 = sumall0(data)
In [5]:
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
}
ll.lst = c()
gm.lst = c()
In [32]:
for(i in seq(0.16,0.2,0.01)){
a = c(rep(.1,3),i)
tmm = proc.time()
gama = optim(a, fr_pi ,control=list("fnscale"=-1, "maxit" = optimization_iteration))$par
gm.lst = rbind(gm.lst,gama)
ll = fr_pi(gama)
ll.lst = c(ll.lst,ll)
proc.time() - tmm
gm.lst
ll.lst
}
In [33]:
gm.lst
ll.lst
In [34]:
which.max(ll.lst)
In [35]:
plot(ll.lst)
