Add pairtools min_mapq option and release 0.2.1

README.rst

@@ -111,6 +111,9 @@ Changelog
 ========= ====================================================================
 Version   Description
 ========= ====================================================================
+0.2.1     Expose pairtools ``min_mapq`` in config/schema and pass ``--min-mapq``
+          to pairtools parse; switch bwa_split options to ``-5SP`` for Hi-C;
+          raise pairtools memory to 8G.
 0.2.0     Production release.
 0.1.0     Migration to modern sequana_pipetools framework (get_shell/get_run,
           schema validation, apptainer support, Python 3.10+).

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api"
 
 [tool.poetry]
 name = "sequana-hic"
-version = "0.2.0"
+version = "0.2.1"
 description = "Hi-C analysis, snakemake, sequana, container, reproducibility"
 authors = ["Sequana Team"]
 license = "BSD-3"

sequana_pipelines/hic/config.yaml

@@ -97,11 +97,16 @@ fastp:
         mem: 8G
 
 
+# 5any vs 5unique: takes 5'-most alignment regardless of uniqueness → repeat reads survive
 pairtools:
     walk_policy: "5unique"
+    # Minimum MAPQ to keep an alignment in pairtools parse (default pairtools value is 1).
+    # Set to 0 to keep multi-mapped reads (BWA assigns MAPQ 0 to multimappers), so they
+    # appear in the Hi-C contact map.
+    min_mapq: 1
     threads: 4
     resources:
-        mem: 4G
+        mem: 8G
 
 
 #############################################################################
@@ -134,7 +139,7 @@ bwa:
 bwa_split:
     nreads: 100000
     index_algorithm: is
-    options: -5SP -M
+    options: -5SP
     threads: 4
     tmp_directory: ./tmp
     resources:

sequana_pipelines/hic/hic.rules

@@ -422,6 +422,7 @@ rule pairtools:
     params:
         # walk_policy could be mask, 5any, 5unique, 3any|3unique|all
         walk_policy=config["pairtools"]["walk_policy"],
+        min_mapq=config["pairtools"].get("min_mapq", 1),
         mem=config["pairtools"]["resources"]["mem"],
     threads: config["pairtools"]["threads"]
     resources:
@@ -431,7 +432,7 @@ rule pairtools:
         # pairtools parse generate 1 large file.  pairtools sort as well. Better to pipe them to get
         # only one file. pairs number is identical of course. same content but sorted so statistics
         # would be the same
-        pairtools parse --nproc-in {threads} --nproc-out {threads} {input.bam} -c {input.chrom_size} --walks-policy {params.walk_policy} |  \
+        pairtools parse --nproc-in {threads} --nproc-out {threads} {input.bam} -c {input.chrom_size} --walks-policy {params.walk_policy} --min-mapq {params.min_mapq} |  \
             pairtools sort --nproc {threads} --memory {params.mem} -o {output.sort} >{log} 2>&1
         """
 

sequana_pipelines/hic/schema.yaml

@@ -163,6 +163,9 @@ mapping:
                 type: str
                 required: True
                 enum: ["mask", "5any", "5unique", "3any", "3unique", "all"]
+            "min_mapq":
+                type: int
+                required: False
             "threads":
                 type: int
                 required: True