1
2
3
4
5
setwd("/storage11_7T/fuy/TADA-A/TADA-A/lib")
library(data.table)
library(doParallel)
library(foreach)
library(data.table)
1
2
3
4
5
6
# f = fread("/storage11_7T/fuy/TADA-A/cell_WES/cd_mis_pli_syn.6954DNM_6430fam_scaling_factors.txt")
# colnames(f)[2] = "z"
# f$scaling_factor = rep(1,nrow(f))
# head(f,2)

# fwrite(f[,c(1,3)],"/storage11_7T/fuy/TADA-A/cell_WES/cd_uniform_scaling_factors.txt",sep="\t")
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
a = TADA_A_DNM_generator(window_file =  '/storage11_7T/fuy/TADA-A/db/MS_data/windows_partition/cd_windows_with_div_score.bed',
                                 mutrate_scaling_files = c("/storage11_7T/fuy/TADA-A/cell_WES/cd_uniform_scaling_factors.txt"),
                                 sample_sizes = 6430,
                                 gene_prior_file = "/storage11_7T/fuy/TADA-A/db/MS_data/prior/uniform_gene_prior.txt",
                                 nonAS_noncoding_annotations = NA,
                                AS_noncoding_annotations = list(
                              c("/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_A.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_C.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_G.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_T.bed")),
 
                                 report_proportion = 1,
                                 chunk_partition_num =2,
                                 node_n = 2,
                    mutrate_ref_files = c(
       "/storage11_7T/fuy/TADA-A/db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_A.mutrate.bw",
       "/storage11_7T/fuy/TADA-A/db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_C.mutrate.bw",
       "/storage11_7T/fuy/TADA-A/db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_G.mutrate.bw",
       "/storage11_7T/fuy/TADA-A/db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_T.mutrate.bw"),                                 
                                    rr = 0,
                                 output_allele_info_files = "3allele.bed",
                                 output_bed_files = "3.bed",
                                 output_risk_genes_file = "3temp.riskgenes",
                                 compact_mut_output = NA,
                                 MPI = 1)
[1] "echo \"Temp files cleaned and mutation generation finished finished!\""
1
2
library(tadaA)
setwd("/storage11_7T/fuy/TADA-A")
1
2
3
4
5
6
7
# TADA_A_adjust_mutrate(window_file = 'db/MS_data/windows_partition/cd_windows_with_div_score.bed',
#                       mut_file = '/storage11_7T/fuy/TADA-A/TADA-A/lib/2.bed',
#                       genes = "all",
#                       sample_size = 6430,
#                       mutrate_mode = "regular",
#                       scaling_file_name = "/storage11_7T/fuy/TADA-A/TADA-A/lib/2bed_sf.txt",
#                       scale_features = c("GC_content", "div_score"))
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
system.time(compact_data <- TADA_A_read_info_by_chunks(
mut_files = "/storage11_7T/fuy/TADA-A/TADA-A/lib/3allele.bed",
    
# c("db/MS_data/DNM/Jiang_cases_DNM_with_allele_info.txt",
# "db/MS_data/DNM/Kong_cases_DNM_with_allele_info.txt",
# "db/MS_data/DNM/Michaelson_cases_DNM_with_allele_info.txt",
# "db/MS_data/DNM/Yuen_NM2015_cases_DNM_with_allele_info.txt",
# "db/MS_data/DNM/Wu_cases_DNM_with_allele_info.txt"),    
    
window_file =  'db/MS_data/windows_partition/cd_windows_with_div_score.bed',
    
# mutrate_scaling_files = c("/storage11_7T/fuy/TADA-A/TADA-A/lib/3bed_sf.txt"),  
mutrate_scaling_files = c("/storage11_7T/fuy/TADA-A/cell_WES/cd_uniform_scaling_factors.txt"),
    
# sample_sizes = 4265,
sample_sizes = 6430,    
    
# mutrate_scaling_files = c("db/MS_data/results/Jiang_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
#                          "db/MS_data/results/Kong_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
#                          "db/MS_data/results/Michaelson_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
#                          "db/MS_data/results/Yuen_NM2015_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
#                          "db/MS_data/results/Wu_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt"),    
    
# sample_sizes = c(32,78,10,162,32),    
    
    
gene_prior_file = "db/MS_data/prior/uniform_gene_prior.txt",
    
nonAS_noncoding_annotations = NA, #c( 
#                                  "/storage11_7T/fuy/TADA-A/annotation/ccr/ccrs.allchrom.gt90.bed"),
                                #"/storage11_7T/fuy/TADA-A/annotation/ccr/chr_ccr_gt95_syn_rm.bed"),

AS_noncoding_annotations = list(
# c("/storage11_7T/fuy/TADA-A/annotation/pLI/pLI_nlt_995_gnomad.v2.1.1.all_lofs_snv_alt_A.bed",
#  "/storage11_7T/fuy/TADA-A/annotation/pLI/pLI_nlt_995_gnomad.v2.1.1.all_lofs_snv_alt_C.bed",
#  "/storage11_7T/fuy/TADA-A/annotation/pLI/pLI_nlt_995_gnomad.v2.1.1.all_lofs_snv_alt_G.bed",
#  "/storage11_7T/fuy/TADA-A/annotation/pLI/pLI_nlt_995_gnomad.v2.1.1.all_lofs_snv_alt_T.bed")
                              c("/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_A.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_C.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_G.bed",
                                "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/cd.window.all_cd_syn_alt_T.bed")
),                    

report_proportion = 18665/18665,
#chunk_partition_num =1,
chunk = 2,
node_n = 2,
mutrate_ref_files = c("db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_A.mutrate.bw",
                   "db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_C.mutrate.bw",
                   "db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_G.mutrate.bw",
                   "db/MS_data/mutrate/Example_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_T.mutrate.bw")
    
# mutrate_ref_files = c("/storage11_7T/data.for.yuwen/new_mutrate/window.hg19.genome.tri2.alt_A.uq.bedGraph.sort.bw" ,
#                       "/storage11_7T/data.for.yuwen/new_mutrate/window.hg19.genome.tri2.alt_C.uq.bedGraph.sort.bw" ,
#                       "/storage11_7T/data.for.yuwen/new_mutrate/window.hg19.genome.tri2.alt_G.uq.bedGraph.sort.bw" ,
#                       "/storage11_7T/data.for.yuwen/new_mutrate/window.hg19.genome.tri2.alt_T.uq.bedGraph.sort.bw" )
                      

))

saveRDS(compact_data,paste0("/storage11_7T/fuy/TADA-A/TADA-A/lib/",Sys.Date(),"_sim_rr0_6166_syn_compact.rds"))
 user   system  elapsed 
1.797    0.501 1324.798

un-calibrated

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
tm = proc.time()
dfn = c()
for(i in 1){
    
RR = TADA_A_RR_estimate(data = compact_data$base_info, selected_annotations = c(i), 
#           gene_prior_file = "db/MS_data/prior/Example_gene_prior.txt", 
        gene_prior_file = "db/MS_data/prior/uniform_gene_prior.txt",
                        optimization_iteration = 2000)
    
dfn = rbind(dfn,RR$rr_report)
    }
proc.time() - tm
# saveRDS(RR,paste0("/storage11_7T/fuy/TADA-A/annotation/results/",Sys.Date(),"_wes_deepsea_brain_RR.rds"))

dfn
[1] "Read in gene prior file. Time consumed: 0.004s."
[1] "Finished optimization. Time consumed: 326.646.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"



    user   system  elapsed 
2143.883    2.392  326.667
A data.frame: 1 × 3
logRRlower_boundupper_bound
<dbl><dbl><dbl>
0.2644308-0.43235260.9612143
 
1
exp(0.2644308)

1.30268927397531

1
saveRDS(dfn,"/storage11_7T/fuy/TADA-A/cell_WES/sim/syn_rr0_6166_rr.rds")
1
2
3
4
5
6
7
8
l = length(compact_data$base_info)
l

s = 0
for(i in 1:l){
    s = s + compact_data$base_info[[i]]$`1`$sum_mut_rate
}
s

67106

1691.73066556607

1
2
3
4
5
6
7
8
l = length(compact_data$base_info)
l

s = 0
for(i in 1:l){
    s = s + compact_data$base_info[[i]]$`1`$sum_mut_count
}
s

67106

1716

1
2
cd /storage11_7T/fuy/TADA-A/TADA-A/lib
wc -l 3allele.bed
6052 3allele.bed