Add print statements and update templates

Author: jacobimarsh

Bvalcalc/templates/AraTha_Cds_Params.py

@@ -11,7 +11,7 @@
 u = 6.95e-9 * x # Mutation rate (all types) per bp, per generation [4]
 g = r * 50 * (1-f) * x # Gene conversion initiation rate per bp, per generation [5]
 k = 553 # Gene conversion tract length (bp) [5]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.2 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [6]
 f1 = 0.14 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [6]
 f2 = 0.44 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [6]
@@ -21,7 +21,10 @@
 time_of_change = 1 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [2]
 ## Advanced DFE parameters 
 h = 0.5 + (f-0.5*f) # Dominance coefficient of selected alleles, NOTE: this is h_eff for h=0.5, replace BOTH 0.5's with your dominance coefficient [Naive value]
-mean, shape, proportion_synonymous = 500, 0.5, 0.3 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (2Ns), shape, proportion synonymous] (!Requires --gamma_dfe) [Naive value]## Literature cited
+mean, shape, proportion_synonymous = 500/(2*Nanc), 0.5, 0.3 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [Naive value]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
+## Literature cited
 # [1] Platt et al. 2010 doi: 10.1371/journal.pgen.1000843
 # [2] Durvasula et al. 2017 doi: 10.1073/pnas.1616736114
 # [3] Rowan et al. 2019 doi: 10.1534/genetics.119.302406

Bvalcalc/templates/AraTha_Phastcons_Params.py

@@ -11,7 +11,7 @@
 u = 6.95e-9 * x # Mutation rate (all types) per bp, per generation [4]
 g = r * 50 * (1-f) * x # Gene conversion initiation rate per bp, per generation [5]
 k = 553 # Gene conversion tract length (bp) [5]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.28 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [6*]
 f1 = 0.33 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [6*]
 f2 = 0.35 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [6*]
@@ -21,7 +21,10 @@
 time_of_change = 1 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [2]
 ## Advanced DFE parameters 
 h = 0.5 + (f-0.5*f) # Dominance coefficient of selected alleles, NOTE: this is h_eff for h=0.5, replace BOTH 0.5's with your dominance coefficient [Naive value]
-mean, shape, proportion_synonymous = 500, 0.5, 0.3 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (2Ns), shape, proportion synonymous] (!Requires --gamma_dfe) [Naive value]## Literature cited
+mean, shape, proportion_synonymous = 500 / (2*Nanc), 0.5, 0.3 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [Naive value]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
+## Literature cited
 # [1] Platt et al. 2010 doi: 10.1371/journal.pgen.1000843
 # [2] Durvasula et al. 2017 doi: 10.1073/pnas.1616736114
 # [3] Rowan et al. 2019 doi: 10.1534/genetics.119.302406

Bvalcalc/templates/DroMel_Cds_Params.py

@@ -10,7 +10,7 @@
 u = 6e-9 * x # Mutation rate (all types) per bp, per generation [1]
 g = 1e-8 * x # Gene conversion initiation rate per bp, per generation [3]
 k = 440 # Gene conversion tract length (bp) [3]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.25 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [4]
 f1 = 0.49 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [4]
 f2 = 0.04 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [4]
@@ -20,7 +20,9 @@
 time_of_change = 0.45 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [6]
 ## Advanced DFE parameters 
 h = 0.5 # Dominance coefficient of selected alleles [Naive value]
-mean, shape, proportion_synonymous = 811/(2*Nanc), 0.347, 0.3 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (s), shape, proportion synonymous] (!Requires --gamma_dfe) [7]
+mean, shape, proportion_synonymous = 811/(2*Nanc), 0.347, 0.3 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [7]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1] Keightley et al 2014  doi: 10.1534/genetics.113.158758
 # [2] Comeron et al 2012 doi: 10.1371/journal.pgen.1002905

Bvalcalc/templates/DroMel_Phastcons_Params.py

@@ -10,7 +10,7 @@
 u = 6e-9 * x # Mutation rate (all types) per bp, per generation [1]
 g = 1e-8 * x # Gene conversion initiation rate per bp, per generation [3]
 k = 440 # Gene conversion tract length (bp) [3]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.25 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [4]
 f1 = 0.49 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [4]
 f2 = 0.04 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [4]
@@ -20,7 +20,9 @@
 time_of_change = 0.45 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [6]
 ## Advanced DFE parameters 
 h = 0.5 # Dominance coefficient of selected alleles [Naive value]
-mean, shape, proportion_synonymous = 37.9/(2*Nanc), 0.371, 0.0 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (s), shape, proportion synonymous] (!Requires --gamma_dfe) [7]
+mean, shape, proportion_synonymous = 37.9/(2*Nanc), 0.371, 0.0 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [7]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1] Keightley et al 2014  doi: 10.1534/genetics.113.158758
 # [2] Comeron et al 2012 doi: 10.1371/journal.pgen.1002905

Bvalcalc/templates/DroMel_Utr_Params.py

@@ -10,7 +10,7 @@
 u = 6e-9 * x # Mutation rate (all types) per bp, per generation [1]
 g = 1e-8 * x # Gene conversion initiation rate per bp, per generation [3]
 k = 440 # Gene conversion tract length (bp) [3]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.25 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [4]
 f1 = 0.49 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [4]
 f2 = 0.04 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [4]
@@ -20,7 +20,9 @@
 time_of_change = 0.45 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [6]
 ## Advanced DFE parameters 
 h = 0.5 # Dominance coefficient of selected alleles [Naive value]
-mean, shape, proportion_synonymous = 93.3/(2*Nanc), 0.24, 0.0 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (s), shape, proportion synonymous] (!Requires --gamma_dfe) [7]
+mean, shape, proportion_synonymous = 93.3/(2*Nanc), 0.24, 0.0 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [7]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1] Keightley et al 2014  doi: 10.1534/genetics.113.158758
 # [2] Comeron et al 2012 doi: 10.1371/journal.pgen.1002905

Bvalcalc/templates/HomSap_Cds_Params.py

@@ -10,7 +10,7 @@
 u = 1.25e-8 * x # Mutation rate (all types) per bp, per generation [2]
 g = 3.27e-8 * x # Gene conversion initiation rate per bp, per generation [3]
 k = 113 # Gene conversion tract length (bp) [3]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.51 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [4]
 f1 = 0.14 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [4]
 f2 = 0.14 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [4]
@@ -20,7 +20,9 @@
 time_of_change = 0.81 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [1]
 ## Advanced DFE parameters 
 h = 0.5 # Dominance coefficient of selected alleles [Naive value]
-mean, shape = 500, 0.5 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (2Ns), shape] (!Requires --gamma_dfe) [Naive value]
+mean, shape, proportion_synonymous = 100 / (2*Nanc), 1, 0.3 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [Naive value]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1] Gutenkunst et al 2009 doi: 10.1371/journal.pgen.1000695
 # [2] Kong et al 2012 doi: 10.1038/nature11396

Bvalcalc/templates/HomSap_Phastcons_Params.py

@@ -10,7 +10,7 @@
 u = 1.25e-8 * x # Mutation rate (all types) per bp, per generation [2]
 g = 3.27e-8 * x # Gene conversion initiation rate per bp, per generation [3]
 k = 113 # Gene conversion tract length (bp) [3]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.59 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [5]
 f1 = 0.18 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [5]
 f2 = 0.18 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [5]
@@ -20,7 +20,9 @@
 time_of_change = 0.81 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [1]
 ## Advanced DFE parameters 
 h = 0.5 # Dominance coefficient of selected alleles [Naive value]
-mean, shape, proportion_synonymous = 16, 0.1165, 0 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (2Ns), shape] (!Requires --gamma_dfe) [5]
+mean, shape, proportion_synonymous = 16 / (2*Nanc), 0.1165, 0 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [5]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1] Gutenkunst et al 2009 doi: 10.1371/journal.pgen.1000695
 # [2] Kong et al 2012 doi: 10.1038/nature11396

Bvalcalc/templates/SelfingParams.py

@@ -11,7 +11,7 @@
 u = 1e-8 * x # Mutation rate (all types) per bp, per generation [Naive value]
 g = 1e-8 * (1-f) * x # Gene conversion initiation rate per bp, per generation [Naive value]
 k = 500 # Gene conversion tract length (bp) [Naive value]
-## DFE parameters for ALL sites in annotated regions (Sum must equal 1)
+## Basic DFE parameters for ALL sites in annotated regions (Sum must equal 1), can overwrite with other --*_dfe flags
 f0 = 0.25 # Proportion of effectively neutral mutations with 0 <= |2Ns| < 1 (Note that 2Ns<5 does not contribute to BGS) [Naive value]
 f1 = 0.25 # Proportion of weakly deleterious mutations with 1 <= |2Ns| < 10 [Naive value]
 f2 = 0.25 # Proportion of moderately deleterious mutations with 10 <= |2Ns| < 100 [Naive value]
@@ -21,7 +21,9 @@
 time_of_change = 1 * Nanc # Time in generations ago that effective population size went from Nanc to Ncur (!Requires --pop_change) [Naive value]
 ## Advanced DFE parameters 
 h = 0.5 + (f-0.5*f) # Dominance coefficient of selected alleles, NOTE: this is h_eff for h=0.5, replace BOTH 0.5's with your dominance coefficient [Naive value]
-mean, shape, proportion_synonymous = 100, 1, 0.3 # Gamma distribution of DFE to discretize and replace f0-f3 [mean (2Ns), shape, proportion synonymous] (!Requires --gamma_dfe) [Naive value]
+mean, shape, proportion_synonymous = 100 / (2*Nanc), 1, 0.3 # Gamma distribution of DFE to discretize into 9 bins  [mean (s), shape, strictly neutral proportion] (!Requires --gamma_dfe) [Naive value]
+s_breaks = 0, 1/(2*Nanc), 10/(2*Nanc), 100/(2*Nanc), 1 # Custom DFE parameter controlling the homozygous selection coefficient (s) breakpoints (!Requires --custom_dfe) [Naive value]
+bin_proportions = 0.25, 0.25, 0.25, 0.25 # Custom DFE parameter controlling the proportion of mutations between each bin by s_breaks, overwriting the f1-f3 values above (!Requires --custom_dfe) [Naive value]
 ## Literature cited
 # [1]
 # [2]

Bvalcalc/utils/dfe_helper.py

@@ -221,6 +221,7 @@ def customDFE_getf0(s_breaks, bin_props, s_cutoff):
     new_props : np.ndarray
         Bin proportions corresponding to new_breaks.
     """
+
     s_breaks = np.array(s_breaks, dtype=float)
     bin_props = np.array(bin_props, dtype=float)
 
@@ -253,6 +254,12 @@ def customDFE_getf0(s_breaks, bin_props, s_cutoff):
         else:
             new_props.append(prop)
             new_breaks.append(right)
+    
+    print(f"Custom DFE parameters active (overwriting f1-f3)")
+    print(f"s_edges, homozygous selection coefficient breakpoint for each bin =", np.array2string(np.array(new_breaks, dtype=float),precision=3))
+    print(f"s_cutoff, minimum selection coefficient for which BGS is calculated (below added to f0) = {s_cutoff:.6f}")
+    print(f"f0, effectively neutral proportion in selected region = {f0:.6f} which is not used in BGS calculations")
+    print(f"f_x, proportion in each remaining bin, excluding f0 =", np.array2string(np.array(new_props, dtype=float),formatter={'float_kind': lambda x: f'{x:.6f}'}))
 
     return ( # Return s clean arrays for exporting to calcB
         f0,
@@ -293,4 +300,6 @@ def legacyDFE_to_bins(params):
         t4
     ], dtype=float)
 
+    print(f"Basic DFE parameters active by default (see --constant_dfe, --gamma_dfe, --custom_dfe)")
+
     return params