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] "Made a tmp folder for tmp files."
[1] "Finished calculating mutation count within per genomic window. Time consumed: 1.955s"
[1] "Finished fitting mutation rate model and calculating scaling factors. Time consumed: 3.976s"
[1] "Removed temporary files."
A matrix: 3 × 4 of type dbl
EstimateStd. Errorz valuePr(>|z|)
(Intercept) 0.0018433270.01345905 0.13695820.89106386
GC_content 0.0034250120.01363202 0.25124750.80162276
div_score-0.0223681650.01276638-1.75211450.07975413
 
1
2
3
4
f = fread("/storage11_7T/fuy/TADA-A/cell_WES/cd_uniform_scaling_factors.txt")
colnames(f)[2] = "z"
f$scaling_factor = f$z * (1465/1692)
fwrite(f[,c(1,3)],"/storage11_7T/fuy/TADA-A/cell_WES/cd_uniform_scaling_factors_1465_d_1692.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
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
system.time(compact_data <- TADA_A_read_info_by_chunks(
mut_files = "/storage11_7T/fuy/TADA-A/cell_WES/DNM/6166SNV.cd.window.cd_mis_pli_syn.7131DNM.allele.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/cell_WES/cd_uniform_scaling_factors_1465_d_1692.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/all_lof_alt_A.05-995.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_C.05-995.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_G.05-995.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_T.05-995.bed") ,
    
c("/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_A.ls05.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_C.ls05.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_G.ls05.bed",
 "/storage11_7T/fuy/TADA-A/annotation/driverMAPS/all_lof_alt_T.ls05.bed") ,  
    
c("/storage11_7T/fuy/TADA-A/annotation/MPC_score/fordist_constraint_official_mpc_values_v2_MPC_gt2_altA.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/fordist_constraint_official_mpc_values_v2_MPC_gt2_altC.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/fordist_constraint_official_mpc_values_v2_MPC_gt2_altG.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/fordist_constraint_official_mpc_values_v2_MPC_gt2_altT.bed"),
    
c("/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC12_alt_A.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC12_alt_C.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC12_alt_G.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC12_alt_T.bed"),
    
c("/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC01_alt_A.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC01_alt_C.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC01_alt_G.bed",
"/storage11_7T/fuy/TADA-A/annotation/MPC_score/chr.MPC01_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(),"_SNV_6166_sf1465_d_1692_pli_mis_compact.rds"))
 user   system  elapsed 
2.876    0.983 2971.767

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:6){
    
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: 364.058.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
[1] "Read in gene prior file. Time consumed: 0.004s."
[1] "Finished optimization. Time consumed: 306.742.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
[1] "Read in gene prior file. Time consumed: 0.005s."
[1] "Finished optimization. Time consumed: 456.712.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
[1] "Read in gene prior file. Time consumed: 0.005s."
[1] "Finished optimization. Time consumed: 308.202.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
[1] "Read in gene prior file. Time consumed: 0.005s."
[1] "Finished optimization. Time consumed: 266.43.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
[1] "Read in gene prior file. Time consumed: 0.005s."
[1] "Finished optimization. Time consumed: 334.862.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"



     user    system   elapsed 
12966.046    15.455  2037.049
A data.frame: 6 × 3
logRRlower_boundupper_bound
<dbl><dbl><dbl>
3.87605892.80508754.9470303
1.42399730.36152832.4864663
0.54422680.30622120.7822325
2.30038072.06305102.5377103
1.39378141.11948561.6680773
1.18819461.04510791.3312813
 
1
2
3
exp(1.1194856)
exp(1.3937814)
exp(1.6680773)

3.0632780477128

4.03006054734171

5.30196390481291

1
exp(4.9470303)

140.756338534684

1
2
3
4
RR = TADA_A_RR_estimate(data = compact_data$base_info, selected_annotations = seq(1,6), 
#           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)
[1] "Read in gene prior file. Time consumed: 0.004s."
[1] "Finished optimization. Time consumed: 10344.042.s"
[1] "Got confidence intervals of RR estimates."
[1] "Finished RR estimation!"
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
options(scipen=200)
tm = proc.time()
BF = TADA_A_get_BFs(data = compact_data$base_info,
               selected_annotations = seq(1,6),
               rr = RR$rr_report$logRR , 
#                additional_BF_file = c("/storage11_7T/fuy/TADA-A/annotation/cd_BF/Example_gene_coding_BF.txt"), 
               additional_BF_file = c("/storage11_7T/fuy/TADA-A/annotation/cd_BF/uniform_cd_BFs.txt"),
               TADA_p0 = 0.94)

g_BF = BF$gene_BF_table
g_BF = g_BF[order(g_BF$FDR_all),]
g_BF2 = g_BF[g_BF$FDR_all <= 0.1,]
proc.time() - tm
nrow(g_BF2)

g_BF2
[1] "Read in additional BF file /storage11_7T/fuy/TADA-A/annotation/cd_BF/uniform_cd_BFs.txt."
[1] "Flagged genes that don't have any bases with any informative annotation."
[1] "Got genenames without bases that have informative rr features."
[1] "Got logBF_noncoding."
[1] "Added coding and non-coding logBF. Time consumed: 1.727s."
[1] "Finished!"



   user  system elapsed 
  8.728   0.000   7.941

21

A data.table: 21 × 9
genenamelogBF_noncodinglogBF_codinglogBF_allBF_noncodingBF_codingBF_allFDR_codingFDR_all
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
SCN2A 24.310465024.31046536132728519.62630136132728519.626300.940.0000000004335866
SLC6A1 15.884952015.884952 7920400.547441 7920400.547440.940.0000009892220579
CHD8 11.234650011.234650 75708.849621 75708.849620.940.0000696229093144
DEAF1 7.7490210 7.749021 2319.300361 2319.300360.940.0017296141750668
KCNQ3 7.5356520 7.535652 1873.665281 1873.665280.940.0030421258699277
KDM5B 7.0559360 7.055936 1159.721921 1159.721920.940.0047565924493450
FOXP1 6.0898020 6.089802 441.333861 441.333860.940.0089744373086258
DPYSL2 5.7372500 5.737250 310.210291 310.210290.940.0138620583244326
ADCY5 5.4090330 5.409033 223.415451 223.415450.940.0196027620954920
STXBP1 5.3880990 5.388099 218.787011 218.787010.940.0243246873250874
PPP2R5D 4.7819880 4.781988 119.341351 119.341350.940.0326626692498743
TRIP12 4.7184260 4.718426 111.991811 111.991810.940.0401677201875189
CHD2 4.6339300 4.633930 102.917781 102.917780.940.0472405112924787
SYNGAP1 4.5585710 4.558571 95.447021 95.447020.940.0539373843374813
SF3B1 4.5293770 4.529377 92.700821 92.700820.940.0599795462870336
TRAF7 4.4497200 4.449720 85.602991 85.602990.940.0658997292383351
SLC6A13 4.3299900 4.329990 75.943551 75.943550.940.0720829452627882
SRCAP 4.3029170 4.302917 73.915061 73.915060.940.0777942733959040
DYNC1H1 4.2093640 4.209364 67.313691 67.313690.940.0836366642041547
EIF3G 4.1233640 4.123364 61.766701 61.766700.940.0895710558994734
GNAI1 4.1057960 4.105796 60.691051 60.691050.940.0950759886127658
 
1
2