WulfRNA

WulfRNA is a packaged CLI for a focused bulk RNA-seq workflow on paired-end (default) or single-end FASTQ input.

Current outputs:

• source-qualified, gene-level expected counts and TPM from Salmon or kallisto
• optional STAR count and WulfRNA-calculated gene-level TPM matrices
• optional coordinate-sorted STAR BAM and BAI
• STAR splice-junction and gene-count outputs
• final MultiQC report

Pipeline steps: FastQC (raw) → Cutadapt → FastQC (trimmed) → optional STAR alignment → transcript quantification (Salmon or kallisto) → gene-level aggregation → MultiQC.

Public-readiness note: WulfRNA v0.2.2 is a lightweight single-server bulk RNA-seq runner intended for local or HPC workstation use.

1) Prerequisites

WulfRNA does not install analysis binaries for you. They must already be installed and visible in PATH at runtime.

Required commands:

• bash
• python3
• fastqc
• 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

2) Install the packaged CLI

Use Python 3.10+ and a virtual environment:

python3 -m venv .venv
source .venv/bin/activate
pip install .

This installs the wulfrna command.

3) Required input layout

Work directory

WORKDIR must contain fastq/ and uses an explicit input layout:

paired_end (default)

• sample_id_R1.fastq.gz
• sample_id_R2.fastq.gz

single_end (enable with --single-end or --SE)

• sample_id.fastq.gz or sample_id_R1.fastq.gz
• *_R2.fastq.gz files are not allowed in single-end mode (run fails if present).

Reference directory

--reference must point to a directory containing backend-specific resources plus a shared tx2gene map:

• shared: combined_tx2gene.tsv
• shared: combined_gene_annotation.tsv, with columns gene_id, GeneName, and positive integer gene_length_bp (the union length of annotated exons, counting overlaps once)
• 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

Build the shared gene annotation

Generate combined_gene_annotation.tsv from the same GTF used to build the STAR index and combined_tx2gene.tsv:

python scripts/build_gene_annotation.py \
  --gtf /path/to/combined.gtf \
  --output /path/to/reference/combined_gene_annotation.tsv

The builder reads exon records, preserves each GTF gene_id exactly, merges overlapping or adjacent exon intervals within a gene, and writes their union length as gene_length_bp. It uses gene_name for GeneName and writes NA when no symbol exists. A gene assigned to multiple chromosomes, inconsistent gene names, malformed records, or invalid coordinates cause a clear failure. Use the identical combined GTF for all reference artifacts so gene IDs match between this table, combined_tx2gene.tsv, and STAR.

4) Run command

wulfrna run WORKDIR --reference /path/to/reference --stranded reverse --threads 16

Arguments:

• WORKDIR (positional): working directory containing fastq/
• --reference PATH (required): prepared reference directory
• --stranded {none|forward|reverse} (required): library strandedness
  • Salmon mapping: none -> IU, forward -> ISF, reverse -> ISR
  • 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,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.

5) Minimal smoke test plan

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 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

6) Main outputs and status markers

Expected primary outputs on full success:

• abundance/<quantifier>_gene_expected_counts.tsv
• abundance/<quantifier>_gene_tpm.tsv
• abundance/gene_expected_counts.tsv (legacy v0.2.x compatibility output without GeneName)
• abundance/gene_tpm.tsv (legacy v0.2.x compatibility output without GeneName)
• 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
• abundance/star_gene_counts/<sample>.star.ReadsPerGene.out.tab (metadata-preserving copy)
• abundance/star_gene_counts.tsv
• abundance/star_gene_tpm.tsv

New matrices begin with gene_id, GeneName, then deterministically sorted sample columns. Salmon/kallisto TPM is summed from transcript-quantifier output. In contrast, star_gene_tpm.tsv is WulfRNA-calculated gene-level TPM from STAR counts using exon-union lengths from combined_gene_annotation.tsv; it is not a native STAR output. STAR's count column is selected from --stranded (none: column 2, forward: column 3, reverse: column 4). STAR TPM values are serialized with six decimal places, then each serialized sample column is required to sum to 1,000,000 within an absolute tolerance of 0.01.

Matrix behavior note:

• Gene-level matrices include only genes observed in the transcript quantification input (unobserved zero-only genes are not emitted).

Status files in WORKDIR/status/:

• during run: RUNNING
• full success: SUCCESS, finished_at.txt, summary.txt
• dry-run success: DRY_RUN_OK, finished_at.txt, summary.txt
• failure: FAILED, failed_step.txt (and optionally failed_sample.txt), summary.txt
• phase checkpoint markers: status/steps/<phase>.done for each completed phase
• run compatibility manifest: status/manifest.json

7) Automatic resume behavior

By default, wulfrna run ... automatically resumes at the phase level (not sample-level).

Phases:

• fastqc_raw
• cutadapt
• fastqc_trimmed
• align (only when --aligner star)
• quant
• aggregate
• multiqc

A phase is skipped only when:

• status/steps/<phase>.done exists, and
• expected outputs for that phase exist and are non-empty, and
• status/manifest.json is compatible with the current run configuration.

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 combined_gene_annotation.tsv fingerprint changed: aggregate and multiqc rerun.
• If only thread count changed: resume is allowed.

To fully disable resume:

wulfrna run WORKDIR --reference REFDIR --stranded reverse --threads 16 --no-resume

To rerun only aggregation/reporting after fixing combined_tx2gene.tsv:

• preferred: --force-from aggregate, or
• remove aggregate + multiqc outputs and rerun.

Example:

wulfrna run WORKDIR --reference REFDIR --stranded reverse --threads 16 --force-from aggregate

8) Citation

If you use WulfRNA in your research, please cite:

Kudoh, F. (2026). WulfRNA (Version 0.2.2) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.21409117

Citation metadata is also available in CITATION.cff.