Add differential peak calling extension

assets/chrom_sizes_stub.sizes

@@ -0,0 +1 @@
+chr1	1

assets/differential_summary_stub.tsv

@@ -0,0 +1 @@
+

assets/ms_coeffs_stub.tsv

@@ -0,0 +1 @@
+sample_id	ms_coeff

assets/multiqc/chipbinner_summary_header.txt

@@ -0,0 +1,10 @@
+#id: 'chipbinner_summary'
+#parent_id: 'differential_analysis'
+#parent_name: 'Differential analysis'
+#parent_description: 'Differential peak calling outputs'
+#section_name: 'ChIPBinner summary'
+#description: 'Per-group ChIPBinner summary metrics.'
+#plot_type: 'table'
+#pconfig:
+#    id: 'chipbinner_summary_table'
+#    title: 'ChIPBinner summary'

assets/multiqc/differential_design_header.txt

@@ -0,0 +1,10 @@
+#id: 'differential_design'
+#parent_id: 'differential_analysis'
+#parent_name: 'Differential analysis'
+#parent_description: 'Differential peak calling outputs'
+#section_name: 'Differential analysis design'
+#description: 'Eligibility and design status for differential comparisons.'
+#plot_type: 'table'
+#pconfig:
+#    id: 'differential_design_table'
+#    title: 'Differential design'

assets/multiqc/differential_summary_header.txt

@@ -0,0 +1,10 @@
+#id: 'differential_summary'
+#parent_id: 'differential_analysis'
+#parent_name: 'Differential analysis'
+#parent_description: 'Differential peak calling outputs'
+#section_name: 'Differential analysis summary'
+#description: 'Summary of differential analysis across methods.'
+#plot_type: 'table'
+#pconfig:
+#    id: 'differential_summary_table'
+#    title: 'Differential summary'

assets/multiqc/span_summary_header.txt

@@ -0,0 +1,10 @@
+#id: 'span_summary'
+#parent_id: 'differential_analysis'
+#parent_name: 'Differential analysis'
+#parent_description: 'Differential peak calling outputs'
+#section_name: 'SPAN summary'
+#description: 'Per-group SPAN summary metrics.'
+#plot_type: 'table'
+#pconfig:
+#    id: 'span_summary_table'
+#    title: 'SPAN summary'

assets/multiqc_config.yml

@@ -64,6 +64,8 @@ report_section_order:
 
 custom_content:
   order:
+    - differential_design
+    - differential_summary
     - primary_peak_counts
     - consensus_peak_counts
     - primary_frip_score
@@ -94,6 +96,9 @@ section_comments:
 
 # Customise the module search patterns to speed up execution time
 sp:
+  custom_content:
+    fn: "*_mqc.tsv"
+
   cutadapt:
     fn: "*trimming_report.txt"
 

assets/multiqc_differential_config.yml

@@ -0,0 +1,18 @@
+report_comment: >
+  Differential analysis report generated by nf-core/cutandrun. The summary table
+  includes method-specific metrics, with dedicated ChIPBinner/SPAN sections for
+  detailed per-method summaries.
+
+run_modules:
+  - custom_content
+
+custom_content:
+  order:
+    - differential_design
+    - differential_summary
+    - chipbinner_summary
+    - span_summary
+
+sp:
+  custom_content:
+    fn: "*_mqc.tsv"

bin/annotate_regions.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+import argparse
+import csv
+import os
+import subprocess
+import tempfile
+
+
+def read_tsv(path):
+    with open(path, "r", newline="") as handle:
+        reader = csv.reader(handle, delimiter="\t")
+        rows = [row for row in reader]
+    return rows
+
+
+def is_header(row):
+    if len(row) < 3:
+        return True
+    try:
+        int(row[1])
+        return False
+    except Exception:
+        return True
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Annotate regions with nearest gene")
+    parser.add_argument("--regions", required=True)
+    parser.add_argument("--gene-bed", required=True)
+    parser.add_argument("--out", required=True)
+    args = parser.parse_args()
+
+    out_dir = os.path.dirname(args.out)
+    if out_dir:
+        os.makedirs(out_dir, exist_ok=True)
+
+    rows = read_tsv(args.regions)
+    if not rows:
+        with open(args.out, "w") as handle:
+            handle.write("")
+        return
+
+    header = None
+    data_rows = rows
+    if is_header(rows[0]):
+        header = rows[0]
+        data_rows = rows[1:]
+
+    annotation_cols = ["nearest_feature_id", "nearest_gene_name", "distance_to_feature"]
+
+    if not data_rows:
+        with open(args.out, "w") as handle:
+            if header:
+                handle.write("\t".join(header + annotation_cols) + "\n")
+        return
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        bed_path = os.path.join(tmpdir, "regions.bed")
+        with open(bed_path, "w") as handle:
+            for idx, row in enumerate(data_rows):
+                handle.write(f"{row[0]}\t{row[1]}\t{row[2]}\t{idx}\n")
+
+        cmd = [
+            "bedtools",
+            "closest",
+            "-d",
+            "-a",
+            bed_path,
+            "-b",
+            args.gene_bed,
+        ]
+        result = subprocess.run(cmd, capture_output=True, text=True, check=True)
+        annotations = {}
+        for line in result.stdout.strip().split("\n"):
+            if not line:
+                continue
+            parts = line.split("\t")
+            idx = int(parts[3])
+            a_cols = 4
+            distance = parts[-1] if parts else "NA"
+            b_cols = len(parts) - a_cols - 1
+            gene_id = "NA"
+            gene_name = "NA"
+            if b_cols >= 4:
+                gene_id = parts[a_cols + 3]
+            if b_cols >= 5:
+                gene_name = parts[a_cols + 4]
+            if gene_id in [".", ""]:
+                gene_id = "NA"
+            if gene_name in [".", ""]:
+                gene_name = "NA"
+            annotations[idx] = (gene_id, gene_name, distance)
+
+    with open(args.out, "w") as handle:
+        if header:
+            handle.write("\t".join(header + annotation_cols) + "\n")
+        for idx, row in enumerate(data_rows):
+            gene_id, gene_name, distance = annotations.get(idx, ("NA", "NA", "NA"))
+            handle.write("\t".join(row + [gene_id, gene_name, distance]) + "\n")
+
+
+if __name__ == "__main__":
+    main()

bin/chipbinner_rots.R

@@ -0,0 +1,89 @@
+#!/usr/bin/env Rscript
+
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) %% 2 != 0) {
+  stop("Arguments must be provided as --key value pairs")
+}
+arg_map <- list()
+for (i in seq(1, length(args), by = 2)) {
+  key <- sub("^--", "", args[i])
+  arg_map[[key]] <- args[i + 1]
+}
+
+matrix_path <- arg_map[["matrix"]]
+samplesheet_path <- arg_map[["samplesheet"]]
+treated_label <- arg_map[["treated"]]
+control_label <- arg_map[["control"]]
+bootstrap <- as.integer(arg_map[["bootstrap"]])
+k_value <- as.integer(arg_map[["k-value"]])
+out_path <- arg_map[["out"]]
+
+if (is.null(matrix_path) || is.null(samplesheet_path) || is.null(out_path)) {
+  stop("Missing required arguments")
+}
+
+suppressMessages(library(ROTS))
+
+mat <- read.table(matrix_path, header = TRUE, sep = "\t", check.names = FALSE, row.names = 1)
+mat <- as.matrix(mat)
+
+samples <- read.csv(samplesheet_path, header = TRUE, stringsAsFactors = FALSE)
+subset <- samples[samples$condition %in% c(treated_label, control_label), ]
+if (nrow(subset) == 0) {
+  stop("No samples match contrast labels")
+}
+
+sample_ids <- subset$sample_id
+missing <- setdiff(sample_ids, colnames(mat))
+if (length(missing) > 0) {
+  stop(paste("Missing samples in matrix:", paste(missing, collapse = ",")))
+}
+
+mat <- mat[, sample_ids, drop = FALSE]
+
+classes <- ifelse(subset$condition == treated_label, 1, 2)
+
+rot <- ROTS(mat, groups = classes, B = bootstrap, K = k_value)
+
+pvals <- NULL
+fdrs <- NULL
+
+if (!is.null(rot$pvalue)) {
+  pvals <- rot$pvalue
+} else if (!is.null(rot$pval)) {
+  pvals <- rot$pval
+}
+
+if (is.null(pvals)) {
+  if ("pvalue" %in% slotNames(rot)) {
+    pvals <- slot(rot, "pvalue")
+  }
+}
+
+if (is.null(pvals)) {
+  if (exists("pvalue", where = asNamespace("ROTS"), mode = "function")) {
+    pvals <- ROTS::pvalue(rot)
+  }
+}
+
+if (!is.null(rot$FDR)) {
+  fdrs <- rot$FDR
+} else if (!is.null(rot$fdr)) {
+  fdrs <- rot$fdr
+}
+
+if (is.null(fdrs) && !is.null(pvals)) {
+  fdrs <- p.adjust(pvals, method = "BH")
+}
+
+if (is.null(pvals)) {
+  stop("Unable to extract p-values from ROTS output")
+}
+
+out <- data.frame(
+  bin_id = rownames(mat),
+  pval = as.numeric(pvals),
+  FDR = as.numeric(fdrs)
+)
+
+write.table(out, file = out_path, sep = "\t", quote = FALSE, row.names = FALSE)