Adding leafcutter2 module

code/script/molecular_phenotypes/calling/leafcutter2_cluster.sh

@@ -0,0 +1,94 @@
+#!/usr/bin/env bash
+# ============================================================
+# leafcutter2_cluster.sh
+# Cluster, classify and quantify junctions with LeafCutter2, for the leafcutter2_2
+# and leafcutter2_cluster steps of splicing_calling_leafcutter2.ipynb.
+#
+# Three checks wrap the leafcutter2 call, because each failure they catch is
+# otherwise silent or surfaces much later:
+#   1. every junction file has the same column count. BED12 (regtools .junc) and
+#      BED6 (star2junc) need different --offset values, so a mixed list cannot be right.
+#   2. junction and annotation chromosome naming agree (chr22 against 22). If they
+#      differ, nothing classifies and every junction simply comes back unannotated.
+#   3. at least one cluster survived the read filters. If none did, clustering still
+#      reports success and the failure appears later, in classification.
+#
+# --keepleafcutter1 is always passed so the LeafCutter1-format tables that
+# leafcutter_ds.R reads are written alongside.
+# ============================================================
+# No pipefail: the checks read only the first lines (zcat | head), which ends zcat early.
+set -eu
+
+JUNCFILES=""
+OUTPREFIX=""
+RUNDIR=""
+ANNOT=""
+GENOME=""
+MAXINTRONLEN=500000
+MINCLUREADS=30
+MINREADS=5
+MINREADSTD=0.5
+MINCLURATIO=0.001
+OFFSET=""
+
+while [[ $# -gt 0 ]]; do
+  case "$1" in
+    --juncfiles)    JUNCFILES="$2"; shift 2 ;;
+    --outprefix)    OUTPREFIX="$2"; shift 2 ;;
+    --rundir)       RUNDIR="$2"; shift 2 ;;
+    --annot)        ANNOT="$2"; shift 2 ;;
+    --genome)       GENOME="$2"; shift 2 ;;
+    --maxintronlen) MAXINTRONLEN="$2"; shift 2 ;;
+    --minclureads)  MINCLUREADS="$2"; shift 2 ;;
+    --minreads)     MINREADS="$2"; shift 2 ;;
+    --minreadstd)   MINREADSTD="$2"; shift 2 ;;
+    --mincluratio)  MINCLURATIO="$2"; shift 2 ;;
+    --offset)       OFFSET="$2"; shift 2 ;;
+    *) echo "Unknown argument: $1" >&2; exit 1 ;;
+  esac
+done
+
+: "${JUNCFILES:?--juncfiles is required}"
+: "${OUTPREFIX:?--outprefix is required}"
+: "${RUNDIR:?--rundir is required}"
+: "${ANNOT:?--annot is required}"
+: "${GENOME:?--genome is required}"
+: "${OFFSET:?--offset is required (0 for BED6, 1 for regtools .junc)}"
+
+# 1. One offset applies to the whole run.
+ncol=$(while read -r f; do [ -n "$f" ] && zcat -f "$f" | head -n1 | awk '{print NF}'; done < "$JUNCFILES" | sort -u | tr '\n' ' ')
+if [ "$(echo $ncol | wc -w)" -gt 1 ]; then
+  echo "ERROR: junction files mix formats (column counts: $ncol). A single --offset cannot suit both." >&2
+  exit 1
+fi
+
+# 2. Chromosome naming.
+jchr=$(zcat -f "$(head -n1 "$JUNCFILES")" | head -n1 | cut -f1)
+achr=$(zcat -f "$ANNOT" | grep -v '^#' | head -n1 | cut -f1)
+case "$jchr" in chr*) jp=prefixed;; *) jp=bare;; esac
+case "$achr" in chr*) ap=prefixed;; *) ap=bare;; esac
+if [ "$jp" != "$ap" ]; then
+  echo "ERROR: chromosome naming differs. Junctions use '$jchr', annotation uses '$achr'." >&2
+  exit 1
+fi
+
+leafcutter2 \
+  --juncfiles "$JUNCFILES" \
+  --outprefix "$OUTPREFIX" \
+  --rundir "$RUNDIR" \
+  --annot "$ANNOT" \
+  --genome "$GENOME" \
+  --maxintronlen "$MAXINTRONLEN" \
+  --minclureads "$MINCLUREADS" \
+  --minreads "$MINREADS" \
+  --minreadstd "$MINREADSTD" \
+  --mincluratio "$MINCLURATIO" \
+  --offset "$OFFSET" \
+  --keepleafcutter1
+
+# 3. An empty result means the filters removed every cluster.
+if [ "$(zcat "$RUNDIR/$OUTPREFIX.cluster_ratios.gz" | head -n2 | wc -l)" -lt 2 ]; then
+  echo "ERROR: no clusters survived filtering. Lower --min-clu-reads (currently $MINCLUREADS)" >&2
+  echo "  or --min-reads (currently $MINREADS), or check that the junction files are not empty." >&2
+  exit 1
+fi

code/script/molecular_phenotypes/calling/leafcutter_ds.R

@@ -0,0 +1,109 @@
+#!/usr/bin/env Rscript
+# leafcutter_ds.R from LeafCutter (https://github.com/davidaknowles/leafcutter,
+# scripts/leafcutter_ds.R; Apache License 2.0), used by the differential_splicing
+# workflow of splicing_calling_leafcutter2.ipynb.
+# One change from upstream: the count table is restricted to the samples present in
+# both it and the groups file, so a groups file naming a sample the table lacks drops
+# that sample instead of stopping the run.
+library(optparse)
+library(leafcutter)
+
+arguments <- parse_args(OptionParser(usage = "%prog [options] counts_file groups_file", description="LeafCutter differential splicing command line tool. Required inputs:\n <counts_file>: Intron usage counts file. Must be .txt or .txt.gz, output from clustering pipeline.\n <groups_file>: Two+K column file: 1. sample names (must match column names in counts_file), 2. groups (currently only two groups, i.e. pairwise, supported. Some samples in counts_file can be missing from this file, in which case they will not be included in the analysis. Additional columns can be used to specify confounders, e.g. batch/sex/age. Numeric columns will be treated as continuous, so use e.g. batch1, batch2, batch3 rather than 1, 2, 3 if you a categorical variable.",option_list=list(
+  make_option(c("-o","--output_prefix"), default = "leafcutter_ds", help="The prefix for the two output files, <prefix>_cluster_significance.txt (containing test status, log likelihood ratio, degree of freedom, and p-value for each cluster) and <prefix>_effect_sizes.txt (containing the effect sizes for each intron)  [default %default]"),
+  make_option(c("-s","--max_cluster_size"), default=Inf, help="Don't test clusters with more introns than this [default %default]"), 
+  make_option(c("-i","--min_samples_per_intron"), default=5, help="Ignore introns used (i.e. at least one supporting read) in fewer than n samples [default %default]") , 
+  make_option(c("-g","--min_samples_per_group"), default=3, help="Require this many samples in each group to have at least min_coverage reads [default %default]"), 
+  make_option(c("-c","--min_coverage"), default=20, help="Require min_samples_per_group samples in each group to have at least this many reads [default %default]"), 
+  make_option(c("-t","--timeout"), default=30, help="Maximum time (in seconds) allowed for a single optimization run [default %default]"),
+  make_option(c("-p","--num_threads"), default=1, help="Number of threads to use [default %default]"),
+  make_option(c("-e","--exon_file"), default=NULL, help="File defining known exons, example in data/gencode19_exons.txt.gz. Columns should be chr, start, end, strand, gene_name. Optional, only just to label the clusters."),
+  make_option(c("--init"), default="smart", help="One of 'smart' (default) or 'random'."), 
+  make_option(c("--seed"), default=12345, help="Random seed if using random initialization."))),
+  positional_arguments = 2)
+
+opt=arguments$opt
+counts_file=arguments$args[1]
+groups_file=arguments$args[2]
+
+cat("Loading counts from",counts_file,"\n")
+if (!file.exists(counts_file)) stop("File ",counts_file," does not exist")
+counts=read.table(counts_file, header=T, check.names = F)
+
+cat("Loading metadata from",groups_file,"\n")
+if (!file.exists(groups_file)) stop("File ",groups_file," does not exist")
+meta=read.table(groups_file, header=F, stringsAsFactors = F)
+colnames(meta)[1:2]=c("sample","group")
+
+common_cols <- intersect(colnames(counts), meta$sample)
+counts <- counts[, common_cols]
+meta <- meta[meta$sample %in% colnames(counts),]
+
+#counts=counts[,meta$sample]
+
+group_names=unique(meta$group) # keep order from groups_file unless numeric
+if (is.numeric(meta$group)) group_names=sort(group_names)
+meta$group=factor(meta$group, group_names)
+
+stopifnot(length(group_names)==2)
+
+cat("Encoding as",group_names[1],"=0,",group_names[2],"=1\n")
+numeric_x=as.numeric(meta$group)-1
+
+confounders=NULL
+if (ncol(meta)>2) {
+    confounders=meta[,3:ncol(meta),drop=F]
+    # scale continuous confounders
+    for (i in seq_len(ncol(confounders)))
+        if (is.numeric(confounders[,i]))
+            confounders[,i]=scale(confounders[,i])
+    # convert factors to one-of-K encoding
+    confounders=model.matrix( ~., data=confounders )
+    confounders=confounders[,2:ncol(confounders),drop=F] # remove intercept
+}
+
+minimum_group_size=min(sum(numeric_x==0),sum(numeric_x==1))
+if (minimum_group_size < opt$min_samples_per_intron)
+  stop("The number of samples in the smallest group is less than min_samples_per_intron, which means no clusters are testable. You can reduce min_samples_per_intron using the -i option, but note that we have only carefully checked the calibration of leafcutter p-values down to n=4 samples per group.")
+if (minimum_group_size < opt$min_samples_per_group)
+  stop("The number of samples in the smallest group is less than min_samples_per_group, which means no clusters are testable. You can reduce min_samples_per_intron using the -g option, but note that we have only carefully checked the calibration of leafcutter p-values down to n=4 samples per group.")
+
+require(doMC)
+registerDoMC(opt$num_threads)
+
+cat("Settings:\n")
+print(opt)
+
+cat("Running differential splicing analysis...\n")
+results <- differential_splicing(counts, numeric_x, confounders=confounders, max_cluster_size=opt$max_cluster_size, min_samples_per_intron=opt$min_samples_per_intron, min_samples_per_group=opt$min_samples_per_group, min_coverage=opt$min_coverage, timeout=opt$timeout, init=opt$init, seed=opt$seed) 
+
+cat("Saving results...\n")
+
+# Make cluster table
+cluster_table          = cluster_results_table(results)
+cluster_table$cluster  = add_chr(cluster_table$cluster)
+
+# Add gene names to clusters if an exon file is available
+if (!is.null(opt$exon_file)) {
+  cat("Loading exons from",opt$exon_file,"\n")
+  if (file.exists(opt$exon_file)) {
+     tryCatch( {
+          exons_table     = read.table(opt$exon_file, header=T, stringsAsFactors = F)
+          intron_meta     = get_intron_meta(rownames(counts))
+          exons_table$chr = add_chr(exons_table$chr)
+          intron_meta$chr = add_chr(intron_meta$chr)
+          clu_gene_map    = map_clusters_to_genes(intron_meta, exons_table)
+          cluster_table   = merge(cluster_table, clu_gene_map, by.x="cluster", by.y="clu", all.x=TRUE)
+     }, error=function(err) warning(as.character(err)) ) # ignore errors here
+  } else warning("File ",opt$exon_file," does not exist")
+} else cat("No exon_file provided.\n")
+write.table( cluster_table, paste0(opt$output_prefix,"_cluster_significance.txt"), quote=F, sep="\t", row.names = F)
+
+# Write effect size table
+effect_size_table                = leaf_cutter_effect_sizes(results)
+colnames(effect_size_table)[3:4] = group_names
+effect_size_table$intron = add_chr(effect_size_table$intron)
+write.table( effect_size_table, paste0(opt$output_prefix,"_effect_sizes.txt"), quote=F, col.names = T, row.names = F, sep="\t")
+
+# Save RData image
+# save.image(paste0(opt$output_prefix,"_cluster_significance.RData"));
+cat("All done, exiting\n")

code/script/molecular_phenotypes/calling/leafcutter_ds_merge.R

@@ -0,0 +1,80 @@
+#!/usr/bin/env Rscript
+# ============================================================
+# leafcutter_ds_merge.R
+# Join per-part leafcutter_ds.R results and recompute the FDR over all clusters,
+# for the differential_splicing_3 step of splicing_calling_leafcutter2.ipynb.
+#
+# Each part was tested separately, so its p.adjust column was corrected only within
+# that part. The p-values themselves do not depend on the split, so the column is
+# recomputed here with p.adjust(p, method = "fdr") over every cluster that has a
+# p-value, as leafcutter_ds.R does in an unsplit run. Clusters without a p-value
+# keep NA.
+#
+# Each <part>_cluster_significance.txt is read with its <part>_effect_sizes.txt
+# sibling. Other paths are ignored, so a step's whole output list can be passed.
+# Empty files (parts with no clusters) are skipped. Writes
+# <output-prefix>_cluster_significance.txt and <output-prefix>_effect_sizes.txt.
+# ============================================================
+
+suppressPackageStartupMessages(library(argparser))
+
+parser <- arg_parser("Join per-part leafcutter_ds.R results and recompute the FDR")
+parser <- add_argument(parser, "--output-prefix", type = "character", help = "<dir>/<prefix> for the joined tables")
+parser <- add_argument(parser, "--files", nargs = Inf, help = "per-part *_cluster_significance.txt files")
+argv <- parse_args(parser)
+if (is.na(argv$output_prefix)) stop("--output-prefix is required")
+if (all(is.na(argv$files))) stop("--files is required")
+
+die <- function(...) { message("ERROR: ", ...); quit(status = 1) }
+
+SIG <- "_cluster_significance.txt"
+EFF <- "_effect_sizes.txt"
+
+# Header and body lines of a table; NULL header for an empty file.
+read_table <- function(path) {
+  lines <- vroom::vroom_lines(path)
+  lines <- lines[nzchar(trimws(lines))]
+  if (!length(lines)) list(header = NULL, body = character(0))
+  else list(header = lines[[1L]], body = lines[-1L])
+}
+
+sig_files <- argv$files[endsWith(argv$files, SIG)]
+if (!length(sig_files)) die("no *", SIG, " files given")
+
+sig_hdr <- eff_hdr <- NULL
+rows <- list()
+eff_rows <- character(0)
+used <- 0L
+for (s in sig_files) {
+  e <- paste0(substr(s, 1L, nchar(s) - nchar(SIG)), EFF)
+  if (!file.exists(e)) die(e, " is missing")
+  sig <- read_table(s)
+  eff <- read_table(e)
+  if (is.null(sig$header)) next
+  if (!is.null(sig_hdr) && sig$header != sig_hdr) die(s, " has a different header")
+  if (!is.null(eff$header) && !is.null(eff_hdr) && eff$header != eff_hdr)
+    die(e, " has a different header (were the groups the same?)")
+  sig_hdr <- sig$header
+  if (!is.null(eff$header)) eff_hdr <- eff$header
+  rows <- c(rows, strsplit(sig$body, "\t", fixed = TRUE))
+  eff_rows <- c(eff_rows, eff$body)
+  used <- used + 1L
+}
+if (is.null(sig_hdr)) die("every part was empty; nothing was tested")
+
+cols <- strsplit(sig_hdr, "\t", fixed = TRUE)[[1]]
+ip <- match("p", cols)
+iq <- match("p.adjust", cols)
+p <- vapply(rows, `[[`, "", ip)
+tested <- !p %in% c("NA", "")
+q <- rep("NA", length(rows))
+adj <- p.adjust(as.numeric(p[tested]), method = "fdr")
+q[tested] <- as.character(adj)
+for (i in seq_along(rows)) rows[[i]][[iq]] <- q[[i]]
+
+vroom::vroom_write_lines(c(sig_hdr, vapply(rows, paste, "", collapse = "\t")),
+                         paste0(argv$output_prefix, SIG))
+vroom::vroom_write_lines(c(if (is.null(eff_hdr)) "" else eff_hdr, eff_rows),
+                         paste0(argv$output_prefix, EFF))
+cat("parts with clusters ", used, " of ", length(sig_files), "; clusters ", length(rows),
+    "; tested ", sum(tested), "; FDR < 0.05 ", sum(adj < 0.05), "\n", sep = "")

code/script/molecular_phenotypes/calling/leafcutter_ds_split.R

@@ -0,0 +1,79 @@
+#!/usr/bin/env Rscript
+# ============================================================
+# leafcutter_ds_split.R
+# Split a LeafCutter intron count table into one table per chromosome, for the
+# differential_splicing_1 step of splicing_calling_leafcutter2.ipynb.
+#
+# leafcutter_ds.R tests every cluster on its own, so the clusters can be tested in
+# separate jobs and the results joined afterwards. Only the multiple-testing
+# correction needs every cluster at once; leafcutter_ds_merge.R recomputes it.
+#
+# Rows go to parts 1-22 and X by their chromosome (a leading "chr" is ignored) and
+# everything else (Y, MT, unplaced contigs) to "other". All 24 parts are always
+# written, with the header, so the set of output files is fixed in advance; a part
+# with no clusters holds the header only.
+#
+# With --groups, the groups file is checked against the count table first, so a
+# naming mismatch stops the run here rather than in every per-part job.
+#
+# Writes <output-prefix>.part_<name>.counts.gz for each name in PARTS.
+# ============================================================
+
+suppressPackageStartupMessages(library(argparser))
+
+parser <- arg_parser("Split a LeafCutter intron count table by chromosome")
+parser <- add_argument(parser, "--counts", type = "character", help = "<prefix>_perind_numers.counts.gz")
+parser <- add_argument(parser, "--output-prefix", type = "character", help = "<dir>/<prefix> for the part files")
+parser <- add_argument(parser, "--groups", type = "character", default = NA,
+                       help = "optional groups file (sample, group, covariates) to check against the table")
+argv <- parse_args(parser)
+if (is.na(argv$counts)) stop("--counts is required")
+if (is.na(argv$output_prefix)) stop("--output-prefix is required")
+
+die <- function(...) { message("ERROR: ", ...); quit(status = 1) }
+
+PARTS <- c(as.character(1:22), "X", "other")
+
+part_of <- function(intron_id) {
+  chrom <- sub("^chr", "", sub(":.*$", "", intron_id))
+  ifelse(chrom %in% PARTS[-length(PARTS)], chrom, "other")
+}
+
+# Fail early on the mistakes leafcutter_ds.R would otherwise hit in every part.
+check_groups <- function(groups, header) {
+  columns <- strsplit(trimws(header), "\\s+")[[1]]
+  lines <- vroom::vroom_lines(groups)
+  rows <- strsplit(trimws(lines[nzchar(trimws(lines))]), "\\s+")
+  if (any(lengths(rows) < 2L))
+    die(groups, " needs at least two columns: sample and group")
+  sample <- vapply(rows, `[[`, "", 1L)
+  group <- vapply(rows, `[[`, "", 2L)
+  used <- sample %in% columns
+  if (!any(used))
+    die("no sample in ", groups, " matches a column of the count table. ",
+        "Column names look like: ", paste(head(sort(columns), 3L), collapse = " "))
+  labels <- sort(unique(group[used]))
+  if (length(labels) != 2L)
+    die(groups, " gives ", length(labels), " groups (", paste(labels, collapse = ", "),
+        ") among the matched samples; leafcutter_ds compares exactly two")
+  n <- table(factor(group[used], labels))
+  cat("groups: ", sum(used), " of ", length(rows), " samples matched; ",
+      paste0(labels, "=", n, collapse = ", "), "\n", sep = "")
+}
+
+lines <- vroom::vroom_lines(argv$counts)
+if (!length(lines) || !nzchar(trimws(lines[[1L]]))) die(argv$counts, " is empty")
+header <- lines[[1L]]
+body <- lines[-1L]
+body <- body[nzchar(trimws(body))]
+if (!is.na(argv$groups)) check_groups(argv$groups, header)
+
+d <- dirname(argv$output_prefix)
+if (nzchar(d)) dir.create(d, showWarnings = FALSE, recursive = TRUE)
+part <- factor(part_of(body), PARTS)
+for (p in PARTS)
+  vroom::vroom_write_lines(c(header, body[part == p]),
+                           paste0(argv$output_prefix, ".part_", p, ".counts.gz"))
+n <- table(part)
+cat("introns per part: ", paste0(PARTS, "=", n[PARTS], collapse = " "), "\n", sep = "")
+cat("total ", length(body), "\n", sep = "")