Add source-qualified annotated gene matrices and STAR-derived TPM output

[Fuki-Kudoh]

Summary

Extend the aggregate phase to produce source-qualified, analysis-ready gene-level matrices for Salmon, kallisto, and STAR.

The new outputs should:

• identify the quantification source in each filename;

• include both Ensembl gene_id and GeneName;

• collect per-sample STAR ReadsPerGene.out.tab files under abundance/;

• aggregate STAR counts into a gene-by-sample matrix;

• calculate a STAR-derived gene-level TPM matrix using gene counts and exon-union gene lengths;

• preserve the current unqualified output files temporarily for backward compatibility.

Suggested release: v0.2.2

Motivation

WulfRNA currently produces:

abundance/gene_expected_counts.tsv
abundance/gene_tpm.tsv

for Salmon or kallisto, while STAR gene counts remain only in:

align/star/<sample_id>/ReadsPerGene.out.tab

This creates several usability problems:

1. The matrix filenames do not indicate whether the values came from Salmon or kallisto.

2. STAR gene counts are not collected or aggregated under abundance/.

3. Gene-level matrices contain gene_id but not an analysis-friendly gene symbol column.

4. STAR runs do not currently produce a TPM matrix.

5. Downstream scripts must manually locate and combine STAR outputs.

Proposed output contract

Salmon

When --quantifier salmon is selected:

abundance/
├── salmon/
│   └── <sample_id>/quant.sf
├── salmon_gene_expected_counts.tsv
└── salmon_gene_tpm.tsv

Kallisto

When --quantifier kallisto is selected:

abundance/
├── kallisto/
│   └── <sample_id>/abundance.tsv
├── kallisto_gene_expected_counts.tsv
└── kallisto_gene_tpm.tsv

STAR

When --aligner star is selected, additionally create:

abundance/
├── star_gene_counts/
│   ├── <sample_id>.star.ReadsPerGene.out.tab
│   └── ...
├── star_gene_counts.tsv
└── star_gene_tpm.tsv

The original STAR outputs must remain unchanged under:

align/star/<sample_id>/

The files under abundance/star_gene_counts/ should be copies of the original STAR files, not replacements.

Matrix format

All new gene-level matrices must use the following leading columns:

gene_id    GeneName    <sample_1>    <sample_2>    ...

Example:

gene_id                 GeneName    sample01    sample02
ENSMUSG00000000001.5    Gnai3       120.0       98.0
ENSMUSG00000000003.14   Pbsn        0.0         4.0

Requirements:

• gene_id is the canonical row identifier.

• Ensembl version suffixes must be preserved when present.

• GeneName is annotation metadata only and must not be used as the unique key.

• Duplicate GeneName values are allowed.

• Missing gene symbols may be written as NA.

• Sample columns must remain deterministically sorted by sample ID.

Shared gene annotation resource

Add a gene-level annotation file to the reference contract:

combined_gene_annotation.tsv

Required columns:

gene_id    GeneName    gene_length_bp

Example:

gene_id                 GeneName    gene_length_bp
ENSMUSG00000000001.5    Gnai3       3264
ENSMUSG00000000003.14   Pbsn        842

Definitions:

• gene_id: must match the gene IDs used by combined_tx2gene.tsv and STAR.

• GeneName: gene symbol; NA is allowed when no symbol is available.

• gene_length_bp: positive integer representing the union length of all annotated exon intervals for that gene.

The exon-union length must count overlapping exonic bases only once within each gene.

combined_tx2gene.tsv should remain the transcript-to-gene mapping resource. Gene-level metadata should be read from combined_gene_annotation.tsv rather than duplicated for every transcript.

Reference validation

Validate combined_gene_annotation.tsv before execution.

The run must fail with a clear reference_check error when:

• the file is missing;

• a required column is missing;

• a gene_id is blank;

• a gene_id occurs more than once;

• gene_length_bp is missing, non-numeric, zero, or negative;

• the same gene cannot be matched consistently across the required reference resources.

A missing gene symbol is not fatal and should be normalized to NA.

Salmon and kallisto aggregation

Rename the new canonical matrices according to the selected quantifier:

abundance/<quantifier>_gene_expected_counts.tsv
abundance/<quantifier>_gene_tpm.tsv

Specifically:

abundance/salmon_gene_expected_counts.tsv
abundance/salmon_gene_tpm.tsv
abundance/kallisto_gene_expected_counts.tsv
abundance/kallisto_gene_tpm.tsv

Add GeneName as the second column using combined_gene_annotation.tsv.

The underlying aggregation behavior should otherwise remain unchanged:

• Salmon expected counts come from NumReads.

• Salmon TPM comes from TPM.

• Kallisto expected counts come from est_counts.

• Kallisto TPM comes from tpm.

• Transcript values are summed by gene_id through combined_tx2gene.tsv.

STAR raw-count collection

During the aggregate phase, for every sample copy:

align/star/<sample_id>/ReadsPerGene.out.tab

to:

abundance/star_gene_counts/<sample_id>.star.ReadsPerGene.out.tab

Requirements:

• preserve the complete native STAR file;

• retain the four STAR summary rows;

• preserve file contents exactly;

• use shutil.copy2() or equivalent metadata-preserving copying;

• fail clearly if an expected source file is missing or empty.

STAR count matrix

Create:

abundance/star_gene_counts.tsv

Select the STAR count column according to WulfRNA --stranded:

WulfRNA option | STAR column | Meaning -- | -- | -- none | 2 | unstranded counts forward | 3 | first-read strand aligned with RNA reverse | 4 | second-read strand aligned with RNA

Exclude the four STAR summary rows:

N_unmapped
N_multimapping
N_noFeature
N_ambiguous

The resulting matrix must contain:

gene_id    GeneName    <sample_1>    <sample_2>    ...

STAR count values should be written as integers.

The aggregate step must fail if:

• a STAR row does not contain four columns;

• a selected count is not an integer;

• duplicate gene IDs occur within one sample;

• sample STAR files contain inconsistent gene sets;

• a STAR gene ID is absent from combined_gene_annotation.tsv.

STAR-derived TPM

Create:

abundance/star_gene_tpm.tsv

STAR itself supplies the gene-level counts. WulfRNA should calculate gene-level TPM using gene_length_bp from combined_gene_annotation.tsv.

For each gene (g) in sample (s):

RPK[g,s] = count[g,s] / (gene_length_bp[g] / 1000)

Then:

TPM[g,s] = RPK[g,s] / sum(RPK[,s]) * 1,000,000

Requirements:

• use the same strandedness-selected STAR counts as star_gene_counts.tsv;

• use exon-union gene length;

• write TPM values with six decimal places;

• retain zero-count genes as 0.000000;

• require a positive gene length for every STAR gene;

• fail clearly if the total RPK for a sample is zero;

• verify that each non-empty sample TPM column sums to approximately 1,000,000 within a documented floating-point tolerance.

The STAR TPM output must be described as WulfRNA-calculated gene-level TPM from STAR counts, not as a native STAR output.

Backward compatibility

Temporarily preserve the existing files:

abundance/gene_expected_counts.tsv
abundance/gene_tpm.tsv

For v0.2.x:

• these files should continue to be generated;

• they should retain the current legacy structure of gene_id followed by sample columns;

• they should contain the same numerical values as the selected Salmon or kallisto canonical matrices;

• they should not include GeneName, to avoid breaking existing downstream scripts;

• documentation should mark them as legacy compatibility outputs.

The source-qualified files should be the new canonical outputs.

Removal of the legacy filenames may be considered separately for v0.3.0.

Aggregate-phase behavior

The aggregate phase should perform the following operations:

1. Load combined_tx2gene.tsv.

2. Load and validate combined_gene_annotation.tsv.

3. Aggregate the selected transcript quantifier.

4. Write source-qualified expected-count and TPM matrices.

5. Write the legacy unqualified matrices.

6. When --aligner star is enabled:

   • copy the native STAR count files;

   • build the STAR count matrix;

   • calculate the STAR TPM matrix.

7. Validate all expected outputs before writing status/steps/aggregate.done.

Resume behavior

Update phase_outputs() so that aggregate completion requires the new outputs.

For Salmon without STAR:

abundance/salmon_gene_expected_counts.tsv
abundance/salmon_gene_tpm.tsv
abundance/gene_expected_counts.tsv
abundance/gene_tpm.tsv

For kallisto without STAR:

abundance/kallisto_gene_expected_counts.tsv
abundance/kallisto_gene_tpm.tsv
abundance/gene_expected_counts.tsv
abundance/gene_tpm.tsv

When STAR is enabled, additionally require:

abundance/star_gene_counts.tsv
abundance/star_gene_tpm.tsv
abundance/star_gene_counts/<sample_id>.star.ReadsPerGene.out.tab

An existing v0.2.1 work directory should therefore rerun from aggregate when the new files are absent, without rerunning trimming, alignment, or transcript quantification.

Changes to combined_gene_annotation.tsv must force rerunning aggregate and downstream phases.

Add its path and fingerprint to:

status/manifest.json

Summary and final validation

Update status/summary.txt and final output validation to report the relevant source-qualified matrices.

When STAR is enabled, include:

abundance/star_gene_counts.tsv
abundance/star_gene_tpm.tsv

Tests

Add unit and integration tests covering at least the following cases.

Annotation parsing

• valid combined_gene_annotation.tsv;

• missing required columns;

• duplicate gene IDs;

• missing GeneName;

• invalid gene length;

• gene ID mismatch.

Salmon

• source-qualified expected-count matrix is created;

• source-qualified TPM matrix is created;

• both include gene_id and GeneName;

• legacy matrices remain available and unchanged in structure.

Kallisto

• source-qualified expected-count matrix is created;

• source-qualified TPM matrix is created;

• both include gene_id and GeneName;

• legacy matrices remain available and unchanged in structure.

STAR

• native STAR files remain under align/star/<sample_id>/;

• copied files are created under abundance/star_gene_counts/;

• copied file contents match the originals;

• count matrix is generated for none, forward, and reverse;

• four STAR summary rows are excluded;

• GeneName is correctly added;

• TPM is calculated from known counts and known exon-union lengths;

• TPM columns sum to approximately 1,000,000;

• zero-count genes remain zero;

• missing or invalid gene lengths produce a clear error.

Resume

• a completed v0.2.1-style work directory reruns aggregate;

• STAR alignment is not rerun when only aggregate outputs are missing;

• changing combined_gene_annotation.tsv reruns aggregate;

• a second unchanged invocation skips the completed aggregate phase.

Documentation

Update:

• README.md;

• docs/spec.md;

• CHANGELOG.md;

• reference preparation documentation;

• output examples;

• the explanation of STAR-derived TPM.

Clearly distinguish:

• transcript-quantifier TPM from Salmon or kallisto;

• WulfRNA-calculated gene-level TPM based on STAR counts and exon-union lengths.

Acceptance criteria

• Source-qualified Salmon matrix filenames are implemented.

• Source-qualified kallisto matrix filenames are implemented.

• Every new matrix includes gene_id and GeneName.

• combined_gene_annotation.tsv is validated and fingerprinted.

• STAR native count files are copied during aggregate.

• STAR counts are aggregated into star_gene_counts.tsv.

• STAR TPM is calculated into star_gene_tpm.tsv.

• STAR strandedness column selection follows --stranded.

• STAR summary rows are excluded from matrices.

• Legacy unqualified Salmon/kallisto matrices remain available.

• Resume behavior recognizes all new aggregate outputs.

• Tests cover Salmon, kallisto, STAR counts, STAR TPM, annotation validation, and resume.

• README, spec, and changelog are updated.