merge from upstream

pydamage/cli.py

@@ -74,6 +74,12 @@ def cli(ctx, outdir):
 @click.option(
     "-f", "--force", is_flag=True, help="Force overwriting of results directory"
 )
+@click.option(
+    "-r",
+    "--count_reverse",
+    is_flag=True,
+    help="Also count GtoA mutations in reverse reads",
+)
 @click.option(
     "-g",
     "--group",

pydamage/damage.py

@@ -16,7 +16,6 @@ def __init__(self, reference, al_handle, wlen):
             al_handle(pysam.AlignmentFile)
             wlen (int): window length
 
-
         """
         self.alignments = al_handle.fetch(reference)
         self.reference = reference
@@ -174,7 +173,7 @@ def check_model_fit(model_dict, wlen, verbose):
     return model_dict
 
 
-def test_damage(ref, bam, mode, wlen, show_al, process, verbose):
+def test_damage(ref, bam, mode, wlen, show_al, count_reverse, process, verbose):
     """Prepare data and run LRtest to test for damage
 
     Args:
@@ -183,6 +182,7 @@ def test_damage(ref, bam, mode, wlen, show_al, process, verbose):
         mode (str): opening mode of alignment file
         wlen (int): window length
         show_al (bool): Show alignment representations
+        count_reverse(bool): count reverse alignments as well
         process (int): Number of process for parallelization
         verbose (bool): Run in verbose mode
     Returns:

pydamage/main.py

@@ -20,26 +20,28 @@ def pydamage_analyze(
     process=1,
     outdir="",
     plot=False,
+    count_reverse=False,
     verbose=False,
     force=False,
     group=False,
 ):
 
     if group:
         pydamage_analyze_group(
-            bam, wlen, show_al, process, outdir, plot, verbose, force
+            bam, wlen, show_al, count_reverse, process, outdir, plot, verbose, force
         )
 
     else:
         pydamage_analyze_multi(
-            bam, wlen, show_al, process, outdir, plot, verbose, force
+            bam, wlen, show_al, count_reverse, process, outdir, plot, verbose, force
         )
 
 
 def pydamage_analyze_multi(
     bam,
     wlen=30,
     show_al=False,
+    count_reverse=False,
     process=1,
     outdir="",
     plot=False,
@@ -62,6 +64,7 @@ def pydamage_analyze_multi(
     """
     if verbose:
         print(f"Pydamage version {__version__}\n")
+    print("force", force)
     utils.makedir(outdir, force=force)
 
     if not verbose:

pydamage/model_fit.py

@@ -69,10 +69,12 @@ def fit_models(
     # Likelihood for model A - Damage Model
     binom_p_damage = model_A.fit(c_sites, **optim_A)
     LA = binom.logpmf(k=binom_k, n=binom_n, p=binom_p_damage)
+    sigma_damage = np.sqrt(binom_n * binom_p_damage * 1 - binom_p_damage)
 
     # Likelihood for model B - Null Model
     binom_p_null = model_B.fit(c_sites, **optim_B)
     LB = binom.logpmf(k=binom_k, n=binom_n, p=binom_p_null)
+    sigma_null = np.sqrt(binom_n * binom_p_null * 1 - binom_p_null)
 
     # Difference of number of paramters between model A and model B
     pdiff = len(model_A.kwds) - len(model_B.kwds)