End-to-end bulk RNA-seq orchestrator — takes raw FASTQ reads through QC and trimming (FastQC, fastp/Trim Galore), alignment and quantification (STAR, Salmon, featureCounts), assembles a gene-level counts matrix, then hands off to differential expression (pydeseq2), pathway/GSEA enrichment (pathway-e

Bulk RNA seq Overview This skill prepares bulk RNA seq reads and quantification output for a reproducible gene level comparison. It owns sample validation, reads to counts recipes, and the count/metadata handoff; downstream statistical fitting and enrichment remain in their specialist skills. "Defensible" means three things, applied throughout: Reproducible — pinned pipeline/tool versions, containers where possible, recorded parameters, fixed random seeds. Quality gated — QC is inspected and acted on before, during, and after quantification, not skipped. Statistically sound — adequate replication, a design that matches the biology, counts handled correctly, and FDR controlled testing. The pipeline is: FastQC/trim → align/quant (STAR/Salmon) → counts → DE (pydeseq2) → enrichment (pathway enrichment) → figures . When to Use This Skill Use this skill when the user wants to: Go from FASTQ files (or a sequencing run) to differentially expressed genes and pathways. Run or configure nf core/rnaseq , or align/quantify with STAR, Salmon, or featureCounts. Turn Salmon/STAR/featureCounts output into a counts matrix ready for DESeq2/PyDESeq2. Design or sanity check a bulk RNA seq experiment (replicates, batch, strandedness) before committing compute. Scope an end to end RNA seq analysis and decide which tools and skills to chain. This is bulk RNA seq (samples = biological specimens). For single cell/nuclei data use scanpy ; for the DE statistics alone use pydeseq2 ; for enrichment alone use pathway enrichment . The Pipeline at a Glance Two Upstream Paths — Pick One The reads → counts stage can be run two ways. Both produce gene level counts, but STAR/featureCounts and Salmon/tximport do not generally give identical or interchangeable values: assignment rules, multimapping, and effective length corrections differ. Choose one quantification route for the entire comparison. Use Path A — nf core/rnaseq when… Use Path B — standalone tools when… You want the field standard, audited, citable pipeline with one command You have a few samples and want to learn/inspect each step Many samples, or you'll scale to HPC/cloud No Nextflow/containers available, or a constrained environment Reproducibility and a full MultiQC report matter most You need a non standard step the pipeline doesn't expose → Drive it through the nextflow skill → Follow references/upstream manual.md When unsure, prefer Path A : nf core/rnaseq already wires together FastQC → trimming → STAR/Salmon → quantification → tximport → MultiQC with sensible, reviewed defaults, which is the most defensible option. Path B exists for transparency and constrained setups. Both paths converge on a gene level counts matrix . Preserve whether counts are raw or length scaled, the transcript to gene mapping release, and any offsets. Length scaled counts must not receive a second transcript length correction. Setup Reviewed against nf core/rnaseq 3.27.0, STAR 2.7.11b documentation, Salmon 2.8.0, pytximport 0.13.0 and PyDESeq2 0.5.4. The bundled Python bridge and a tiny Salmon index/quant run were executed; the full Nextflow/STAR/trimming recipes are illustrative, not an end to end validation. Salmon 2.x cannot read older C++ indexes: rebuild them. The bridge's length scaled Salmon route is for full length bulk RNA seq; 3′ counting assays need original counts without transcript length correction. Record the exact versions you use (pipeline revision, tool versions, reference genome + annotation release) — they belong in the methods section and make the analysis reproducible. Version 2.0 makes STAR strandedness explicit and rejects ambiguous samples, gene sets, transcript mappings and fractional featureCounts data that older bridge versions accepted. Quick Start Path A — nf core/rnaseq (recommended) nf core/rnaseq runs tximport internally, so gene counts come out already merged — no bridge script needed. Use results/star salmon/salmon.merged.gene counts length scaled.tsv for DE. Samplesheet format, aligner choice, and outputs: references/upstream nfcore.md . For engine/HPC/cloud/container detail, use the nextflow skill. Path B — standalone STAR/Salmon (abbreviated) Full recipes (FastQC, fastp/Trim Galore, STAR index+align+ quantMode GeneCounts , Salmon decoy aware index, featureCounts, strandedness): references/upstream manual.md . Counts → DE → enrichment (both paths) Stage by Stage Workflow Work top to bottom. Each stage names the skill or file that owns the detail. Don't skip the design/QC stages — they are where bulk RNA seq studies most often go wrong. 1. Design & sample sheet. Confirm ≥3 biological replicates per group, identify batch/confounders, and choose the comparison(s). Build the samplesheet and validate it with scripts/validate samplesheet.py . Rationale and rules: references/design and qc.md . 2. Raw read QC. FastQC per file; aggregate with MultiQC. Check per base quality, adapter content, duplication, and over representation. Thresholds: references/design and qc.md . 3. Trimming. Remove adapters and low quality tails (via fastp or Trim Galore ). Re run FastQC to confirm. Recipes: references/upstream manual.md (Path A does this for you). 4. Align / quantify. STAR (genome alignment + quantMode GeneCounts ) and/or Salmon (decoy aware selective alignment). Determine strandedness — the wrong convention can silently discard most assigned reads. Detail: references/upstream manual.md ; pipeline params: references/upstream nfcore.md . 5. Build the counts matrix. Turn quant output into a gene × sample integer matrix and a metadata template ( scripts/build counts matrix.py ). The estimated count and gene ID mapping nuances live in references/counts and handoff.md . 6. Differential expression → pydeseq2 skill. Load counts.csv + metadata.csv , set the design (e.g. ~batch + condition ), check full rank and residual degrees of freedom, fit, and test an explicit contrast (e.g. treated versus control) with FDR control. Inspect the PCA and p value histogram as QC. 7. Enrichment → pathway enrichment skill. For GSEA, rank the full gene list by the DESeq2 stat ; for ORA, pass the thresholded hit list (padj < 0.05, optionally log2FC 1). Match identifiers to the selected library; retain an auditable mapping and an assay specific ORA background. 8. Figures → scientific visualization skill. Volcano, MA, sample distance heatmap, PCA, and enrichment dotplots, plus the MultiQC report for the QC narrative. The counts → DE bridge (the key glue) This is the one stage with no upstream/downstream skill, so this skill owns it. scripts/build counts matrix.py converts quant output into exactly what pydeseq2 expects: Salmon ( from salmon ): aggregates per sample quant.sf to gene level with pytximport using counts from abundance="length scaled tpm" (an offset free choice for full length gene level DE), needs a tx2gene map. STAR ( from star ): reads each ReadsPerGene.out.tab , selecting the column for your strandedness (unstranded/forward/reverse). featureCounts ( from featurecounts ): parses the combined featureCounts matrix. It writes counts.csv (genes × samples, integers), counts provenance.json (input hashes, import mode and sample order), and metadata template.csv (one row per sample) for you to fill in. Salmon/RSEM counts are estimates (non integer); they are rounded to integers because PyDESeq2 requires integer counts — see references/counts and handoff.md for why this is acceptable with length scaled tpm and how it differs from the offset based DESeq2+tximport route. That reference also covers identifier mapping (when required by the selected enrichment library) and the exact orientation PyDESeq2 wants. Common Pitfalls These cause most wrong or irreproducible bulk RNA seq results: 1. Too few replicates. Plan biological replication from expected variability and effect size; three per group is a starting point, not a power guarantee. Technical lanes do not increase biological sample size. 2. Confounded batch and condition. If every treated sample was processed on a different day/lane than controls, the effect is unrecoverable. Randomize, and model known batches ( ~batch + condition ). See references/design and qc.md . 3. Wrong strandedness. Choosing the wrong STAR column or featureCounts s /Salmon library type can discard most assigned reads; there is no universal 50% loss. Use Salmon l A or infer strandedness, and verify the assigned reads fraction. 4. Feeding TPM/FPKM to DESeq2. DESeq2 needs raw (or length scaled) counts , never TPM/FPKM/normalized values. The bridge handles this. 5. Non integer counts. The bridge rounds length scaled Salmon estimates only; it rejects fractional featureCounts values rather than truncating them. Its Salmon route is for full length assays, not 3′ tag counts. 6. Gene ID mismatch into enrichment. Match gene IDs to the selected organism and gene set release. Many GMT libraries use symbols; g:Profiler can accept Ensembl IDs directly. Do not collapse ambiguous mappings silently. 7. Skipping post quant QC. Always look at the PCA and sample distance heatmap before trusting DE — they expose swapped labels, outliers, and hidden batches. 8. Mixing aligners across samples. Quantify every sample with the same tool, version, reference, and parameters. 9. Unpinned versions. "latest" pipelines/genomes make results unreproducible; pin r , tool versions, and the genome/annotation release. Integration with Other Skills Upstream execution: nextflow (runs nf core/rnaseq , Path A; HPC/cloud/containers). Reference data / gene IDs: gget ( gget ref for genome+GTF, gget info / gget search for ID mapping), database lookup (Ensembl/NCBI), biopython / pysam (FASTA/BAM handling). Differential expression: pydeseq2 (the DE engine this skill hands counts to). Enrichment: pathway enrichment (ORA + GSEA; its scripts/run enrichment.py reads a DESeq2 results CSV directly). Figures & reporting: scientific visualization , matplotlib , seaborn ; scientific writing for the methods/results narrative. Related but distinct: scanpy (single cell), statistical analysis (multiple testing depth). Reference Files Read the relevant file when you need depth — each is self contained: references/upstream nfcore.md — Path A: samplesheet format, aligner / pseudo aligner choice, key params, the salmon.merged.gene counts .tsv outputs, MultiQC, and what to hand to pydeseq2 . references/upstream manual.md — Path B: FastQC, fastp/Trim Galore, STAR genome index + alignment + quantMode GeneCounts , Salmon decoy aware index + quant , featureCounts, and how to determine strandedness. references/counts and handoff.md — turning quant output into PyDESeq2 ready counts.csv / metadata.csv (pytximport, STAR column selection, featureCounts), the integer/estimated count nuance, Ensembl→symbol mapping, and the DE→enrichment rank/hit list recipe. references/design and qc.md — experimental design (replication, batch, confounding, design formulas) and QC metric interpretation (mapping rate, duplication, rRNA, complexity, PCA/outliers) — the defensible pipeline backbone. Resources nf core/rnaseq: https://nf co.re/rnaseq · STAR: https://github.com/alexdobin/STAR · Salmon: https://combine lab.github.io/salmon/ fastp: https://github.com/OpenGene/fastp · Trim Galore: https://github.com/FelixKrueger/TrimGalore · MultiQC: https://multiqc.info pytximport: https://pytximport.complextissue.com · featureCounts (Subread): https://subread.sourceforge.net Method background: Love et al. 2014 (DESeq2) DOI 10.1186/s13059 014 0550 8 · Soneson et al. 2015 (tximport) DOI 10.12688/f1000research.7563.2 Citing Scientific Agent Skills This skill is part of Scientific

bulk-rnaseq