From 90146a34bae2beb6c6f0aa9415566af9256afaf9 Mon Sep 17 00:00:00 2001 From: maxibor Date: Thu, 25 Jul 2024 11:17:35 +0200 Subject: [PATCH] feat: limit rescaling to window size --- pydamage/main.py | 1 + pydamage/rescale.py | 11 +++++------ 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/pydamage/main.py b/pydamage/main.py index 9d86c88..0f9e6e0 100644 --- a/pydamage/main.py +++ b/pydamage/main.py @@ -163,6 +163,7 @@ def pydamage_analyze( bam=bam, threshold=threshold, alpha=0.05, + wlen=wlen, damage_dict=df.to_dict(), read_dict=read_dict, grouped=group, diff --git a/pydamage/rescale.py b/pydamage/rescale.py index 0bcc0ab..7febc57 100644 --- a/pydamage/rescale.py +++ b/pydamage/rescale.py @@ -62,7 +62,7 @@ def rescale_qual(read_qual, dmg_pmf, damage_bases, reverse): def rescale_bam( - bam, threshold, alpha, damage_dict, read_dict, grouped, outname, threads + bam, threshold, alpha, wlen, damage_dict, read_dict, grouped, outname, threads ): """Rescale quality scores in BAM file using damage model @@ -70,6 +70,7 @@ def rescale_bam( bam (str): Path to BAM file threshold (float): Predicted accuracy threshold alpha (float): Q-value threshold + wlen (int): Window size for damage model damage_dict (dict): Damage model parameters read_dict (dict): Dictionary of read names grouped (bool): Grouped analysis @@ -95,15 +96,13 @@ def rescale_bam( pydam_ref = ref dmg = damage_model() if pydam_ref in read_dict: - pass_filter = False if ( threshold - and threshold <= damage_dict["predicted_accuracy"][pydam_ref] + and + threshold <= damage_dict["predicted_accuracy"][pydam_ref] ) and (alpha and alpha >= damage_dict["qvalue"][pydam_ref]): - pass_filter = True - if pass_filter: dmg_pmf = dmg.fit( - x=np.arange(400), + x=np.arange(wlen), p=damage_dict["damage_model_p"][pydam_ref], pmin=damage_dict["damage_model_pmin"][pydam_ref], pmax=damage_dict["damage_model_pmax"][pydam_ref],