@@ -1919,7 +1919,8 @@ def compute_fold_change_with_combined_p_values_plot_heatmaps(combine_p_values_me
1919
1919
combine_p_values_method ,
1920
1920
num_of_real_data_avg_overlap ,
1921
1921
nucleosome_file ,
1922
- epigenomics_dna_elements )
1922
+ epigenomics_dna_elements ,
1923
+ log_file )
1923
1924
1924
1925
if plot_detailed_epigemomics_heatmaps :
1925
1926
# Plot heatmaps using step2 under epigenomics_occupancy/heatmaps/detailed
@@ -1942,7 +1943,8 @@ def compute_fold_change_with_combined_p_values_plot_heatmaps(combine_p_values_me
1942
1943
combine_p_values_method ,
1943
1944
num_of_real_data_avg_overlap ,
1944
1945
nucleosome_file ,
1945
- epigenomics_dna_elements )
1946
+ epigenomics_dna_elements ,
1947
+ log_file )
1946
1948
1947
1949
step3_sbs_signature2dna_element2avg_fold_change_dict , \
1948
1950
step3_dbs_signature2dna_element2avg_fold_change_dict , \
@@ -2104,16 +2106,16 @@ def calculate_fold_change_real_over_sim(center,
2104
2106
epigenomics_biosamples ,
2105
2107
occupancy_type ):
2106
2108
2107
- avg_real_signal = None
2108
- avg_sim_signal = None
2109
+ avg_real_signal = None
2110
+ avg_sim_signal = None
2109
2111
fold_change = None
2110
- min_sim_signal = None
2111
- max_sim_signal = None
2112
- pvalue = None
2113
- num_of_sims = None
2114
- num_of_sims_with_not_nan_avgs = None
2115
- real_data_avg_count = None
2116
- sim_avg_count = None
2112
+ min_sim_signal = None
2113
+ max_sim_signal = None
2114
+ pvalue = None
2115
+ num_of_sims = None
2116
+ num_of_sims_with_not_nan_avgs = None
2117
+ real_data_avg_count = None
2118
+ sim_avg_count = None
2117
2119
simulationsHorizontalMeans = None
2118
2120
2119
2121
biosample = None
@@ -2138,7 +2140,7 @@ def calculate_fold_change_real_over_sim(center,
2138
2140
else :
2139
2141
real_data_avg_signal_array = readData (None , signature , SIGNATUREBASED , output_dir , jobname , occupancy_type , dna_element_to_be_read ,AVERAGE_SIGNAL_ARRAY )
2140
2142
2141
- if real_data_avg_signal_array is not None :
2143
+ if real_data_avg_signal_array is not None and ( len ( real_data_avg_signal_array ) > 0 ) and ( not np . isnan ( real_data_avg_signal_array [ start : end ]). all ()) :
2142
2144
#If there is nan in the list np.mean returns nan.
2143
2145
avg_real_signal = np .nanmean (real_data_avg_signal_array [start :end ])
2144
2146
@@ -2148,7 +2150,7 @@ def calculate_fold_change_real_over_sim(center,
2148
2150
else :
2149
2151
real_data_accumulated_count_array = readData (None , signature , SIGNATUREBASED , output_dir , jobname , occupancy_type , dna_element_to_be_read ,ACCUMULATED_COUNT_ARRAY )
2150
2152
2151
- if real_data_accumulated_count_array is not None :
2153
+ if real_data_accumulated_count_array is not None and ( len ( real_data_accumulated_count_array ) > 0 ) and ( not np . isnan ( real_data_accumulated_count_array [ start : end ]). all ()) :
2152
2154
#If there is nan in the list np.mean returns nan.
2153
2155
real_data_avg_count = np .nanmean (real_data_accumulated_count_array [start :end ])
2154
2156
@@ -2223,7 +2225,9 @@ def calculate_fold_change_real_over_sim(center,
2223
2225
2224
2226
2225
2227
if (avg_real_signal is not None ) and (avg_sim_signal is not None ):
2226
- fold_change = avg_real_signal / avg_sim_signal
2228
+
2229
+ if avg_sim_signal > 0 :
2230
+ fold_change = avg_real_signal / avg_sim_signal
2227
2231
2228
2232
if (simulationsHorizontalMeans is not None ):
2229
2233
zstat , pvalue = calculate_pvalue_teststatistics (avg_real_signal , simulationsHorizontalMeans )
@@ -2345,7 +2349,6 @@ def step1_calculate_p_value(fold_change_window_size,
2345
2349
signature2Biosample2DNAElement2PValueDict = {}
2346
2350
plusorMinus = fold_change_window_size // 2
2347
2351
2348
-
2349
2352
def update_dictionary (complete_list ):
2350
2353
# complete_list:[jobname, signature, biosample, dna_element, avg_real_signal, avg_sim_signal, fold_change, min_sim_signal,
2351
2354
# max_sim_signal, pvalue, num_of_sims, num_of_sims_with_not_nan_avgs, real_data_avg_count, sim_avg_count,
@@ -2415,7 +2418,8 @@ def step2_combine_p_value(signature2Biosample2DNAElement2PValueDict,
2415
2418
combine_p_values_method ,
2416
2419
num_of_real_data_avg_overlap ,
2417
2420
nucleosome_file ,
2418
- epigenomics_dna_elements ):
2421
+ epigenomics_dna_elements ,
2422
+ log_file ):
2419
2423
2420
2424
# Fill and return this dictionary
2421
2425
signature2biosample2pooled_dna_element2combined_p_value_list_dict = {}
@@ -2498,7 +2502,9 @@ def step2_combine_p_value(signature2Biosample2DNAElement2PValueDict,
2498
2502
try :
2499
2503
test_statistic , combined_p_value = scipy .stats .combine_pvalues (p_values_array , method = combine_p_values_method , weights = None )
2500
2504
except FloatingPointError :
2501
- print ('signature:%s dna_element:%s fold_change_list:%s p_value_list:%s' % (signature ,dna_element ,fold_change_list ,p_value_list ))
2505
+ log_out = open (log_file , 'a' )
2506
+ print ('signature:%s dna_element:%s fold_change_list:%s p_value_list:%s' % (signature , dna_element , fold_change_list , p_value_list ), file = log_out )
2507
+ log_out .close ()
2502
2508
if len (p_value_list )> 0 :
2503
2509
combined_p_value = p_value_list [0 ]
2504
2510
@@ -2537,7 +2543,8 @@ def step3_combine_p_value(signature2Biosample2DNAElement2PValueDict,
2537
2543
combine_p_values_method ,
2538
2544
num_of_real_data_avg_overlap ,
2539
2545
nucleosome_file ,
2540
- epigenomics_dna_elements ):
2546
+ epigenomics_dna_elements ,
2547
+ log_file ):
2541
2548
2542
2549
# Fill and return this dictionary
2543
2550
signature2dna_element2combined_p_value_list_dict = {}
@@ -2613,7 +2620,9 @@ def step3_combine_p_value(signature2Biosample2DNAElement2PValueDict,
2613
2620
try :
2614
2621
test_statistic , combined_p_value = scipy .stats .combine_pvalues (p_values_array , method = combine_p_values_method , weights = None )
2615
2622
except FloatingPointError :
2616
- print ('signature:%s dna_element:%s fold_change_list:%s p_value_list:%s' % (signature ,dna_element ,fold_change_list ,p_value_list ))
2623
+ log_out = open (log_file , 'a' )
2624
+ print ('signature:%s dna_element:%s fold_change_list:%s p_value_list:%s' % (signature ,dna_element ,fold_change_list ,p_value_list ), file = log_out )
2625
+ log_out .close ()
2617
2626
if len (p_value_list )> 0 :
2618
2627
combined_p_value = p_value_list [0 ]
2619
2628
0 commit comments