Add draft fihspond rule

Author: guillaumecharbonnier

src/r/scripts/fishpond_tximeta.R

@@ -0,0 +1,68 @@
+# fishpond_tximeta.R
+# Usage: Rscript fishpond_tximeta.R <coldata_tsv> <outdir> [bfc_path]
+#
+# coldata_tsv : TSV file with columns 'name' (sample label) and 'file'
+#               (path to quant.sf produced by salmon with Gibbs samples)
+# outdir      : output directory; receives transcripts.rds and genes.rds
+# bfc_path    : (optional) BiocFileCache directory for tximeta reference
+#               downloads; default: out/TximetaBFC
+
+suppressPackageStartupMessages({
+    library(tximeta)
+    library(fishpond)
+    library(org.Hs.eg.db)
+    library(SummarizedExperiment)
+    library(data.table)
+})
+
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) < 2) {
+    stop("Usage: Rscript fishpond_tximeta.R <coldata_tsv> <outdir> [bfc_path]")
+}
+
+coldata_tsv <- args[1]
+outdir      <- args[2]
+bfc_path    <- if (length(args) >= 3) args[3] else "out/TximetaBFC"
+
+dir.create(outdir,   showWarnings = FALSE, recursive = TRUE)
+dir.create(bfc_path, showWarnings = FALSE, recursive = TRUE)
+
+setTximetaBFC(bfc_path)
+
+coldata <- as.data.frame(fread(coldata_tsv, sep = "\t"))
+
+if (!all(c("name", "file") %in% names(coldata))) {
+    stop("coldata TSV must contain 'name' and 'file' columns")
+}
+
+# tximeta expects 'names' and 'files' columns
+coldata$names <- coldata$name
+coldata$files <- coldata$file
+
+missing <- !file.exists(coldata$files)
+if (any(missing)) {
+    warning(
+        sum(missing), " quant.sf file(s) not found and will be skipped:\n",
+        paste(coldata$files[missing], collapse = "\n")
+    )
+    coldata <- coldata[!missing, ]
+}
+
+# Transcript-level SummarizedExperiment (GRCh38 Ensembl; release detected
+# automatically from the salmon index metadata embedded in quant.sf)
+se <- tximeta(coldata)
+
+# Add gene-level SYMBOL annotations
+se <- addIds(se, "SYMBOL", gene = TRUE)
+
+# Scale inferential replicates (Gibbs samples) prior to downstream testing
+se <- scaleInfReps(se)
+
+saveRDS(se, file.path(outdir, "transcripts.rds"))
+message("Saved: ", file.path(outdir, "transcripts.rds"))
+
+# Gene-level summary
+gse <- summarizeToGene(se)
+
+saveRDS(gse, file.path(outdir, "genes.rds"))
+message("Saved: ", file.path(outdir, "genes.rds"))

src/snakemake/rules/fishpond_tximeta.smk

@@ -0,0 +1,57 @@
+"""
+Aim:
+    Rules to import salmon quant.sf files with tximeta (human GRCh38 / Ensembl)
+    and run fishpond preprocessing, producing:
+      - transcript-level SummarizedExperiment (transcripts.rds)
+      - gene-level SummarizedExperiment      (genes.rds)
+
+Usage:
+    Register a coldata TSV in mwconf['ids'] under a chosen <coldata_id> key.
+    The TSV must have two tab-separated columns:
+        name    sample label (used as column name in the SE)
+        file    path to the corresponding salmon quant.sf
+
+    Then request:
+        out/fishpond/tximeta_GRCh38_ensembl/<coldata_id>/transcripts.rds
+        out/fishpond/tximeta_GRCh38_ensembl/<coldata_id>/genes.rds
+"""
+
+
+def _fishpond_tximeta_quant_files(wildcards):
+    """Return the list of quant.sf paths declared in the coldata TSV."""
+    coldata_path = mwconf['ids'][wildcards.coldata_id]
+    with open(coldata_path) as fh:
+        reader = csv.DictReader(fh, delimiter='\t')
+        return [row['file'] for row in reader]
+
+
+rule fishpond_tximeta_GRCh38_ensembl:
+    """
+    Created:
+        2026-02-22
+    Aim:
+        Import salmon quant.sf files (human GRCh38, Ensembl release detected
+        automatically from salmon index metadata) with tximeta, add SYMBOL IDs,
+        scale inferential replicates, and summarise to gene level.
+    Input coldata_id:
+        Key in mwconf['ids'] that points to a TSV with 'name' and 'file' columns.
+    Test:
+        out/fishpond/tximeta_GRCh38_ensembl/my_coldata_id/transcripts.rds
+        out/fishpond/tximeta_GRCh38_ensembl/my_coldata_id/genes.rds
+    """
+    input:
+        coldata  = lambda wildcards: mwconf['ids'][wildcards.coldata_id],
+        quant_sf = _fishpond_tximeta_quant_files,
+        script   = "src/r/scripts/fishpond_tximeta.R"
+    output:
+        transcripts = "out/fishpond/tximeta_GRCh38_ensembl/{coldata_id}/transcripts.rds",
+        genes       = "out/fishpond/tximeta_GRCh38_ensembl/{coldata_id}/genes.rds"
+    params:
+        outdir   = "out/fishpond/tximeta_GRCh38_ensembl/{coldata_id}",
+        bfc_path = "out/TximetaBFC"
+    wildcard_constraints:
+        coldata_id = r"[-a-zA-Z0-9_]+"
+    conda:
+        "../envs/r_fishpond.yaml"
+    shell:
+        "Rscript {input.script} {input.coldata} {params.outdir} {params.bfc_path}"