{
  "skill": "alterlab-rnaseq-quant",
  "evals": [
    {
      "id": "salmon-decoy-quant",
      "prompt": "I have paired-end bulk RNA-seq FASTQs for 6 samples and a human transcriptome plus the genome FASTA. Can you set up salmon to quantify them properly? I want it done the accurate way.",
      "expected_output": "Invokes alterlab-rnaseq-quant on the salmon path. Builds a DECOY-AWARE gentrome index first (decoys.txt from the genome sequence names, gentrome = transcripts concatenated before genome) and runs `salmon index`, then quantifies each sample with `salmon quant -l A --gcBias` (selective alignment is the default in salmon 2.x; `--validateMappings` is accepted but ignored). Explicitly warns that a salmon 2.x index must be rebuilt because 2.0 cannot read C++/pufferfish indices, and checks the inferred library type in lib_format_counts.json. Stops at per-sample quant.sf and offers to summarize to a gene matrix.",
      "assertions": [
        {
          "type": "should_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "decoy"
        },
        {
          "type": "behavior",
          "value": "Builds a decoy-aware gentrome index and runs salmon quant with --gcBias -l A, and warns that salmon 2.x cannot read a C++/pufferfish index so it must be rebuilt."
        }
      ]
    },
    {
      "id": "tximport-to-deseq2-handoff",
      "prompt": "I already ran salmon and I have a folder of quant.sf files, one per sample. I need to turn these into a gene-level count matrix I can feed into DESeq2. Can you help me aggregate them?",
      "expected_output": "Invokes alterlab-rnaseq-quant in the import/aggregation mode: builds a tx2gene map from the annotation (scripts/build_tx2gene.py) and summarizes the per-sample quant.sf files to a gene-level count matrix via tximport semantics (countsFromAbundance lengthScaledTPM), e.g. scripts/import_quant.py. Produces gene_counts.tsv and then hands the differential-expression step off to alterlab-pydeseq2 rather than calling DEGs itself.",
      "assertions": [
        {
          "type": "should_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "tximport"
        },
        {
          "type": "behavior",
          "value": "Aggregates quant.sf to a gene-level matrix using a tx2gene map and a lengthScaledTPM import (tximport semantics), then routes the actual differential-expression analysis to alterlab-pydeseq2."
        }
      ]
    },
    {
      "id": "kallisto-kb-pseudoalign",
      "prompt": "I'd rather use kallisto than salmon for speed on my bulk RNA-seq. Can you show me how to build the index and run kb count on my paired-end reads?",
      "expected_output": "Invokes alterlab-rnaseq-quant on the kallisto / kb-python path: builds a kallisto index (or `kb ref`) from the transcriptome and runs `kb count -x bulk` (or `kallisto quant`) on the paired-end reads. Notes that kallisto uses a plain transcriptome index (no decoy/gentrome construction), points to lr-kallisto `--long` for long reads, and explains that abundance.tsv/.h5 feeds tximport identically to salmon before handoff to alterlab-pydeseq2.",
      "assertions": [
        {
          "type": "should_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "kb count"
        },
        {
          "type": "behavior",
          "value": "Drives the kallisto/kb-python pseudoalignment workflow (kb ref / kb count -x bulk or kallisto index/quant) and notes kallisto needs no decoy construction."
        }
      ]
    },
    {
      "id": "near-miss-pydeseq2-de",
      "prompt": "I already have a gene-level count matrix (genes x samples) and a sample sheet with treated vs control. Run the differential expression analysis and give me the significant DEGs with a volcano plot.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-pydeseq2. The user already has the count matrix and wants differential expression (size-factor normalization, dispersion, Wald tests, BH-FDR, volcano plot) — that is pydeseq2's job. alterlab-rnaseq-quant only produces the abundance/count matrix from FASTQs; it does not call DEGs.",
      "assertions": [
        {
          "type": "should_not_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "alterlab-pydeseq2"
        }
      ]
    },
    {
      "id": "near-miss-single-cell-alevin",
      "prompt": "I have 10x Chromium single-cell RNA-seq FASTQs with cell barcodes and UMIs. Can I quantify them with salmon alevin to get a cell-by-gene matrix?",
      "expected_output": "Does NOT invoke this skill's bulk path. Explains that `salmon alevin` was REMOVED in salmon 2.x and that single-cell quantification now uses the piscem + alevin-fry pipeline, then routes downstream single-cell analysis to alterlab-scanpy / alterlab-scvi-tools. This skill covers bulk RNA-seq only and must not force barcoded single-cell data through the bulk salmon/kallisto path.",
      "assertions": [
        {
          "type": "should_not_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "alevin-fry"
        }
      ]
    },
    {
      "id": "near-miss-sarek-variants",
      "prompt": "I have whole-exome sequencing FASTQs and I need to align them and call germline SNPs and indels to get a VCF for variant annotation.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-nf-core-sarek. The user wants FASTQ-to-VCF germline variant calling (alignment, BQSR, variant calling), not transcript abundance quantification. alterlab-rnaseq-quant is for RNA-seq expression quantification with salmon/kallisto and does not align to a genome for variant calling.",
      "assertions": [
        {
          "type": "should_not_trigger",
          "value": true
        },
        {
          "type": "output_contains",
          "value": "alterlab-nf-core-sarek"
        }
      ]
    }
  ]
}

# Decoy-aware salmon index (gentrome)

## Why decoys

Selective alignment scores how well a read maps to the transcriptome. Some reads
actually originate from intronic / intergenic genomic sequence and would
otherwise be force-assigned to a transcript. Including the **genome as a decoy**
lets salmon recognize "this read aligns better to the genome than to any
transcript" and decline the spurious transcript assignment. This is the
recommended construction for accurate bulk quantification.

## Construction

A **gentrome** is the transcripts FASTA concatenated with the genome FASTA. The
`decoys.txt` file lists the genome sequence (chromosome/scaffold) names so salmon
knows which gentrome entries are decoys.

```bash
# decoys = the genome's sequence names (strip '>' and any description after a space)
grep "^>" genome.fa | sed 's/^>//; s/ .*//' > decoys.txt

# gentrome: transcripts FIRST, then genome (order matters — decoys go last)
cat transcripts.fa genome.fa > gentrome.fa

salmon index -t gentrome.fa -d decoys.txt -i salmon_index -k 31 -p 8
```

### Gotchas

- **Order matters.** Transcripts must come before the genome in the gentrome;
  `decoys.txt` names must match the genome FASTA headers exactly.
- **Header trimming.** Many genome FASTAs have headers like
  `>1 dna:chromosome ...`; the `sed 's/ .*//'` keeps only `1`. salmon matches on
  the first whitespace-delimited token, so trim consistently.
- **Rebuild for salmon 2.x.** The Rust rewrite uses a new index format and rejects
  C++/pufferfish indices outright — rebuild with the binary you will quantify with
  (see tool_versions.md).
- Use a transcriptome FASTA and genome from the **same assembly/annotation
  release** (e.g. matching Ensembl/GENCODE versions). Fetch references with
  `alterlab-gget` if needed.

## Helper

`scripts/make_decoys.py` automates the two derived files and adds a guard: it
extracts genome names into `decoys.txt`, concatenates the gentrome in the correct
order, and warns if the transcript and genome FASTAs look swapped or empty. It
shells out nothing — pure stdlib file IO — so it runs in a bare `uv` env.

```bash
uv run python scripts/make_decoys.py \
  --transcripts transcripts.fa \
  --genome genome.fa \
  --out-gentrome gentrome.fa \
  --out-decoys decoys.txt
```

Then run `salmon index -t gentrome.fa -d decoys.txt -i salmon_index`.

# kallisto & kb-python (pseudoalignment path)

Targets standalone kallisto **v0.52.0** and kb-python (the `kb` CLI). See
../references/tool_versions.md for the version facts.

## Plain kallisto

```bash
# index
kallisto index -i kallisto_index.idx transcripts.fa

# paired-end quant
kallisto quant -i kallisto_index.idx -o quants_kallisto/sampleA \
  sampleA_R1.fastq.gz sampleA_R2.fastq.gz

# single-end (must supply fragment length mean/SD)
kallisto quant -i kallisto_index.idx -o quants_kallisto/sampleA \
  --single -l 200 -s 20 sampleA.fastq.gz
```

Outputs per sample: `abundance.tsv` (`target_id`, `length`, `eff_length`,
`est_counts`, `tpm`), `abundance.h5`, and `run_info.json`. Feed `abundance.h5`
(or `abundance.tsv`) to `tximport` with `type="kallisto"`.

## kb-python (`kb`) — wraps kallisto | bustools

```bash
# build references (index + t2g) from a GENOME FASTA + GTF
# positional order: <GENOME_FASTA> <GTF> (genome first, annotation second)
kb ref -i index.idx -g t2g.txt -f1 cdna.fa genome.fa annotation.gtf

# bulk quantification
kb count -i index.idx -g t2g.txt -x bulk -o quants_kb/sampleA \
  sampleA_R1.fastq.gz sampleA_R2.fastq.gz
```

- `kb ref` extracts the cDNA from the genome+GTF, builds the kallisto index, and
  writes the transcript-to-gene map (`t2g.txt`) in one step.
- `-f1 cdna.fa` is the **output** cDNA FASTA that `kb ref` writes (not an input
  transcriptome). The two trailing positionals are the genome FASTA **then** the
  GTF — genome first, in that order.
- `-x` is the technology string; `bulk` for bulk RNA-seq.
- kb-python emits tidy count outputs and a run log.

## Long reads (lr-kallisto)

kb-python exposes **lr-kallisto** through the `--long` flag and supports k-mer
sizes greater than 31, for long-read cDNA (ONT/PacBio). kb-python v0.29.1
introduced this option. Use it instead of the short-read defaults for long-read
libraries; consult `kb count --help` in your installed version for the exact
`--long` usage and the recommended `-k`.

## Notes

- kallisto/kb does **not** use the salmon gentrome/decoy construction — that is
  salmon-specific. A plain transcriptome index is correct here.
- The bundled kallisto inside a given kb-python build may lag the standalone
  kallisto release; check your install.
- Downstream, kallisto output aggregates to gene level via `tximport`
  identically to salmon — see ../references/tximport_handoff.md — then hand the
  gene matrix to `alterlab-pydeseq2`.

# salmon quant — flags, library type, QC

Targets salmon **2.x** (the Rust rewrite; selective alignment by default). See
tool_versions.md for the version facts (2.0 cannot read C++/pufferfish indices;
`salmon alevin` removed; `--validateMappings` accepted but ignored).

## Mapping-based quant command

```bash
salmon quant \
  -i salmon_index \
  -l A \
  -1 sampleA_R1.fastq.gz -2 sampleA_R2.fastq.gz \
  --gcBias \
  -p 8 \
  -o quants/sampleA
```

| Flag | Meaning |
|------|---------|
| `-i` | path to the (decoy-aware, freshly built) index |
| `-l A` | **auto-detect** library type / strandedness; verify the result afterwards |
| `-1` / `-2` | paired-end read files (use `-r` for single-end) |
| `--gcBias` | correct fragment-level GC bias (recommended for DE) |
| `--seqBias` | optional: correct 5'/3' sequence-specific bias |
| `--posBias` | optional: correct positional (5'/3' coverage) bias |
| `--sketch` | opt **out** of selective alignment into faster alignment-free pseudoalignment (2.x) |
| `--ignoreTxVersion` | with `-g/--geneMap`, match transcript IDs ignoring the trailing `.N` (2.x) |
| `--numBootstraps` / `--numGibbsSamples` | inferential replicates, written in the same format fishpond/swish expect |
| `-p` | threads |
| `-o` | per-sample output directory |

Selective alignment is the default in 2.x, so there is no flag that enables it.
`--validateMappings` parses and warns; `--mimicBT2`, `--mimicStrictBT2`,
`--minAssignedFrags`, `--numBiasSamples` and `--alternativeInitMode` are removed and
now error.

## Outputs (per sample)

- `quant.sf` — transcript-level table: `Name`, `Length`, `EffectiveLength`,
  `TPM`, `NumReads`. This is the file `tximport` reads.
- `lib_format_counts.json` — the **inferred library type** and compatible
  fragment counts. Always check this when using `-l A`.
- `logs/salmon_quant.log` and `aux_info/meta_info.json` — overall mapping rate
  and run metadata; record the mapping rate as a QC metric.

## Library type (`-l`)

`-l A` lets salmon infer strandedness from the data. The inferred code (e.g.
`ISR`, `ISF`, `IU`) appears in `lib_format_counts.json`. Only override with an
explicit code when you have a documented, trustworthy protocol; a wrong manual
`-l` silently biases counts. If the inferred type is inconsistent across samples
that should share a protocol, investigate before proceeding.

## QC checks before handoff

- **Mapping rate** (`aux_info/meta_info.json` → `percent_mapped`): unexpectedly
  low rates suggest a contaminating organism, wrong reference, or adapter/quality
  issues upstream.
- **Consistent inferred library type** across replicates.
- **Effective length** sanity: very short effective lengths flag fragment-length
  distribution problems for single-end data (set `--fldMean`/`--fldSD` if you
  must quantify SE without a distribution).

## Single-end reads

```bash
salmon quant -i salmon_index -l A -r sampleA.fastq.gz \
  --gcBias -p 8 -o quants/sampleA
```

For SE data salmon cannot empirically learn the fragment-length distribution;
provide `--fldMean` and `--fldSD` if known.

# Single-cell RNA-seq: `salmon alevin` is removed

## The change

As of salmon **2.x**, **`salmon alevin` has been removed** (upstream release
notes; see ../references/tool_versions.md). Single-cell / droplet RNA-seq
quantification is no longer a salmon subcommand. Any instruction that writes
`salmon alevin ...` is outdated and will fail.

## The replacement: piscem + alevin-fry

COMBINE-lab directs former alevin users to the **piscem + alevin-fry** pipeline,
which the release notes describe as offering improved memory efficiency and
throughput and being actively maintained:

- **piscem** — the indexing + mapping tool: https://github.com/COMBINE-lab/piscem
- **alevin-fry** — single-cell quantification (collation, barcode/UMI
  resolution, count-matrix generation): https://github.com/COMBINE-lab/alevin-fry

The typical flow is: build a piscem index of the (spliced+intron, if RNA
velocity is desired) reference, map the single-cell reads to produce a RAD file,
then run alevin-fry `generate-permit-list` → `collate` → `quant` to produce the
cell × gene count matrix.

## Routing

This skill (`alterlab-rnaseq-quant`) covers **bulk** RNA-seq only. For
single-cell:

1. Quantify with piscem + alevin-fry (their docs are the source of truth for the
   exact subcommand syntax — verify before running, as the CLI evolves).
2. Take the resulting count matrix downstream to:
   - `alterlab-scanpy` — single-cell analysis (QC, clustering, UMAP).
   - `alterlab-scvi-tools` — probabilistic single-cell models.
   - `alterlab-anndata` — to load/manipulate the matrix as an AnnData object.

Do not attempt to force single-cell barcoded data through the bulk salmon/kallisto
path; the barcode/UMI structure requires the alevin-fry tooling.

# Pinned tool versions & upstream facts

All version-specific claims in this skill trace to the upstream release pages
below. Re-verify against the linked release notes before changing a pin.

## salmon — 2.x (the Rust rewrite)

- Source: COMBINE-lab/salmon — https://github.com/COMBINE-lab/salmon, with the
  breaking changes catalogued in the repo's `MIGRATION.md`.
- Version in bioconda at review time (2026-09-23): **2.7.0** (uploaded 2026-08-30).
- **salmon 2.0 is a from-scratch Rust rewrite.** Same workflow
  (`salmon index` -> `salmon quant` -> `quant.sf`), same downstream output formats,
  single portable binary with no Boost/compiler dependency. The final C++ release is
  salmon **1.12.0**, kept on the upstream `cpp` branch and packaged as `salmon-cpp`.
- **Index break.** 2.0 uses a new index format and *cannot read C++ (pufferfish)
  indices*; the mismatch is detected and rejected with a clear error in both
  directions. Rebuild with the binary you quantify with.
- **Outputs are stable.** `quant.sf` is unchanged, and inferential replicates
  (`aux_info/bootstrap/...` from `--numBootstraps` and `--numGibbsSamples`) keep the
  C++ format, so tximport / tximeta / fishpond / swish work unmodified. The bias
  diagnostic dumps in `aux_info/` moved to a documented Rust format; no standard R
  package reads them.
- **Removed subcommand:** `salmon alevin`. It prints a redirect and exits. Single-cell
  moved to the **piscem + alevin-fry** ecosystem
  (https://github.com/COMBINE-lab/piscem, https://github.com/COMBINE-lab/alevin-fry).
- **Removed options (now error):** `--features` (index); `--mimicBT2`,
  `--mimicStrictBT2`, `--minAssignedFrags`, `--alternativeInitMode`,
  `--bootstrapReproject`, `--noGammaDraw`, `--numBiasSamples` (quant);
  `--auxTargetFile`, `--writeOrphanLinks` (quant -a).
- **Accepted but ignored (parse + warn):** `--validateMappings` (selective alignment is
  the default), `--eqclasses`, `--noFragLengthDist`, `--noSingleFragProb`,
  `--mismatchSeedSkip`, `--disableChainingHeuristic`, `--hitFilterPolicy`,
  `--maxRecoverReadOcc`, `--filterSize` (index), and several `quant -a` options.
- **New in 2.0:** `--sketch` (alignment-free pseudoalignment mode),
  `--sketchStrictOrphans`, `--allowDovetail` honored in sketch mode, and
  `--ignoreTxVersion` for `-g/--geneMap` matching. With `-g`, unmatched transcripts are
  still emitted as single-transcript genes (nothing is dropped), but 2.x warns once
  with a count and writes the names to `aux_info/genemap_unmatched_txps.json` instead
  of warning per transcript.

## kallisto — v0.52.0

- Source: pachterlab/kallisto releases — https://github.com/pachterlab/kallisto/releases
- Latest release at authoring time: **v0.52.0** (released 25 Feb). This release
  restores features (pseudobam, genomebam, fusion) that were missing after the
  index-structure rework; v0.51.1 and v0.51.0 precede it.

## kb-python (the `kb` CLI)

- Source: pachterlab/kb_python releases — https://github.com/pachterlab/kb_python/releases
- Latest release at authoring time: **v0.30.2** (released 19 May).
- **lr-kallisto / `--long`.** kb-python v0.29.1 release notes: *"Added
  lr-kallisto (--long) option, and enabling k>31"* and shipped kallisto/bustools
  binaries with and without long k-mer support. That release upgraded the bundled
  kallisto to **0.51.1** and bustools to **0.44.1**. (The bundled kallisto in a
  given kb-python build may lag the standalone kallisto release above; check
  `kb info` / your install for the exact bundled version.)

## nf-core/rnaseq — v3.26.0 (turnkey alternative)

- Source: nf-core/rnaseq releases — https://github.com/nf-core/rnaseq/releases
- Latest release at authoring time: **v3.26.0** ("Chromium Cuttlefish", released
  7 May).
- **Default route is `--aligner star_salmon`.** From the v3.26.0 usage docs:
  *"By default, the pipeline uses STAR (i.e. `--aligner star_salmon`) to map the
  raw FastQ reads to the reference genome, project the alignments onto the
  transcriptome and to perform the downstream BAM-level quantification with
  Salmon."* (https://nf-co.re/rnaseq/3.26.0/docs/usage)

## tximport / tximeta

- Bioconductor packages used for aggregating transcript-level estimates to the
  gene level (`tximport`) with full transcriptome provenance (`tximeta`,
  `linkedTxome`). They are the canonical R import layer for salmon/kallisto
  output feeding DESeq2. The Python helper in this skill reproduces the
  `countsFromAbundance = "lengthScaledTPM"` computation so the workflow can stay
  in `uv`; use the R packages when you need full `tximeta` metadata.

# tximport / tximeta handoff to pydeseq2

## What tximport does

salmon `quant.sf` and kallisto `abundance.tsv/.h5` are **transcript-level**.
DESeq2/PyDESeq2 work at the **gene level**. `tximport` summarizes transcript
estimates to genes using a transcript→gene map (`tx2gene`) and, importantly,
generates an **average-transcript-length offset** so the model accounts for
changes in transcript usage between samples.

The recommended import for DESeq2 uses
`countsFromAbundance = "lengthScaledTPM"` (or `"scaledTPM"`), which produces
length-scaled, library-size-scaled counts that correct for between-sample
differences in average transcript length, suitable to feed directly into the
differential-expression model. (GC/sequence-bias correction is a separate step
done earlier by salmon's `--gcBias`/`--seqBias`, not by tximport.)

## tx2gene map

A two-column table mapping `transcript_id` → `gene_id`. Build it from the same
GTF/GFF3 annotation used to make the transcriptome FASTA. `scripts/build_tx2gene.py`
parses `transcript_id` and `gene_id` attributes from a GTF and writes a TSV.

```bash
uv run python scripts/build_tx2gene.py annotation.gtf --out tx2gene.tsv
```

Caveat: transcript IDs in the `quant.sf`/`abundance.tsv` `Name`/`target_id`
column must match the IDs in `tx2gene` (watch for version suffixes like
`ENST00000456328.2` vs `ENST00000456328` — strip with `--strip-version` if your
annotation and reference disagree).

## Python aggregation helper

`scripts/import_quant.py` reproduces the tximport `lengthScaledTPM` computation
in pure stdlib (no pandas) so the workflow can stay in `uv`:

```bash
uv run python scripts/import_quant.py \
  --quants quants \
  --tx2gene tx2gene.tsv \
  --out-counts gene_counts.tsv \
  --out-tpm gene_tpm.tsv
```

It reads every `<sample>/quant.sf` (auto-detects kallisto `abundance.tsv` too)
and applies tximport's `makeCountsFromAbundance(..., "lengthScaledTPM")` at the
transcript level before summarizing to genes:

1. `newCounts[tx, s] = TPM[tx, s] * mean_over_samples(effLength[tx, :])`
2. rescale each sample column so its total equals that sample's original
   mapped-read count: `* (sum NumReads[:, s] / sum newCounts[:, s])`
3. sum the scaled transcript counts per gene and round to integers.

The result is length-corrected and library-size-scaled (per-sample column totals
equal the salmon/kallisto library size) — **not** a plain sum of raw `NumReads`.
The gene × TPM matrix is the per-gene sum of transcript TPMs. The only
differences from Bioconductor `tximport` are integer rounding and the lack of
`tximeta` provenance metadata; for exact parity or provenance use the R path
below.

## When to use the R path instead

Prefer the Bioconductor `tximport` + `tximeta` packages when you need:

- Full **`tximeta` provenance** — `linkedTxome` records the exact transcriptome
  checksum, source, and release, attaching it to the imported object.
- The officially validated `dtuScaledTPM` mode for differential transcript usage,
  or exact parity with a published DESeq2/tximport pipeline.

In R the flow is `tximport(files, type="salmon", tx2gene=..., countsFromAbundance="lengthScaledTPM")`
→ `DESeqDataSetFromTximport(...)`. The Python helper here is a convenience that
covers the common count-matrix case.

## Handoff

The resulting `gene_counts.tsv` (gene × sample integer matrix) plus a sample
metadata sheet (sample → condition) are exactly the inputs **`alterlab-pydeseq2`**
expects. This skill stops here; pydeseq2 owns normalization, dispersion, Wald
tests, FDR, and plots.

#!/usr/bin/env python3
"""build_tx2gene.py — make a transcript->gene (tx2gene) map from a GTF/GFF.

tximport needs a two-column transcript_id -> gene_id table to summarize
transcript-level salmon/kallisto estimates up to genes. This parses the GTF
attribute column for transcript_id and gene_id and writes a deduplicated TSV.

Pure stdlib (no pandas needed). Handles both GTF style
  gene_id "ENSG..."; transcript_id "ENST...";
and GFF3 style
  ID=transcript:ENST...;Parent=gene:ENSG...
attribute fields, preferring GTF; falls back to GFF3 keys when GTF keys are absent.

Usage:
  uv run python build_tx2gene.py annotation.gtf --out tx2gene.tsv
  uv run python build_tx2gene.py annotation.gtf --strip-version --out tx2gene.tsv

Output: a TSV with header `transcript_id\tgene_id`, one row per unique transcript.

Exit codes: 0 = wrote map; 2 = bad input / no transcript records found.
"""

from __future__ import annotations

import argparse
import re
import sys
from pathlib import Path

# GTF: key "value"; GFF3: key=value within a ;-separated attribute string.
_GTF_ATTR = re.compile(r'(\w+)\s+"([^"]*)"')
_GFF_ATTR = re.compile(r"(\w+)=([^;]+)")


def _strip_ver(x: str) -> str:
    """ENST00000456328.2 -> ENST00000456328 (drop a trailing .<digits>)."""
    return re.sub(r"\.\d+$", "", x)


def parse_attrs(attr_field: str) -> dict[str, str]:
    attrs = {k: v for k, v in _GTF_ATTR.findall(attr_field)}
    if "transcript_id" not in attrs and "gene_id" not in attrs:
        # try GFF3 form
        for k, v in _GFF_ATTR.findall(attr_field):
            attrs.setdefault(k, v)
    return attrs


def transcript_gene(attrs: dict[str, str]) -> tuple[str, str] | None:
    tx = attrs.get("transcript_id")
    gene = attrs.get("gene_id")
    # GFF3 fallbacks
    if tx is None:
        raw = attrs.get("ID", "")
        if raw.startswith("transcript:"):
            tx = raw.split(":", 1)[1]
    if gene is None:
        parent = attrs.get("Parent", "")
        if parent.startswith("gene:"):
            gene = parent.split(":", 1)[1]
    if tx and gene:
        return tx, gene
    return None


def main(argv: list[str] | None = None) -> int:
    p = argparse.ArgumentParser(description=__doc__.split("\n")[0])
    p.add_argument("gtf", type=Path, help="GTF/GFF annotation (matching the transcriptome)")
    p.add_argument("--out", required=True, type=Path)
    p.add_argument(
        "--strip-version",
        action="store_true",
        help="drop trailing .N version suffixes from IDs (match Ensembl reference vs annotation)",
    )
    p.add_argument("--force", action="store_true", help="overwrite existing output")
    args = p.parse_args(argv)

    if not args.gtf.is_file():
        print(f"ERROR: annotation not found: {args.gtf}", file=sys.stderr)
        return 2
    if args.out.exists() and not args.force:
        print(f"ERROR: {args.out} exists (use --force)", file=sys.stderr)
        return 2

    seen: dict[str, str] = {}
    with args.gtf.open("r") as fh:
        for line in fh:
            if not line or line.startswith("#"):
                continue
            cols = line.rstrip("\n").split("\t")
            if len(cols) < 9:
                continue
            feature = cols[2]
            # only feature rows that carry transcript_id (transcript/exon/etc.)
            attrs = parse_attrs(cols[8])
            tg = transcript_gene(attrs)
            if tg is None:
                continue
            tx, gene = tg
            if args.strip_version:
                tx, gene = _strip_ver(tx), _strip_ver(gene)
            seen.setdefault(tx, gene)  # first mapping wins; one row per transcript
            _ = feature  # feature kind not needed once transcript_id is present

    if not seen:
        print(
            "ERROR: no transcript_id/gene_id pairs found. Is this a valid GTF/GFF3 "
            "with transcript_id and gene_id attributes?",
            file=sys.stderr,
        )
        return 2

    with args.out.open("w") as out:
        out.write("transcript_id\tgene_id\n")
        for tx, gene in sorted(seen.items()):
            out.write(f"{tx}\t{gene}\n")

    print(f"Wrote {args.out}: {len(seen)} transcript->gene rows.", file=sys.stderr)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

#!/usr/bin/env python3
"""import_quant.py — aggregate salmon/kallisto transcript estimates to a gene matrix.

A Python convenience that reproduces tximport's gene-level summarization so the
RNA-seq quant workflow can stay in `uv` before handing off to alterlab-pydeseq2.
It reads every per-sample quant file under a directory, maps transcripts to genes
via a tx2gene TSV, and writes:

  * a gene x sample COUNT matrix (lengthScaledTPM, integer-rounded) for DESeq2
  * a gene x sample TPM matrix

Auto-detects salmon (`<sample>/quant.sf`) and kallisto (`<sample>/abundance.tsv`)
outputs. Pure stdlib — no pandas required, so it runs in a bare `uv run` env.

SEMANTICS — countsFromAbundance="lengthScaledTPM": this implements tximport's
`makeCountsFromAbundance(..., countsFromAbundance="lengthScaledTPM")` at the
transcript level, then sums to genes (matching how tximport applies the scaling
before gene summarization). For each transcript t and sample s:

    newCounts[t, s] = TPM[t, s] * mean_over_samples(effLength[t, :])

then each sample column is rescaled so its total equals that sample's original
mapped-read total:

    counts[t, s] = newCounts[t, s] * (sum_t NumReads[t, s] / sum_t newCounts[t, s])

Finally counts are summed per gene and rounded. This is length-corrected and
library-size-scaled — NOT a plain sum of raw NumReads. For full tximeta
provenance or exact parity with a published DESeq2/tximport pipeline, use the
Bioconductor tximport/tximeta packages (see ../references/tximport_handoff.md).

Usage:
  uv run python import_quant.py --quants quants --tx2gene tx2gene.tsv \\
      --out-counts gene_counts.tsv --out-tpm gene_tpm.tsv

Exit codes: 0 = wrote matrices; 2 = bad input / no quant files found.
"""

from __future__ import annotations

import argparse
import sys
from pathlib import Path

# salmon quant.sf cols:   Name  Length  EffectiveLength  TPM  NumReads
# kallisto abundance.tsv: target_id  length  eff_length  est_counts  tpm


def load_tx2gene(path: Path, strip_version: bool) -> dict[str, str]:
    mapping: dict[str, str] = {}
    with path.open("r") as fh:
        for i, line in enumerate(fh):
            parts = line.rstrip("\n").split("\t")
            if len(parts) < 2:
                continue
            tx, gene = parts[0], parts[1]
            if i == 0 and tx.lower() in {"transcript_id", "target_id", "name"}:
                continue  # header
            if strip_version:
                tx = tx.split(".")[0]
            mapping[tx] = gene
    return mapping


def find_quant_files(root: Path) -> list[tuple[str, Path, str]]:
    """Return (sample_name, file_path, kind) for each per-sample subdir."""
    found: list[tuple[str, Path, str]] = []
    for sub in sorted(p for p in root.iterdir() if p.is_dir()):
        salmon = sub / "quant.sf"
        kallisto = sub / "abundance.tsv"
        if salmon.is_file():
            found.append((sub.name, salmon, "salmon"))
        elif kallisto.is_file():
            found.append((sub.name, kallisto, "kallisto"))
    return found


def parse_quant(path: Path, kind: str, strip_version: bool):
    """Yield (transcript_id, eff_length, tpm, num_reads) rows."""
    with path.open("r") as fh:
        header = fh.readline().rstrip("\n").split("\t")
        idx = {name: i for i, name in enumerate(header)}
        if kind == "salmon":
            c_id, c_eff, c_tpm, c_reads = (
                idx["Name"], idx["EffectiveLength"], idx["TPM"], idx["NumReads"],
            )
        else:  # kallisto
            c_id, c_eff, c_tpm, c_reads = (
                idx["target_id"], idx["eff_length"], idx["tpm"], idx["est_counts"],
            )
        for line in fh:
            cols = line.rstrip("\n").split("\t")
            if len(cols) <= max(c_id, c_eff, c_tpm, c_reads):
                continue
            tx = cols[c_id].split(".")[0] if strip_version else cols[c_id]
            try:
                yield tx, float(cols[c_eff]), float(cols[c_tpm]), float(cols[c_reads])
            except ValueError:
                continue


def main(argv: list[str] | None = None) -> int:
    p = argparse.ArgumentParser(description=__doc__.split("\n")[0])
    p.add_argument("--quants", required=True, type=Path, help="dir of per-sample quant subdirs")
    p.add_argument("--tx2gene", required=True, type=Path)
    p.add_argument("--out-counts", required=True, type=Path)
    p.add_argument("--out-tpm", required=True, type=Path)
    p.add_argument("--strip-version", action="store_true")
    args = p.parse_args(argv)

    if not args.quants.is_dir():
        print(f"ERROR: --quants dir not found: {args.quants}", file=sys.stderr)
        return 2
    if not args.tx2gene.is_file():
        print(f"ERROR: --tx2gene not found: {args.tx2gene}", file=sys.stderr)
        return 2

    tx2gene = load_tx2gene(args.tx2gene, args.strip_version)
    if not tx2gene:
        print("ERROR: tx2gene map is empty.", file=sys.stderr)
        return 2

    samples = find_quant_files(args.quants)
    if not samples:
        print(
            f"ERROR: no quant.sf or abundance.tsv under {args.quants}/*/ — "
            "did salmon/kallisto write per-sample subdirs?",
            file=sys.stderr,
        )
        return 2

    sample_names: list[str] = [s[0] for s in samples]

    # First pass: read every transcript row from every sample. We need the full
    # per-transcript effective-length and TPM/NumReads tables BEFORE we can apply
    # lengthScaledTPM, because the length factor is averaged over samples and the
    # rescale factor depends on per-sample totals.
    #   tx_eff[tx][sample]   = effective length
    #   tx_tpm[tx][sample]   = TPM
    #   tx_reads[tx][sample] = NumReads (raw, used only to recover library size)
    tx_eff: dict[str, dict[str, float]] = {}
    tx_tpm: dict[str, dict[str, float]] = {}
    tx_reads: dict[str, dict[str, float]] = {}
    reads_total: dict[str, float] = {s: 0.0 for s in sample_names}
    unmapped = 0
    seen_tx: set[str] = set()

    for sample, path, kind in samples:
        for tx, eff, tpm, reads in parse_quant(path, kind, args.strip_version):
            reads_total[sample] += reads  # library size = sum of ALL mapped reads
            if tx2gene.get(tx) is None:
                if tx not in seen_tx:
                    unmapped += 1
                    seen_tx.add(tx)
                continue
            seen_tx.add(tx)
            tx_eff.setdefault(tx, {})[sample] = eff
            tx_tpm.setdefault(tx, {})[sample] = tpm
            tx_reads.setdefault(tx, {})[sample] = reads

    if not tx_tpm:
        print(
            "ERROR: no transcripts mapped to genes. Check that tx2gene IDs match "
            "the quant file IDs (try --strip-version).",
            file=sys.stderr,
        )
        return 2

    # lengthScaledTPM, transcript level (tximport::makeCountsFromAbundance):
    #   newCounts[tx, s] = TPM[tx, s] * mean_over_samples(effLength[tx, :])
    # then rescale each sample column so its total matches the sample's original
    # mapped-read count (countsSum / newSum). Effective lengths <= 0 are skipped
    # in the row mean (salmon/kallisto can emit 0 for unexpressed transcripts).
    avg_eff: dict[str, float] = {}
    for tx, per_sample in tx_eff.items():
        vals = [v for v in per_sample.values() if v > 0]
        avg_eff[tx] = (sum(vals) / len(vals)) if vals else 0.0

    new_counts: dict[str, dict[str, float]] = {}
    new_sum: dict[str, float] = {s: 0.0 for s in sample_names}
    for tx, tpm_row in tx_tpm.items():
        L = avg_eff.get(tx, 0.0)
        row = {}
        for s in sample_names:
            nc = tpm_row.get(s, 0.0) * L
            row[s] = nc
            new_sum[s] += nc
        new_counts[tx] = row

    # Per-sample rescale factor countsSum/newSum (countsSum = library size).
    scale: dict[str, float] = {}
    for s in sample_names:
        scale[s] = (reads_total[s] / new_sum[s]) if new_sum[s] > 0 else 0.0

    # Summarize scaled transcript counts (and raw TPM) up to gene level.
    gene_counts: dict[str, dict[str, float]] = {}
    gene_tpm: dict[str, dict[str, float]] = {}
    for tx, row in new_counts.items():
        gene = tx2gene[tx]
        gc = gene_counts.setdefault(gene, {})
        gt = gene_tpm.setdefault(gene, {})
        tpm_row = tx_tpm[tx]
        for s in sample_names:
            gc[s] = gc.get(s, 0.0) + row[s] * scale[s]
            gt[s] = gt.get(s, 0.0) + tpm_row.get(s, 0.0)

    genes = sorted(gene_counts)

    def write_matrix(out: Path, data: dict[str, dict[str, float]], integer: bool) -> None:
        with out.open("w") as fh:
            fh.write("gene_id\t" + "\t".join(sample_names) + "\n")
            for g in genes:
                row = data.get(g, {})
                vals = []
                for s in sample_names:
                    v = row.get(s, 0.0)
                    vals.append(str(int(round(v))) if integer else f"{v:.4f}")
                fh.write(g + "\t" + "\t".join(vals) + "\n")

    write_matrix(args.out_counts, gene_counts, integer=True)
    write_matrix(args.out_tpm, gene_tpm, integer=False)

    print(
        f"Wrote {args.out_counts} and {args.out_tpm}: "
        f"{len(genes)} genes x {len(sample_names)} samples "
        f"({unmapped} transcript rows had no gene mapping).",
        file=sys.stderr,
    )
    print("Next: hand gene_counts.tsv + a condition sheet to alterlab-pydeseq2.", file=sys.stderr)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

#!/usr/bin/env python3
"""make_decoys.py — build a decoy-aware salmon gentrome + decoys.txt.

A decoy-aware salmon index needs two derived files:
  * decoys.txt  — the genome's sequence (chromosome/scaffold) names, one per line
  * gentrome.fa — transcripts FASTA concatenated with genome FASTA (transcripts FIRST)

salmon then treats the genome entries as decoys, so reads that align better to the
genome than to any transcript are not force-assigned to a transcript.

This helper is pure stdlib (streaming file IO — never loads a whole FASTA into
memory) and adds guards against the two common silent mistakes:
  * empty / unreadable input FASTAs
  * a swapped argument order (genome passed as transcripts) — detected heuristically
    by comparing header counts (transcriptomes have far more records than genomes)

Usage:
  uv run python make_decoys.py --transcripts transcripts.fa --genome genome.fa \\
      --out-gentrome gentrome.fa --out-decoys decoys.txt

Then:
  salmon index -t gentrome.fa -d decoys.txt -i salmon_index -k 31 -p 8

Exit codes: 0 = wrote both files; 2 = bad input/usage.
"""

from __future__ import annotations

import argparse
import sys
from pathlib import Path


def _first_token(header_line: str) -> str:
    """'>1 dna:chromosome ...' -> '1'. Strip '>' and anything after first space."""
    name = header_line[1:].strip()
    # salmon matches on the first whitespace-delimited token
    return name.split()[0] if name else ""


def collect_decoys(genome_path: Path) -> list[str]:
    """Stream the genome FASTA and return its sequence names (for decoys.txt)."""
    names: list[str] = []
    with genome_path.open("r") as fh:
        for line in fh:
            if line.startswith(">"):
                tok = _first_token(line)
                if tok:
                    names.append(tok)
    return names


def count_headers(path: Path) -> int:
    n = 0
    with path.open("r") as fh:
        for line in fh:
            if line.startswith(">"):
                n += 1
    return n


def append_fasta(src: Path, dst) -> int:
    """Append src FASTA to an open dst file handle. Returns header count."""
    n = 0
    with src.open("r") as fh:
        for line in fh:
            if line.startswith(">"):
                n += 1
            dst.write(line)
    return n


def main(argv: list[str] | None = None) -> int:
    p = argparse.ArgumentParser(description=__doc__.split("\n")[0])
    p.add_argument("--transcripts", required=True, type=Path, help="transcriptome FASTA")
    p.add_argument("--genome", required=True, type=Path, help="genome FASTA (decoys)")
    p.add_argument("--out-gentrome", required=True, type=Path)
    p.add_argument("--out-decoys", required=True, type=Path)
    p.add_argument("--force", action="store_true", help="overwrite existing outputs")
    args = p.parse_args(argv)

    for f in (args.transcripts, args.genome):
        if not f.is_file():
            print(f"ERROR: input not found: {f}", file=sys.stderr)
            return 2

    for out in (args.out_gentrome, args.out_decoys):
        if out.exists() and not args.force:
            print(f"ERROR: {out} exists (use --force to overwrite)", file=sys.stderr)
            return 2

    tx_headers = count_headers(args.transcripts)
    decoy_names = collect_decoys(args.genome)
    if tx_headers == 0:
        print(f"ERROR: no FASTA records in transcripts file {args.transcripts}", file=sys.stderr)
        return 2
    if not decoy_names:
        print(f"ERROR: no FASTA records in genome file {args.genome}", file=sys.stderr)
        return 2

    # Heuristic swap guard: a transcriptome normally has many more records than a
    # genome assembly. If transcripts has fewer headers than the genome, warn loudly.
    if tx_headers < len(decoy_names):
        print(
            "WARNING: the transcripts FASTA has fewer records "
            f"({tx_headers}) than the genome ({len(decoy_names)}). "
            "Did you swap --transcripts and --genome? Continuing anyway.",
            file=sys.stderr,
        )

    # decoys.txt
    args.out_decoys.write_text("\n".join(decoy_names) + "\n")

    # gentrome.fa = transcripts FIRST, then genome (decoys last)
    with args.out_gentrome.open("w") as dst:
        append_fasta(args.transcripts, dst)
        append_fasta(args.genome, dst)

    print(
        f"Wrote {args.out_decoys} ({len(decoy_names)} decoy names) and "
        f"{args.out_gentrome} (transcripts then genome).",
        file=sys.stderr,
    )
    print(
        "Next: salmon index -t "
        f"{args.out_gentrome} -d {args.out_decoys} -i salmon_index -k 31 -p 8",
        file=sys.stderr,
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())