Add --ancientness_threshold flag

Author: merszym

conf/process.config

@@ -128,7 +128,7 @@ process {
     //
 
     withName:"BAM_DEAM_BEST" {
-        ext.args = "only_stats ${params.doublestranded ? 'doublestranded' : ''}"
+        ext.args = "${params.ancientness_threshold} only_stats ${params.doublestranded ? 'doublestranded' : ''}"
         publishDir = [
             [
                 path: "${outdir}/stats",
@@ -173,7 +173,7 @@ process {
                 pattern: "confidence.tsv",
             ],
         ]
-       ext.args = "${params.doublestranded ? 'doublestranded' : ''}"
+       ext.args = "${params.ancientness_threshold} ${params.doublestranded ? 'doublestranded' : ''}"
     }
 
     withName:"PLOT_DEAM" {

modules/local/bam_deam_stats/resources/usr/bin/bam_deam_stats.py

@@ -15,7 +15,7 @@ def binomial_ci(x, n, alpha=0.05):
     return f"{round(lower*100,1):.1f},{round(upper*100,1):.1f}"
 
 
-def main(bamfile, positions=3, stats_only=False, doublestranded=False):
+def main(bamfile, threshold=threshold, positions=3, stats_only=False, doublestranded=False):
     #store all reference bases
     all_first = []
     all_last = []
@@ -155,9 +155,9 @@ def main(bamfile, positions=3, stats_only=False, doublestranded=False):
     try:
         test51,test31 = float(p_deam51_95ci.split(',')[0]),float(p_deam31_95ci.split(',')[0])
         if test51 and test31:
-            if test51 > 9.5 or test31 > 9.5:
+            if test51 > threshold or test31 > threshold:
                 ancientness = '+'
-            if test51 > 9.5 and test31 > 9.5:
+            if test51 > threshold and test31 > threshold:
                 ancientness = '++'
     except ValueError:
         pass
@@ -190,10 +190,11 @@ def pprint(var):
 
 if __name__ == "__main__":
     bamfile = sys.argv[1]
+    threshold = float(sys.argv[2])
     only_stats = 'only_stats' in sys.argv
     doublestranded = 'doublestranded' in sys.argv
 
     # only_stats - dont write deaminated reads to file
     # doublestranded - use G>A rates at 5' end instead of C>T as in singlestranded library prep
 
-    main(bamfile, stats_only=only_stats, doublestranded=doublestranded)
+    main(bamfile, threshold=threshold, stats_only=only_stats, doublestranded=doublestranded)

nextflow.config

@@ -62,6 +62,7 @@ params {
     mapbwa_quality_cutoff     = 25    // filter out mapped reads with qualityscore < 25
     krakenuniq_min_kmers      = 129   // filter out families with < 129 kmers assined
     krakenuniq_min_reads      = 3     // filter out families with < 3 reads assigned
+    ancientness_threshold     = 9.5   // set the threshold for reporting a family as Ancientness '++' or '+'
     reportfilter_percentage   = 0.5   // for the 'filtered_report.tsv', Specify FamPercentage threshold
     reportfilter_breadth      = 0.5   // for the 'filtered_report.tsv', Specify ProportionExpectedBreadth threshold