Add Torino module for annotation-free isoform quantification

code/script/molecular_phenotypes/calling/torino_count_matrix.R

@@ -0,0 +1,155 @@
+#!/usr/bin/env Rscript
+# ============================================================
+# torino_count_matrix.R
+# Per-gene read-depth matrices, the input Torino factorizes, for the torino_2 step
+# of isoform_calling_torino.ipynb. Replaces prepare_count_tables.py from the
+# Torino paper repository (cfbuenabadn/torino_paper, code/scripts).
+#
+# For each gene in --genes the window [start - flank, end + flank) is read from
+# every sample's coverage bedGraph and written as one CSV, samples as rows and
+# bases as columns: a Sample_ID column, then one column per base named
+# <chrom>:<position>, 1-based. This is the layout torino::run_torino_gene reads.
+#
+# Coverage files are found by the names torino_coverage.sh gives them:
+#   <coverage-dir>/<sample_id>.coverage.bedgraph.gz              strand_missing
+#   <coverage-dir>/<sample_id>.{plus,minus}.coverage.bedgraph.gz rf, fr
+# A stranded sample contributes the file matching the gene's strand. Unstranded
+# libraries cannot separate genes that overlap on opposite strands; the Torino
+# paper drops such regions for GTEx, and this script leaves that to the gene list.
+#
+# A gene is written only if at least --min-samples samples have a summed depth
+# (the sum over the window of per-base depth, the quantity Torino's own
+# min_counts filters on) of at least --min-counts. Genes failing the filter are
+# listed, with the reason, in --skipped.
+#
+# --chunk i --n-chunks n processes the i-th of n contiguous slices of --genes, so
+# a cohort-wide run can be split across jobs.
+#
+# Writes <output-dir>/<gene_id>.csv.gz for each kept gene, --output-list (one CSV
+# path per line) and --skipped (gene_id, gene_name, reason; header always written).
+# ============================================================
+
+suppressPackageStartupMessages({
+    library(argparser)
+    library(data.table)
+})
+
+parser <- arg_parser("Per-gene depth matrices for Torino")
+parser <- add_argument(parser, "--samples", type = "character",
+                       help = "sample manifest with sample_id and strand columns")
+parser <- add_argument(parser, "--coverage-dir", type = "character",
+                       help = "directory holding torino_coverage.sh output")
+parser <- add_argument(parser, "--genes", type = "character",
+                       help = "gene BED: chrom, start (0-based), end, gene_id, gene_name, strand")
+parser <- add_argument(parser, "--output-dir", type = "character", help = "directory for the per-gene CSVs")
+parser <- add_argument(parser, "--output-list", type = "character", help = "file listing the CSVs written")
+parser <- add_argument(parser, "--skipped", type = "character", help = "file listing genes not written")
+parser <- add_argument(parser, "--flank", type = "integer", default = 50L,
+                       help = "bases added on each side of the gene")
+parser <- add_argument(parser, "--min-counts", type = "integer", default = 100L,
+                       help = "summed depth a sample needs over the window to count towards --min-samples")
+parser <- add_argument(parser, "--min-samples", type = "integer", default = 10L,
+                       help = "samples that must reach --min-counts for the gene to be written")
+parser <- add_argument(parser, "--chunk", type = "integer", default = 1L, help = "slice of --genes to process")
+parser <- add_argument(parser, "--n-chunks", type = "integer", default = 1L, help = "number of slices")
+argv <- parse_args(parser)
+for (a in c("samples", "coverage_dir", "genes", "output_dir", "output_list", "skipped"))
+    if (is.na(argv[[a]])) stop("--", gsub("_", "-", a), " is required")
+
+die <- function(...) { message("ERROR: ", ...); quit(status = 1) }
+if (argv$n_chunks < 1 || argv$chunk < 1 || argv$chunk > argv$n_chunks)
+    die("--chunk must be between 1 and --n-chunks (got ", argv$chunk, " of ", argv$n_chunks, ")")
+if (argv$flank < 0) die("--flank must be zero or more")
+if (Sys.which("tabix") == "") die("tabix is not on PATH")
+
+# ---- samples and their coverage files
+manifest <- fread(argv$samples, sep = "\t", colClasses = "character")
+for (col in c("sample_id", "strand"))
+    if (!col %in% names(manifest)) die(argv$samples, " has no '", col, "' column")
+if (anyDuplicated(manifest$sample_id)) die("duplicated sample_id in ", argv$samples)
+# the regtools codes other splicing manifests use are accepted as aliases
+alias <- c(XS = "strand_missing", RF = "rf", FR = "fr")
+manifest[strand %in% names(alias), strand := alias[strand]]
+bad <- setdiff(unique(manifest$strand), c("strand_missing", "rf", "fr"))
+if (length(bad)) die("unknown strand value(s) ", paste(bad, collapse = ", "), " in ", argv$samples)
+
+cov_file <- function(id, strand, gene_strand) {
+    tag <- if (strand == "strand_missing") "" else if (gene_strand == "-") ".minus" else ".plus"
+    file.path(argv$coverage_dir, paste0(id, tag, ".coverage.bedgraph.gz"))
+}
+needed <- unlist(lapply(seq_len(nrow(manifest)), function(i)
+    if (manifest$strand[i] == "strand_missing") cov_file(manifest$sample_id[i], "strand_missing", "+")
+    else c(cov_file(manifest$sample_id[i], "rf", "+"), cov_file(manifest$sample_id[i], "rf", "-"))))
+missing <- needed[!file.exists(needed) | !file.exists(paste0(needed, ".tbi"))]
+if (length(missing))
+    die(length(missing), " coverage file(s) or index(es) missing, e.g. ", missing[1])
+
+# ---- genes, this chunk's slice
+gene_lines <- readLines(gzfile(argv$genes))                 # plain or gzipped
+gene_lines <- gene_lines[nzchar(gene_lines) & !startsWith(gene_lines, "#")]
+if (!length(gene_lines)) die(argv$genes, " lists no genes")
+genes <- fread(text = gene_lines, sep = "\t", header = FALSE,
+               colClasses = list(character = c(1, 4, 5, 6), integer = c(2, 3)))
+if (ncol(genes) < 6) die(argv$genes, " needs six columns: chrom, start, end, gene_id, gene_name, strand")
+setnames(genes, 1:6, c("chrom", "start", "end", "gene_id", "gene_name", "strand"))
+if (anyDuplicated(genes$gene_id)) die("duplicated gene_id in ", argv$genes)
+if (!all(genes$strand %in% c("+", "-"))) die("gene strand must be + or - in ", argv$genes)
+per_chunk <- ceiling(nrow(genes) / argv$n_chunks)
+genes <- genes[(seq_len(.N) - 1L) %/% per_chunk + 1L == argv$chunk]
+
+# Reconcile "chr22" and "22" naming against what the coverage files hold.
+seqs <- system2("tabix", c("-l", shQuote(needed[1])), stdout = TRUE)
+genes[, query_chrom := chrom]
+genes[!chrom %in% seqs & paste0("chr", chrom) %in% seqs, query_chrom := paste0("chr", chrom)]
+genes[!chrom %in% seqs & sub("^chr", "", chrom) %in% seqs, query_chrom := sub("^chr", "", chrom)]
+if (nrow(genes) && !any(genes$query_chrom %in% seqs))
+    die("no gene chromosome of this chunk (", paste(head(unique(genes$chrom), 3), collapse = ", "),
+        ") is in the coverage files, which hold ", paste(head(seqs, 3), collapse = ", "), " ...")
+
+# ---- depth of one sample over [ws, we), 0-based half-open
+window_depth <- function(file, chrom, ws, we) {
+    depth <- integer(we - ws)
+    region <- sprintf("%s:%d-%d", chrom, ws + 1L, we)
+    hits <- system2("tabix", c(shQuote(file), shQuote(region)), stdout = TRUE)
+    if (!length(hits)) return(depth)
+    bg <- fread(text = hits, header = FALSE, col.names = c("chrom", "s", "e", "d"))
+    for (r in seq_len(nrow(bg))) {
+        a <- max(bg$s[r], ws); b <- min(bg$e[r], we)
+        if (b > a) depth[(a - ws + 1L):(b - ws)] <- as.integer(bg$d[r])
+    }
+    depth
+}
+
+dir.create(argv$output_dir, recursive = TRUE, showWarnings = FALSE)
+written <- character(0)
+skipped <- data.table(gene_id = character(0), gene_name = character(0), reason = character(0))
+for (g in seq_len(nrow(genes))) {
+    gn <- genes[g]
+    if (!gn$query_chrom %in% seqs) {
+        skipped <- rbind(skipped, data.table(gn$gene_id, gn$gene_name, "chromosome not in coverage files"), use.names = FALSE)
+        next
+    }
+    ws <- max(0L, as.integer(gn$start) - argv$flank)
+    we <- as.integer(gn$end) + argv$flank
+    m <- vapply(seq_len(nrow(manifest)), function(i)
+        window_depth(cov_file(manifest$sample_id[i], manifest$strand[i], gn$strand), gn$query_chrom, ws, we),
+        integer(we - ws))
+    m <- matrix(m, nrow = we - ws)                           # bases x samples, even for one sample
+    passing <- sum(colSums(m) >= argv$min_counts)
+    if (passing < argv$min_samples) {
+        skipped <- rbind(skipped, data.table(gn$gene_id, gn$gene_name,
+            sprintf("%d sample(s) with depth >= %d, %d needed", passing, argv$min_counts, argv$min_samples)),
+            use.names = FALSE)
+        next
+    }
+    dt <- data.table(Sample_ID = manifest$sample_id, t(m))
+    setnames(dt, c("Sample_ID", paste0(gn$chrom, ":", seq.int(ws + 1L, we))))
+    out <- file.path(argv$output_dir, paste0(gn$gene_id, ".csv.gz"))
+    fwrite(dt, out, compress = "gzip")
+    written <- c(written, out)
+}
+
+writeLines(written, argv$output_list)
+fwrite(skipped, argv$skipped, sep = "\t")
+message(sprintf("chunk %d/%d: %d gene(s) written, %d skipped", argv$chunk, argv$n_chunks,
+                length(written), nrow(skipped)))

code/script/molecular_phenotypes/calling/torino_coverage.sh

@@ -0,0 +1,73 @@
+#!/usr/bin/env bash
+# ============================================================
+# torino_coverage.sh
+# Per-base read coverage of one BAM as bgzipped, tabix-indexed bedGraph, for the
+# torino_1 step of isoform_calling_torino.ipynb.
+#
+# Coverage is `bedtools genomecov -bga -split`: spliced reads add depth to the
+# exons they cover and none to the intron they skip, and every base of every
+# reference sequence is reported, zero-depth runs included. -bga rather than -bg
+# lets torino_count_matrix.R tell a zero-coverage stretch from a region absent
+# from the file.
+#
+# --strandness is the manifest value:
+#   strand_missing  one file, <prefix>.coverage.bedgraph.gz, all reads.
+#   rf | fr         two files, <prefix>.plus.coverage.bedgraph.gz and
+#                   <prefix>.minus.coverage.bedgraph.gz, holding the coverage of
+#                   transcripts on that genomic strand. rf is a dUTP library (first
+#                   read antisense to the transcript), fr the opposite. -du flips
+#                   the second mate so both mates of a pair count on one strand.
+#
+# --min-mapq N (default 0, off) keeps only reads with mapping quality >= N; STAR
+# gives uniquely mapped reads 255. Secondary and supplementary alignments are
+# always dropped so each read is counted once.
+# ============================================================
+set -euo pipefail
+
+BAM=""
+PREFIX=""
+STRANDNESS=""
+MIN_MAPQ=0
+
+while [[ $# -gt 0 ]]; do
+    case "$1" in
+        --bam)          BAM="$2"; shift 2 ;;
+        --output-prefix) PREFIX="$2"; shift 2 ;;
+        --strandness)   STRANDNESS="$2"; shift 2 ;;
+        --min-mapq)     MIN_MAPQ="$2"; shift 2 ;;
+        *) echo "Unknown argument: $1" >&2; exit 1 ;;
+    esac
+done
+
+: "${BAM:?--bam is required}"
+: "${PREFIX:?--output-prefix is required}"
+: "${STRANDNESS:?--strandness is required (strand_missing | rf | fr)}"
+case "$MIN_MAPQ" in ""|*[!0-9]*) echo "--min-mapq must be a whole number, got '$MIN_MAPQ'" >&2; exit 1 ;; esac
+case "$STRANDNESS" in
+    strand_missing|rf|fr) ;;
+    *) echo "--strandness must be strand_missing, rf or fr, got '$STRANDNESS'" >&2; exit 1 ;;
+esac
+[ -s "$BAM" ] || { echo "BAM not found or empty: $BAM" >&2; exit 1; }
+
+mkdir -p "$(dirname "$PREFIX")"
+FILTERED="${PREFIX}.filtered.bam"
+trap 'rm -f "$FILTERED"' EXIT
+# 0x100 secondary, 0x800 supplementary, 0x4 unmapped
+samtools view -b -F 0x904 -q "$MIN_MAPQ" -o "$FILTERED" "$BAM"
+
+coverage() {   # coverage <output.bedgraph.gz> [bedtools strand options]
+    local out="$1"; shift
+    bedtools genomecov -bga -split "$@" -ibam "$FILTERED" | bgzip -c > "$out"
+    tabix -f -p bed "$out"
+}
+
+if [ "$STRANDNESS" = strand_missing ]; then
+    coverage "${PREFIX}.coverage.bedgraph.gz"
+else
+    # bedtools -strand selects reads by their own strand; with -du that is the first
+    # mate's strand for a pair. For rf the first read is antisense, so plus-strand
+    # transcripts are read from minus-strand alignments.
+    if [ "$STRANDNESS" = rf ]; then plus_reads="-"; minus_reads="+"; else plus_reads="+"; minus_reads="-"; fi
+    coverage "${PREFIX}.plus.coverage.bedgraph.gz" -strand "$plus_reads" -du
+    coverage "${PREFIX}.minus.coverage.bedgraph.gz" -strand "$minus_reads" -du
+fi

code/script/molecular_phenotypes/calling/torino_factorize.R

@@ -0,0 +1,50 @@
+#!/usr/bin/env Rscript
+# ============================================================
+# torino_factorize.R
+# Torino factorization of one gene's depth matrix, for the torino_3 step of
+# isoform_calling_torino.ipynb. Replaces run_torino_factorization.R from the
+# Torino package (cfbuenabadn/torino, inst/scripts), whose positional arguments
+# and TORINO_FILTER_SAMPLES variable do not fit the pipeline's named-flag wrappers.
+#
+# A thin call to torino::run_torino_gene: samples whose summed depth is below
+# --min-counts are dropped, consecutive bases are summed in groups of --merge-n to
+# bound memory, and the matrix is factorized into --K coverage components by
+# empirical Bayes Poisson matrix factorization.
+#
+# --merge-n auto applies the package's rule for the matrix width: 3 above 15 kb,
+# 9 above 60 kb, 30 above 1 Mb, otherwise 1.
+#
+# Writes --output, an RDS list with $gene, $torino_K<K> and $coords_all (see
+# ?torino::run_torino_gene). Torino returns NULL when the factorization fails,
+# for instance when fewer samples than --K pass --min-counts; that NULL is saved
+# with a warning rather than stopping the run, and torino_usage.R skips it.
+# ============================================================
+
+suppressPackageStartupMessages(library(argparser))
+
+parser <- arg_parser("Torino factorization of one gene")
+parser <- add_argument(parser, "--counts", type = "character", help = "<gene_id>.csv.gz from torino_count_matrix.R")
+parser <- add_argument(parser, "--output", type = "character", help = "RDS to write")
+parser <- add_argument(parser, "--gene", type = "character", default = NA,
+                       help = "gene ID stored in the RDS [default: --counts file name]")
+parser <- add_argument(parser, "--K", type = "integer", default = 3L, help = "number of coverage components")
+parser <- add_argument(parser, "--min-counts", type = "integer", default = 100L,
+                       help = "summed depth below which a sample is dropped")
+parser <- add_argument(parser, "--merge-n", type = "character", default = "auto",
+                       help = "bases summed per column, or auto")
+parser <- add_argument(parser, "--seed", type = "integer", default = 1L, help = "random seed")
+argv <- parse_args(parser)
+for (a in c("counts", "output")) if (is.na(argv[[a]])) stop("--", a, " is required")
+
+die <- function(...) { message("ERROR: ", ...); quit(status = 1) }
+if (is.na(argv$K) || argv$K < 1) die("--K must be a positive integer")
+merge_n <- argv$merge_n
+if (merge_n != "auto" && !grepl("^[1-9][0-9]*$", merge_n)) die("--merge-n must be a positive integer or auto")
+if (!requireNamespace("torino", quietly = TRUE))
+    die("the torino R package is not installed (https://github.com/cfbuenabadn/torino)")
+
+gene <- if (is.na(argv$gene)) sub("\\.csv(\\.gz)?$", "", basename(argv$counts)) else argv$gene
+fit <- torino::run_torino_gene(gene_name = gene, K = argv$K, counts_path = argv$counts,
+                               out_path = argv$output, min_counts = argv$min_counts,
+                               n = if (merge_n == "auto") "auto" else as.integer(merge_n),
+                               seed = argv$seed)