Dev bowtie2

.agents/skills/monsda/reference/config-guide.md

@@ -109,6 +109,20 @@ Each step listed in `WORKFLOWS` needs:
 - `ENV`/`BIN` at the condition level override the step-level default (useful
   to benchmark two installs of the same tool, or point at a local binary for
   tools without a conda package, e.g. `guppy`, `dorado`).
+- `SEQUENCING` can likewise be set at the step level (or nested inside a
+  tool's own block) to override the global `SETTINGS.SEQUENCING` for just
+  that step, same `"single"/"paired"/"singlecell"[,rf|fr]` syntax. Useful for
+  e.g. QuantSeq-style 3' protocols: keep the global value `"paired,rf"` so
+  TRIMMING/DEDUP see both reads (read2 carries the UMI), but set
+  `"SEQUENCING": "single,rf"` under `MAPPING` (and `COUNTING` if needed) so
+  only read1 is actually mapped/counted. Steps downstream of MAPPING (PEAKS,
+  TRACKS, COUNTING, CIRCS, FUSIONS, DE, DEU, DAS) default to the resolved
+  MAPPING pairedness/strandedness, not the global one, when they don't set
+  their own — a MAPPING-level override already propagates to DE/DEU/DAS
+  automatically, e.g. so their featureCounts calls use the right `-p`/`-s`
+  flags. DTU is the exception: its default salmon path quantifies from the
+  trimmed FASTQ directly, independent of MAPPING, so it falls back to the
+  global value instead.
 - Inside each tool's block, `OPTIONS` holds one key per sub-step (e.g.
   `INDEX` vs `MAP` for a mapper) with a **raw command-line string** of extra
   flags for that tool — exactly as you'd type them, minus anything MONSDA

.agents/skills/monsda/reference/workflows.md

@@ -10,7 +10,7 @@ can put under that step's `TOOLS` key. Env names match a `.yaml` in `envs/`.
 | `BASECALL` | FAST5/POD5 -> FASTQ | `guppy: $PATH_TO_LOCAL_BIN` (no conda pkg), `dorado: $PATH_TO_LOCAL_BIN` (no conda pkg) |
 | `QC` | QC of FASTQ/BAM (runs standalone pre-processing, or again post-processing for any later step) | `fastqc: fastqc` (FASTQ/BAM, includes MultiQC), `rustqc: rustqc` (BAM only, includes MultiQC, has a strandedness guard — aborts the run if inferred strandedness conflicts with `SEQUENCING`) |
 | `TRIMMING` | Adapter/quality trimming | `trimgalore: trim_galore`, `cutadapt: cutadapt`, `bbduk: bbduk`, `fastp: fastp` |
-| `MAPPING` | Align to reference | `hisat2: hisat2`, `star: STAR` (also used for `STARsolo` single-cell mode), `rustar: rustar` (Rust STAR reimplementation), `segemehl2`/`segemehl3`/`segemehl2bisulfite`/`segemehl3bisulfite: segemehl.x`, `bwa: bwa mem`, `bwa2: bwa2-mem`, `bwameth: bwameth.py`, `minimap: minimap2` (long-read, use `-x map-ont`/`-x map-pb` in OPTIONS), `salmonalign: salmon` (alignment mode), `rammap: rammap` (long-read, not on bioconda), `piscem: piscem` (single-cell only, emits RAD for `alevinfry`) |
+| `MAPPING` | Align to reference | `hisat2: hisat2`, `star: STAR` (also used for `STARsolo` single-cell mode), `rustar: rustar` (Rust STAR reimplementation), `segemehl2`/`segemehl3`/`segemehl2bisulfite`/`segemehl3bisulfite: segemehl.x`, `bwa: bwa mem`, `bwa2: bwa2-mem`, `bwameth: bwameth.py`, `bowtie2: bowtie2` (no native stranded-protocol flag, so the strandedness part of `SEQUENCING` is not forwarded to it), `minimap: minimap2` (long-read, use `-x map-ont`/`-x map-pb` in OPTIONS), `salmonalign: salmon` (alignment mode), `rammap: rammap` (long-read, not on bioconda), `piscem: piscem` (single-cell only, emits RAD for `alevinfry`) |
 | `DEDUP` | Remove/mark duplicates | `umitools: umi_tools`, `fgumi: fgumi` (UMI extraction + dedup), `picard: picard` (position-based) |
 | `COUNTING` | Quantify reads/transcripts | `countreads: featureCounts`, `salmon: salmon` (FASTQ or, via `OPTIONS.BAM = "transcriptome"\|"genome"`, directly from BAM), `kallisto: kallisto`, `oarfish: oarfish` (long-read, from transcriptome BAM), `alevinfry: alevin-fry` (needs a piscem RAD file), `simpleaf: simpleaf` (self-contained single-cell index+map+quant) |
 | `TRACKS` | trackdb.txt + BIGWIG for UCSC/genome browsers | `ucsc: ucsc` |

AGENTS.md

@@ -14,7 +14,7 @@ State clearly what was actually verified and what was not.
 `snakemake` and `nextflow` live in the `monsda` conda env, they are *not* on the
 default PATH:
 
-    export PATH=/home/fall/anaconda3/envs/monsda/bin:$PATH
+    export PATH=${HOME}/anaconda3/envs/monsda/bin:$PATH
 
 Confirmed there: snakemake 9.20.0, nextflow 26.04.0.12031. `multiqc` is **not**
 installed in that env, MultiQC behaviour cannot be verified locally.

MONSDA/Params.py

@@ -1628,6 +1628,36 @@ def checkstranded(sample: list, config: dict) -> str:
     return stranded
 
 
+def resolve_sequencing(seqval: str, default_paired: str, default_stranded: str):
+    """Resolve a per-workflow-step SEQUENCING override (e.g. 'single,rf') into
+    a (paired, stranded) tuple, falling back to the given defaults (usually
+    the global SETTINGS.SEQUENCING derived values, or the resolved MAPPING
+    values for steps downstream of MAPPING) when the step itself does not
+    define its own override. This allows e.g. a quantseq MAPPING step to run
+    in 'single' mode while TRIMMING/DEDUP still see the global 'paired' mode
+    needed for UMI extraction from read2.
+
+    Parameters
+    ----------
+    seqval : str
+        SEQUENCING value configured for the step, e.g. 'single,rf', or falsy
+        if not configured for this step.
+    default_paired : str
+        paired value to fall back to when seqval is not set.
+    default_stranded : str
+        stranded value to fall back to when seqval is not set.
+
+    Returns
+    -------
+    tuple(str, str)
+        (paired, stranded)
+    """
+    if seqval and seqval != "":
+        parts = str(seqval).split(",")
+        return parts[0], (parts[1] if len(parts) > 1 else "")
+    return default_paired, default_stranded
+
+
 @check_run
 def set_pairing(samples: list, config: dict) -> list:
     """get pairs of samples

MONSDA/Workflows.py

@@ -652,11 +652,25 @@ def create_subworkflow(config, subwork, conditions, envs=None, stage=None):
                             tempconf["RUNDEDUP"] = "enabled"
                         continue
                     if mu.sub_dict(config[subwork], condition).get("ANNOTATION"):
-                        tempconf[subwork]["ANNOTATION"] = mu.sub_dict(
+                        subwork_annotation = mu.sub_dict(
                             config[subwork], condition
                         )["ANNOTATION"]
                     elif config[subwork].get("ANNOTATION"):
-                        tempconf[subwork]["ANNOTATION"] = config[subwork]["ANNOTATION"]
+                        subwork_annotation = config[subwork]["ANNOTATION"]
+                    else:
+                        subwork_annotation = None
+                    # per-step SEQUENCING override (e.g. 'single,rf'), same
+                    # condition-level lookup as ANNOTATION above; a tool-nested
+                    # override (alongside a tool's OPTIONS) is already carried
+                    # through via the per-tool subset_dict merge below.
+                    if mu.sub_dict(config[subwork], condition).get("SEQUENCING"):
+                        subwork_sequencing = mu.sub_dict(
+                            config[subwork], condition
+                        )["SEQUENCING"]
+                    elif config[subwork].get("SEQUENCING"):
+                        subwork_sequencing = config[subwork]["SEQUENCING"]
+                    else:
+                        subwork_sequencing = None
                     if ("TOOLS" in config[key] and env == "" and exe == "") and len(
                         mu.get_from_dict(config[key], condition)
                     ) != 0:  # env and exe overrule TOOLS
@@ -674,6 +688,14 @@ def create_subworkflow(config, subwork, conditions, envs=None, stage=None):
                                 tempconf[key] = mu.subset_dict(config[key], condition)
                     else:
                         tempconf[key] = mu.subset_dict(config[key], condition)
+                    # (Re-)apply the ANNOTATION/SEQUENCING overrides: the TOOLS
+                    # handling above can replace tempconf[key] wholesale (e.g.
+                    # when no envs filter is given, as for DE/DEU/DAS/DTU),
+                    # which would otherwise silently drop them again.
+                    if subwork_annotation:
+                        tempconf[subwork]["ANNOTATION"] = subwork_annotation
+                    if subwork_sequencing:
+                        tempconf[subwork]["SEQUENCING"] = subwork_sequencing
 
                     if stage and stage == "POST":
                         tempconf["SAMPLES"] = mp.get_samples_postprocess(config, key)
@@ -1968,6 +1990,9 @@ def nf_fetch_params(
         retconf["PAIRED"] = paired
     if stranded:
         retconf["STRANDED"] = stranded
+    # Per-step pairedness/strandedness, default to the global values above and
+    # overridden below wherever a workflow step defines its own SEQUENCING key.
+    MAPPING_PAIRED, MAPPING_STRANDED = paired, stranded
     if rundedup:
         retconf["RUNDEDUP"] = rundedup
     if prededup:
@@ -1986,6 +2011,10 @@ def nf_fetch_params(
         REF = mapper_dict.get("REFERENCE", MAPCONF.get("REFERENCE"))
         MANNO = mapper_dict.get("ANNOTATION", MAPCONF.get("ANNOTATION"))
         MDECOY = mapper_dict.get("DECOY", MAPCONF.get("DECOY"))
+        MSEQ = mapper_dict.get("SEQUENCING", MAPCONF.get("SEQUENCING"))
+        MAPPING_PAIRED, MAPPING_STRANDED = mp.resolve_sequencing(
+            MSEQ, paired, stranded
+        )
         if REF:
             REFERENCE = REF
             REFDIR = str(os.path.dirname(REFERENCE))
@@ -2070,6 +2099,8 @@ def nf_fetch_params(
         retconf["MAPPINGUIDXNAME"] = UIDXNAME
         retconf["MAPPINGUIDX2NAME"] = UIDX2NAME
         retconf["MAPPINGPREFIX"] = PREFIX
+        retconf["MAPPINGPAIRED"] = MAPPING_PAIRED
+        retconf["MAPPINGSTRANDED"] = MAPPING_STRANDED
 
     # Peak Calling Variables
     if "PEAKS" in config:
@@ -2085,6 +2116,13 @@ def nf_fetch_params(
             ANNOTATION = (
                 ANNO.get("GTF") if "GTF" in ANNO else ANNO.get("GFF")
             )  # by default GTF forma
+        # downstream steps consume the MAPPING output, so default to its pairedness
+        PSEQ = PEAKCONF.get("SEQUENCING")
+        PEAKS_PAIRED, PEAKS_STRANDED = mp.resolve_sequencing(
+            PSEQ, MAPPING_PAIRED, MAPPING_STRANDED
+        )
+        retconf["PEAKPAIRED"] = PEAKS_PAIRED
+        retconf["PEAKSTRANDED"] = PEAKS_STRANDED
         if not IP:
             IP = mp.check_IP(SAMPLES, config)
         log.debug(logid + "Running Peak finding for " + IP + " protocol")
@@ -2133,9 +2171,29 @@ def nf_fetch_params(
             XCONF = mu.sub_dict(config[x], SETUP)
             XENV = XENV.split("_")[0]
             log.debug(logid + "XCONFIG: " + str(SETUP) + "\t" + str(XCONF))
-            REF = XCONF[XENV].get("REFERENCE", XCONF.get("REFERENCE"))
-            XANNO = XCONF[XENV].get("ANNOTATION", XCONF.get("ANNOTATION"))
-            XDECOY = XCONF[XENV].get("DECOY", XCONF.get("DECOY"))
+            # Subconfigs from the non-combo nf_make_sub path have the toolenv
+            # key stripped out (same as MAPPING above); fall back to XCONF
+            # itself when the toolenv key is absent, otherwise REFERENCE/
+            # ANNOTATION/DECOY/SEQUENCING would silently raise KeyError/be lost.
+            x_dict = XCONF.get(XENV) if mu.dict_inst(XCONF) else None
+            if not mu.dict_inst(x_dict):
+                x_dict = XCONF
+            REF = x_dict.get(
+                "REFERENCE", XCONF.get("REFERENCE", config[x].get("REFERENCE"))
+            )
+            XANNO = x_dict.get(
+                "ANNOTATION", XCONF.get("ANNOTATION", config[x].get("ANNOTATION"))
+            )
+            XDECOY = x_dict.get("DECOY", XCONF.get("DECOY", config[x].get("DECOY")))
+            XSEQ = x_dict.get(
+                "SEQUENCING", XCONF.get("SEQUENCING", config[x].get("SEQUENCING"))
+            )
+            # TRACKS/COUNTING consume the MAPPING output, so default to its pairedness
+            x_paired, x_stranded = mp.resolve_sequencing(
+                XSEQ, MAPPING_PAIRED, MAPPING_STRANDED
+            )
+            retconf[x + "PAIRED"] = x_paired
+            retconf[x + "STRANDED"] = x_stranded
             if XANNO:
                 ANNOTATION = XANNO
             else:
@@ -2194,9 +2252,36 @@ def nf_fetch_params(
             XBIN, XENV = mp.env_bin_from_config(config, x)
             XENV = XENV.split("_")[0]
             log.debug(logid + "XCONFIG: " + str(SETUP) + "\t" + str(XCONF))
-            REF = XCONF[XENV].get("REFERENCE", XCONF.get("REFERENCE"))
-            XANNO = XCONF[XENV].get("ANNOTATION", XCONF.get("ANNOTATION"))
-            XDECOY = XCONF[XENV].get("DECOY", XCONF.get("DECOY"))
+            # Subconfigs from the non-combo nf_make_sub path have the toolenv
+            # key stripped out (same as MAPPING/TRACKS/COUNTING above); fall
+            # back to XCONF itself when the toolenv key is absent, otherwise
+            # REFERENCE/ANNOTATION/DECOY/SEQUENCING would silently raise
+            # KeyError/be lost.
+            x_dict = XCONF.get(XENV) if mu.dict_inst(XCONF) else None
+            if not mu.dict_inst(x_dict):
+                x_dict = XCONF
+            REF = x_dict.get(
+                "REFERENCE", XCONF.get("REFERENCE", config[x].get("REFERENCE"))
+            )
+            XANNO = x_dict.get(
+                "ANNOTATION", XCONF.get("ANNOTATION", config[x].get("ANNOTATION"))
+            )
+            XDECOY = x_dict.get("DECOY", XCONF.get("DECOY", config[x].get("DECOY")))
+            XSEQ = x_dict.get(
+                "SEQUENCING", XCONF.get("SEQUENCING", config[x].get("SEQUENCING"))
+            )
+            if x == "DTU":
+                # DTU's default salmon path quantifies directly from
+                # TRIMMED_FASTQ, independent of the MAPPING step, so fall
+                # back to the global layout
+                x_paired, x_stranded = mp.resolve_sequencing(XSEQ, paired, stranded)
+            else:
+                # DE/DEU/DAS consume the MAPPING output, so default to its pairedness
+                x_paired, x_stranded = mp.resolve_sequencing(
+                    XSEQ, MAPPING_PAIRED, MAPPING_STRANDED
+                )
+            retconf[x + "PAIRED"] = x_paired
+            retconf[x + "STRANDED"] = x_stranded
             if XANNO:
                 ANNOTATION = XANNO
             else:
@@ -2319,9 +2404,16 @@ def nf_fetch_params(
                 ANNO.get("GTF") if "GTF" in ANNO else ANNO.get("GFF")
             )  # by default GTF format will be used
 
+        CSEQ = CIRCCONF.get("SEQUENCING")
+        CIRCS_PAIRED, CIRCS_STRANDED = mp.resolve_sequencing(
+            CSEQ, MAPPING_PAIRED, MAPPING_STRANDED
+        )
+
         retconf["CIRCSREF"] = REFERENCE
         retconf["CIRCSREFDIR"] = REFDIR
         retconf["CIRCSANNO"] = ANNOTATION
+        retconf["CIRCSPAIRED"] = CIRCS_PAIRED
+        retconf["CIRCSSTRANDED"] = CIRCS_STRANDED
 
     # FUSIONS Variables
     if "FUSIONS" in config:
@@ -2340,9 +2432,16 @@ def nf_fetch_params(
                 ANNO.get("GTF") if "GTF" in ANNO else ANNO.get("GFF")
             )  # by default GTF format will be used
 
+        FSEQ = FUSCONF.get("SEQUENCING")
+        FUSIONS_PAIRED, FUSIONS_STRANDED = mp.resolve_sequencing(
+            FSEQ, MAPPING_PAIRED, MAPPING_STRANDED
+        )
+
         retconf["FUSIONSREF"] = REFERENCE
         retconf["FUSIONSREFDIR"] = REFDIR
         retconf["FUSIONSANNO"] = ANNOTATION
+        retconf["FUSIONSPAIRED"] = FUSIONS_PAIRED
+        retconf["FUSIONSSTRANDED"] = FUSIONS_STRANDED
         retconf["FUSIONSLIB"] = FUSCONF.get(FUSENV, {}).get("OPTIONS", {}).get(
             "INDEX", ""
         ) or os.path.join(REFDIR, "CTAT", FUSENV)

configs/template.json

@@ -163,6 +163,7 @@
                      "hisat2": "hisat2",
                      "bwa": "bwa mem",
                      "bwameth": "bwameth",
+                     "bowtie2": "bowtie2",
                      "minimap": "minimap2"
                  },
                  "id": {
@@ -219,6 +220,15 @@
                                     "MAP": "",
                                     "EXTENSION" : ""
                                 }
+                              },
+                              "bowtie2": {
+                                "OPTIONS":
+                                {
+                                    "INDEX" : "",
+                                    "MAP": "",
+                                    "EXTENSION" : ""
+                                },
+                                "SEQUENCING": "single,rf"  # per-step override, e.g. map only read1 for quantseq-style UMI protocols while TRIMMING/DEDUP keep the global paired SEQUENCING from SETTINGS
                               }
                          }
                      }