Add optional STAR alignment with isolated BAM and junction outputs

.github/workflows/ci.yml

@@ -23,11 +23,14 @@ jobs:
           python-version: ${{ matrix.python-version }}
 
       - name: Install package
-        run: pip install -e .
+        run: pip install -e .[test]
 
       - name: Compile package
         run: python -m compileall src/wulfrna
 
+      - name: Run pytest
+        run: pytest -q
+
       - name: CLI help
         run: |
           wulfrna --help

README.md

@@ -25,6 +25,7 @@ Required commands:
 - `cutadapt`
 - `multiqc`
 - `salmon` (default backend) or `kallisto` (optional backend, when selected)
+- `STAR` and `samtools` (optional, required only when `--aligner star` is selected)
 
 Useful standard utilities expected on Linux nodes:
 - `gzip`, `zcat`, `awk`, `sed`, `grep`, `sort`, `head`, `tail`
@@ -64,6 +65,7 @@ This installs the `wulfrna` command.
 - shared: `combined_tx2gene.tsv`
 - Salmon backend: `salmon_index/`
 - kallisto backend: `kallisto_index/combined_transcripts.kidx`
+- STAR aligner (when `--aligner star`): `star_index/` containing non-empty `Genome`, `SA`, `SAindex`, `genomeParameters.txt`, `chrName.txt`, `chrLength.txt`, and `chrNameLength.txt`
 
 ## 4) Run command
 
@@ -79,13 +81,14 @@ Arguments:
   - kallisto mapping: `none -> (unstranded default)`, `forward -> --fr-stranded`, `reverse -> --rf-stranded`
 - `--threads N` (required): total threads (`N >= 1`)
 - `--quantifier {salmon|kallisto}` (optional, default `salmon`)
+- `--aligner {none|star}` (optional, default `none`): run STAR alignment from trimmed FASTQs and write isolated alignment outputs under `align/star/<sample>/`; Salmon/kallisto quantification and abundance matrices are still produced normally
 - `--single-end`, `--SE` (optional): switch layout from `paired_end` to `single_end`
 - `--fragment-length FLOAT` and `--fragment-sd FLOAT` (single-end kallisto only; both required and must be `> 0`)
 - `--min-mapping-rate FLOAT` (optional, default `0.90`): minimum acceptable tx2gene transcript mapping rate per sample (`0.0-1.0`)
 - `--dry-run` (optional): validate tools/reference/inputs, write metadata/status, then exit
 - `--genome NAME` (optional): resolve references as `<reference>/<NAME>/...`
 - `--no-resume` (optional): disable phase-level resume and rerun all phases
-- `--force-from {fastqc_raw,cutadapt,fastqc_trimmed,quant,aggregate,multiqc}` (optional): force rerun from the selected phase onward
+- `--force-from {fastqc_raw,cutadapt,fastqc_trimmed,align,quant,aggregate,multiqc}` (optional): force rerun from the selected phase onward
 
 Backward-compatibility note: legacy invocation without explicit `run` is still accepted, but `wulfrna run ...` is the intended interface.
 
@@ -96,10 +99,11 @@ pip install -e .
 wulfrna --help
 wulfrna run <workdir> --reference <reference_dir> --stranded reverse --threads 4 --dry-run
 wulfrna run <workdir> --reference <reference_dir> --stranded reverse --threads 4 --quantifier kallisto --dry-run
+wulfrna run <workdir> --reference <reference_dir> --stranded reverse --threads 4 --aligner star --dry-run
 ```
 
 Notes:
-- For dry-run to pass, required binaries must be in `PATH`, references must be complete for the selected quantifier, and `<workdir>/fastq` must contain valid FASTQ inputs for the selected layout (`paired_end` or `single_end`).
+- For dry-run to pass, required binaries must be in `PATH`, references must be complete for the selected quantifier and optional aligner, and `<workdir>/fastq` must contain valid FASTQ inputs for the selected layout (`paired_end` or `single_end`).
 - Single-end examples:
   - Salmon: `wulfrna run WORKDIR --reference REFDIR --stranded reverse --threads 4 --single-end --dry-run`
   - kallisto: `wulfrna run WORKDIR --reference REFDIR --stranded reverse --threads 4 --single-end --quantifier kallisto --fragment-length 200 --fragment-sd 20 --dry-run`
@@ -112,6 +116,15 @@ Expected primary outputs on full success:
 - `multiqc/multiqc_report.html`
 - `logs/tx2gene_mapping_stats.tsv`
 
+Additional outputs when `--aligner star` is selected (per sample):
+- `align/star/<sample>/Aligned.sortedByCoord.out.bam`
+- `align/star/<sample>/Aligned.sortedByCoord.out.bam.bai`
+- `align/star/<sample>/SJ.out.tab`
+- `align/star/<sample>/ReadsPerGene.out.tab`
+- `align/star/<sample>/Log.final.out`
+
+STAR outputs coexist with Salmon/kallisto outputs. STAR gene counts are not aggregated into `abundance/gene_expected_counts.tsv` or `abundance/gene_tpm.tsv`.
+
 Matrix behavior note:
 - Gene-level matrices include only genes observed in the transcript quantification input (unobserved zero-only genes are not emitted).
 
@@ -131,6 +144,7 @@ Phases:
 - `fastqc_raw`
 - `cutadapt`
 - `fastqc_trimmed`
+- `align` (only when `--aligner star`)
 - `quant`
 - `aggregate`
 - `multiqc`
@@ -144,6 +158,7 @@ Conservative compatibility rules:
 - If sample IDs changed: resume is blocked with an error.
 - If input layout changed (`paired_end` vs `single_end`): automatic resume is blocked with an error.
 - If `quantifier`, `reference_dir`, or `stranded` changed: `quant` and downstream phases rerun.
+- If `aligner` or STAR index content changed: `align` and downstream phases rerun when STAR alignment is enabled.
 - If `combined_tx2gene.tsv` fingerprint changed: `aggregate` and `multiqc` rerun.
 - If only thread count changed: resume is allowed.
 

docs/spec.md

@@ -8,17 +8,18 @@ Paired-end FASTQ input is the default, and single-end FASTQ input is supported w
 Primary command:
 
 ```bash
-wulfrna run WORKDIR --reference REFDIR --stranded {none|forward|reverse} --threads N [--quantifier {salmon|kallisto}] [--single-end|--SE] [--fragment-length FLOAT] [--fragment-sd FLOAT] [--dry-run] [--genome NAME]
+wulfrna run WORKDIR --reference REFDIR --stranded {none|forward|reverse} --threads N [--quantifier {salmon|kallisto}] [--aligner {none|star}] [--single-end|--SE] [--fragment-length FLOAT] [--fragment-sd FLOAT] [--dry-run] [--genome NAME]
 ```
 
 - `WORKDIR` must contain `fastq/` using the selected input layout.
 - `--reference` points to a reference root.
 - `--genome` (optional) resolves references under `<reference>/<genome>/`.
 - `--quantifier` defaults to `salmon` and selects the transcript quantification backend.
+- `--aligner` defaults to `none`; `--aligner star` runs STAR from trimmed FASTQs and writes isolated outputs under `align/star/<sample>/` while preserving Salmon/kallisto quantification outputs.
 - `--single-end` / `--SE` enables single-end input; paired-end is the default.
 - `--fragment-length` and `--fragment-sd` are required only for single-end kallisto runs and must be positive.
 - `--no-resume` disables automatic phase-level resume.
-- `--force-from` forces rerun from one phase onward (`fastqc_raw`, `cutadapt`, `fastqc_trimmed`, `quant`, `aggregate`, `multiqc`).
+- `--force-from` forces rerun from one phase onward (`fastqc_raw`, `cutadapt`, `fastqc_trimmed`, `align`, `quant`, `aggregate`, `multiqc`).
 
 ## Required input layout
 `WORKDIR/fastq/` must contain one of these layouts:
@@ -38,6 +39,7 @@ All must be in `PATH`:
 - multiqc
 - salmon (if `--quantifier salmon`)
 - kallisto (if `--quantifier kallisto`)
+- STAR and samtools (if `--aligner star`)
 
 ## Required reference files (resolved reference directory)
 Shared:
@@ -49,6 +51,17 @@ Salmon backend:
 kallisto backend:
 - `kallisto_index/combined_transcripts.kidx`
 
+STAR aligner (`--aligner star`):
+- `star_index/Genome`
+- `star_index/SA`
+- `star_index/SAindex`
+- `star_index/genomeParameters.txt`
+- `star_index/chrName.txt`
+- `star_index/chrLength.txt`
+- `star_index/chrNameLength.txt`
+
+All required STAR index files must be non-empty.
+
 ## Pipeline behavior
 1. Validate tools, references, and FASTQ inputs for the selected layout.
 2. Record metadata (`samples.tsv`, run parameters, versions).
@@ -57,6 +70,7 @@ kallisto backend:
    - FastQC on raw FASTQ
    - Cutadapt trimming
    - FastQC on trimmed FASTQ
+   - STAR alignment per sample when `--aligner star` is selected
    - Transcript quantification per sample (`salmon quant` or `kallisto quant`)
    - Aggregate transcript-level estimates to gene-level using `combined_tx2gene.tsv`
    - MultiQC over the work directory
@@ -68,6 +82,7 @@ Completed phases write:
 - `status/steps/fastqc_raw.done`
 - `status/steps/cutadapt.done`
 - `status/steps/fastqc_trimmed.done`
+- `status/steps/align.done` (only when `--aligner star`)
 - `status/steps/quant.done`
 - `status/steps/aggregate.done`
 - `status/steps/multiqc.done`
@@ -81,13 +96,14 @@ Resume is enabled by default. A phase is skipped only if:
 
 Machine-readable manifest:
 - `status/manifest.json`
-- Includes `workdir`, `reference_dir`, `quantifier`, `stranded`, `layout`, `sample_ids`, reference files used, tx2gene fingerprint, `created_at`, `updated_at`.
+- Includes `workdir`, `reference_dir`, `quantifier`, `aligner`, `stranded`, `layout`, `sample_ids`, reference files used, tx2gene fingerprint, STAR index fingerprint when applicable, `created_at`, `updated_at`.
 
 Manifest compatibility rules:
 - sample set change => hard error (no silent resume),
 - missing `layout` in old manifests is treated as `paired_end`,
 - layout mismatch between existing and current manifest => hard error (prevents silent paired-end/single-end mixing),
 - quantifier/reference_dir/stranded change => rerun quant and downstream,
+- aligner or STAR index path/content change => rerun align and downstream when STAR alignment is enabled,
 - `combined_tx2gene.tsv` fingerprint change => reuse quant, rerun aggregate + multiqc,
 - threads-only changes are resumable.
 
@@ -104,6 +120,15 @@ Primary outputs:
 - `multiqc/multiqc_report.html`
 - `logs/tx2gene_mapping_stats.tsv`
 
+STAR outputs when `--aligner star`:
+- `align/star/<sample>/Aligned.sortedByCoord.out.bam`
+- `align/star/<sample>/Aligned.sortedByCoord.out.bam.bai`
+- `align/star/<sample>/SJ.out.tab`
+- `align/star/<sample>/ReadsPerGene.out.tab`
+- `align/star/<sample>/Log.final.out`
+
+STAR outputs coexist with Salmon/kallisto outputs and do not replace or overwrite the gene-level abundance matrices.
+
 Status markers (`WORKDIR/status/`):
 - `RUNNING` while active
 - `SUCCESS` on complete run success

pyproject.toml

@@ -16,6 +16,9 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Bio-Informatics",
 ]
 
+[project.optional-dependencies]
+test = ["pytest"]
+
 [project.scripts]
 wulfrna = "wulfrna.cli:main"
 

src/wulfrna/cli.py

@@ -25,14 +25,20 @@ def build_parser() -> argparse.ArgumentParser:
     )
     run_parser.add_argument("--threads", required=True, type=int, help="Total threads")
     run_parser.add_argument("--quantifier", choices=["salmon", "kallisto"], default="salmon", help="Transcript quantification backend")
+    run_parser.add_argument(
+        "--aligner",
+        choices=["none", "star"],
+        default="none",
+        help="Optional alignment backend; 'star' runs STAR alignment and writes isolated outputs under align/star/<sample>/",
+    )
     run_parser.add_argument("--single-end", "--SE", dest="single_end", action="store_true", help="Use single-end FASTQ input layout")
     run_parser.add_argument("--fragment-length", type=float, help="Single-end fragment length for kallisto")
     run_parser.add_argument("--fragment-sd", type=float, help="Single-end fragment length standard deviation for kallisto")
     run_parser.add_argument("--dry-run", action="store_true", help="Validate inputs, write metadata, and exit without running analysis tools")
     run_parser.add_argument("--no-resume", action="store_true", help="Disable automatic phase-level resume and rerun all phases")
     run_parser.add_argument(
         "--force-from",
-        choices=["fastqc_raw", "cutadapt", "fastqc_trimmed", "quant", "aggregate", "multiqc"],
+        choices=["fastqc_raw", "cutadapt", "fastqc_trimmed", "align", "quant", "aggregate", "multiqc"],
         help="Force rerun from the selected phase onward",
     )
 

src/wulfrna/io.py

@@ -4,22 +4,54 @@
 import socket
 import subprocess
 import shutil
+from importlib import metadata
 from pathlib import Path
+from typing import Optional
 
 
 def now_iso() -> str:
     return dt.datetime.now().astimezone().strftime("%Y-%m-%d %H:%M:%S %Z")
 
 
-def git_commit() -> str:
-    proc = subprocess.run(["git", "rev-parse", "--short", "HEAD"], capture_output=True, text=True)
-    return proc.stdout.strip() if proc.returncode == 0 else "unknown"
+def package_version() -> str:
+    try:
+        return metadata.version("wulfrna")
+    except metadata.PackageNotFoundError:
+        return "unknown"
+
+
+def find_repo_root(package_file: Optional[Path] = None) -> Optional[Path]:
+    start = Path(__file__ if package_file is None else package_file).resolve()
+    for parent in start.parents:
+        if (parent / ".git").exists():
+            return parent
+    return None
+
+
+def git_commit() -> Optional[str]:
+    repo_root = find_repo_root()
+    if repo_root is None:
+        return None
+    proc = subprocess.run(["git", "rev-parse", "--short", "HEAD"], cwd=repo_root, capture_output=True, text=True)
+    return proc.stdout.strip() if proc.returncode == 0 and proc.stdout.strip() else None
+
+
+def pipeline_version(package_ver: Optional[str] = None, commit: Optional[str] = None) -> str:
+    resolved_package_ver = package_version() if package_ver is None else package_ver
+    resolved_commit = git_commit() if commit is None else commit
+    if resolved_commit is not None:
+        return f"package:{resolved_package_ver} git:{resolved_commit}"
+    return f"package:{resolved_package_ver}"
 
 
 def capture_versions_file(out: Path) -> None:
     out.parent.mkdir(parents=True, exist_ok=True)
+    package_ver = package_version()
+    commit = git_commit()
     with out.open("w", encoding="utf-8") as f:
-        f.write(f"pipeline_version: git:{git_commit()}\n")
+        f.write(f"pipeline_version: {pipeline_version(package_ver, commit)}\n")
+        f.write(f"package_version: {package_ver}\n")
+        f.write(f"git_commit: {commit if commit is not None else 'not available'}\n")
         f.write(f"datetime: {now_iso()}\n")
         f.write(f"hostname: {socket.gethostname()}\n")
 
@@ -33,6 +65,8 @@ def capture_versions_file(out: Path) -> None:
             "salmon": ["salmon", "--version"],
             "kallisto": ["kallisto", "version"],
             "multiqc": ["multiqc", "--version"],
+            "STAR": ["STAR", "--version"],
+            "samtools": ["samtools", "--version"],
         }
         for name, cmd in version_cmds.items():
             if shutil.which(cmd[0]) is None: