diff --git a/linkage_analysis/README.md b/linkage_analysis/README.md index ec8e574..b3b360e 100644 --- a/linkage_analysis/README.md +++ b/linkage_analysis/README.md @@ -1,109 +1,67 @@ -# Prerequisites -### Clone projects and restore environment: -```r -git clone https://github.com/stuart-lab/metapeak-analysis.git -cd metapeak-analysis/linkage_analysis -renv::restore() +# REMO analysis + +This repository contains code to reproduce results shown in *Regulatory element modules as universal features for single-cell chromatin analysis*: + +Chrysania Lim, Javen Tan Yih Ruay, and Tim Stuart, 2025. Regulatory element modules as universal features for single-cell chromatin analysis. *bioRxiv*. https://doi.org/10.64898/2025.12.10.692786 + +For code used to generate the REMO modules, please see this repository: https://github.com/stuart-lab/REMO + +For the R package containing REMO modules and associated metadata, please see this repository: https://github.com/stuart-lab/REMO.v1.GRCh38 + +## Expected inputs + +Recommended project layout: + +```text +resources/ + gencode.v32.basic.annotation.gtf.gz + GTEx_49tissues_release1.tsv.gz + multiome/ + _multiome_linkpeaks.rds +scripts/ + config/ + objects_and_links.csv + link_filtering_conds.csv +data/ + multiome_links/ + / + chr1.tsv + chr2.tsv + ... + all_chromosomes.tsv + combined_metadata/ + chr1.tsv + chr2.tsv + ... + filtered_links/ + / + pgl_metadata/pgls_metadata_all_chr.tsv + genes_per_peak/genes_per_peak_all_chr.tsv + gtex_validation/ ``` -For restoring renv refer to https://rstudio.github.io/renv/articles/renv.html#collaboration -### Prepare resources folder -See below for expected structure: -```bash -./linkage_analysis/resources -├── example -│   └── pgls_metadata_all_chr.tsv.gz (from supplementary data) -├── gencode.v32.basic.annotation.gtf.gz (from https://www.gencodegenes.org/human/release_32.html) -└── multiome (your Seurat multiome datasets[.rds]) -``` - -# Script setup -##### Modify slurm parameters as needed, e.g.(`cpus-per-task`, `mem-per-cpu`) -- `/scripts/01_link_peaks.sh` -- `/scripts/02_run_aggregate_links.sh` -- `/scripts/03_run_probability_analysis.sh` - -##### copy pwd to `PROJECT_ROOT` of .sh files -```bash -pwd -$ /metapeak-analysis/linkage_analysis/ -``` - -# Reproduce results -In `/scripts/03_run_probability_analysis.sh`, change `SIG_PGL_FPATH` variable to `"${PROJECT_ROOT}/resources/example/pgl_metadata_all_chr.tsv"` - -Run the following script -```bash -./scripts/03_run_probability_analysis.sh -``` - -# Run full analysis workflow on user-supplied data -### 1. Replace examples files in `/resources/multiome/` with your intended Seurat multiome objects (`.rds`) - -### 2. Modify the following variables in `01_link_peaks.sh` according to each seurat multiome object -- `OUTPUT_DIR` (links for each chromosome are stored here, in `/data/multiome_links/dataset_name`) -- `SEURAT_RDS` (filepath for Seurat multiome object, in `/resources/multiome`) - -E.g. -```bash -OUTPUT_DIR="${PROJECT_ROOT}/data/multiome_links/jejunum" -SEURAT_RDS="${PROJECT_ROOT}/resources/multiome/jejunum_multiome_linkpeaks.rds" -``` - -### 3. Run the follow script: -```bash -./01_link_peaks.sh -``` -Repeat for each Seurat multiome object. +`objects_and_links.csv` is used by step 2 for both aggregation and GTEx/SGL validation. It must contain a `links_dir` column pointing to each dataset's `data/multiome_links/` directory. A `dataset` column is optional; if omitted, the dataset label is inferred from the link directory name. -### 4. Create the following config files in `/scripts/configs` -##### `objects_and_links.csv` (modify example table of seurat obj filepaths and their mutiome_link dirs) -```csv -seurat_obj_fpath,links_dir -/resources/multiome/jejunum_multiome_linkpeaks.rds,/data/multiome_links/jejunum/ -/resources/multiome/bile_duct_multiome_linkpeaks.rds,/data/multiome_links/bile_duct/ -``` - -##### `link_filtering_conds.csv` (specify p-value and score thresholds to filter peak-gene links for downstream analysis) -```csv -condition,pvalue,score -sig_pvalue,0.05,0 -hsig_pvalue_high_coeff,0.01,0.1 -``` +The GTEx/SGL fine-mapping file used by the validation part of step 2 should be placed at `resources/GTEx_49tissues_release1.tsv.gz` by default. It is expected to contain at least the columns `gene`, `chromosome`, `end`, `tissue`, and `pip`. -### 5. Run the following script: -```bash -./scripts/02_run_aggregate_links.sh -``` +## Workflow -### 6. Modify `/scripts/03_run_probability_analysis.sh` -Specify the `SIG_PGL_FPATH` variable as the path to the combined pgl_metadata for the condition of choice (e.g. `sig_pvalue`), example below: -```bash -${PROJECT_ROOT}/data/filtered_links/sig_pvalue/pgl_metadata/pgls_metadata_all_chr.tsv -``` - -### 7. Run the following script: -```bash -./scripts/03_run_probability_analysis.sh -``` +1. `01_link_peaks.R` / `01_link_peaks.sh` + - Run Signac V2-compatible peak-gene linking per dataset. + - The supplied GENCODE GTF is attached to the peak assay. + - Genes with non-unique symbols are excluded before candidate link discovery. + - Per-chromosome link TSVs, `all_chromosomes.tsv`, and a linked Seurat RDS are written. -# Results +2. `02_aggregate_links.R`, `02a_validate_links_against_gtex.R` / `02_aggregate_and_validate_links.sh` + - Run aggregation first, then GTEx/SGL validation using the aggregated outputs. + - `02_aggregate_links.R` combines per-chromosome link TSVs across datasets, adds REMO metadata, creates per-condition filtered PGL metadata, and builds proximal/distal/all gene lists per peak. + - `02a_validate_links_against_gtex.R` validates link metrics against GTEx/SGL fine-mapping positives. + - Uses `scripts/config/objects_and_links.csv` for per-dataset validation. + - Uses `data/combined_metadata/chr*.tsv`, written by the aggregation step, for all-dataset validation. + - Writes per-dataset outputs under `data/gtex_validation/per_dataset//`. + - Writes all-dataset outputs under `data/gtex_validation/all_datasets/`. + - For each dataset and for all aggregated links, writes ROC plots, AUC tables, GTEx-labelled link tables, and diagnostics. + - Metrics evaluated: `score`, `zscore`, `exp(-distance / 200 kb)`, and `score * exp(-distance / 200 kb)`. The redundant `1 / distance` metric is not used. -Output structure: -```bash -./linkage_analysis/data -├── combined_metadata -├── filtered_links -├── multiome_links -└── probability_analysis -``` - -A table of conditional probabilities can be found at: -``` -data/probability_analysis/conditional_prob_unlinked_singleCRE.tsv -``` - -The per-chromosome odds ratio for the likelihood that CRE pairs linked to the same gene are within the same REMO module can be found at: -``` -data/probability_analysis/logistic_regression_results/glm_results_all.tsv -``` \ No newline at end of file +3. `03_get_conditional_probs.R`, `03a_get_glm_coeffs.R`, `03b_build_glm_summary.R` / `03_run_probability_analysis.sh` + - Run conditional probability and GLM analyses on a selected filtered PGL metadata file. diff --git a/linkage_analysis/scripts/01_link_peaks.R b/linkage_analysis/scripts/01_link_peaks.R index 13fa775..d019d22 100644 --- a/linkage_analysis/scripts/01_link_peaks.R +++ b/linkage_analysis/scripts/01_link_peaks.R @@ -4,13 +4,19 @@ suppressPackageStartupMessages({ library(argparse) library(Signac) library(Seurat) + library(SeuratObject) library(GenomicRanges) + library(GenomeInfoDb) + library(IRanges) + library(S4Vectors) + library(InteractionSet) library(Matrix) library(data.table) library(rtracklayer) library(future) library(future.apply) library(pbapply) + library(lifecycle) library(REMO.v1.GRCh38) }) @@ -78,6 +84,11 @@ parser$add_argument( help = "Absolute correlation cutoff for peaks before background z-scoring [default: 0]" ) +parser$add_argument( + "--tss_match_strategy", default = "transcript_nearest", + help = "TSS matching strategy: transcript_nearest or most5 [default: transcript_nearest]" +) + parser$add_argument( "--workers", type = "integer", default = 1, help = "Number of parallel workers (future) [default: 1]" @@ -92,14 +103,32 @@ peak_assay <- args$peak_assay expression_assay <- args$expression_assay peak_layer <- args$peak_layer expression_layer <- args$expression_layer +tss_match_strategy <- match.arg( + args$tss_match_strategy, + choices = c("transcript_nearest", "most5") +) +link_key <- "linkpeaks" dir.create(output_dir, showWarnings = FALSE, recursive = TRUE) +combined_links_tsv <- file.path(output_dir, "all_chromosomes.tsv") +seurat_rds_base <- basename(seurat_rds) +seurat_rds_stem <- sub("\\.[^.]+$", "", seurat_rds_base) +linked_seurat_rds <- file.path( + output_dir, + paste0( + seurat_rds_stem, + "_with_", link_key, "_", tss_match_strategy, "_zscore.rds" + ) +) + message("Input Seurat RDS: ", seurat_rds) message("GENCODE GTF: ", gencode_gtf) message("Output directory: ", output_dir) message("Peak assay: ", peak_assay, " | layer: ", peak_layer) message("Expression assay: ", expression_assay, " | layer: ", expression_layer) +message("TSS strategy: ", tss_match_strategy) +message("Z-score mode: enabled; background matching will be run") # -------------------------- # # Optional renv activation @@ -128,315 +157,1456 @@ if (workers > 1) { } # -------------------------- # -# Helper: CollapseDuplicateColumns +# Signac v2-compatible link helpers # -------------------------- # -#' Collapse duplicate columns of a dgCMatrix using max aggregation. -#' -#' @param mat A dgCMatrix with possibly duplicated column names. -#' @return A dgCMatrix with unique column names, where each entry is the max -#' across all original duplicate columns. -CollapseDuplicateColumns <- function(mat) { - if (!inherits(mat, "dgCMatrix")) { - stop("Input must be of class dgCMatrix as sparse matrix operations are used") - } - columns <- colnames(mat) - if (is.null(columns)) { - stop("Matrix must have column names.") - } - - u_colnames <- unique(columns) - - # Coerce to triplet form for efficient grouped aggregation - T <- as(mat, "dgTMatrix") - - # If there are no non-zero entries, we still want to - # collapse to one column per unique name, but everything is zero. - if (length(T@x) == 0L) { - out <- Matrix::sparseMatrix( - i = integer(0), - j = integer(0), - x = numeric(0), - dims = c(nrow(mat), length(u_colnames)) - ) - rownames(out) <- rownames(mat) - colnames(out) <- u_colnames +DistanceToTSS <- getFromNamespace("DistanceToTSS", "Signac") +Extend <- getFromNamespace("Extend", "Signac") +corSparse <- getFromNamespace("corSparse", "Signac") +SparseSpearmanCor <- getFromNamespace("SparseSpearmanCor", "Signac") + +make_empty_links <- function() { + InteractionSet::GInteractions( + GenomicRanges::GRanges(), + GenomicRanges::GRanges() + ) +} + + +# -------------------------- # +# Signac v2 / GTF annotation helpers +# -------------------------- # + +get_layer_compat <- function(assay, layer) { + out <- tryCatch( + LayerData(object = assay, layer = layer), + error = function(e) NULL + ) + if (!is.null(out)) { return(out) } - # Map each non-zero entry's column to its unique-name index - group <- match(columns[T@j + 1L], u_colnames) + if (layer %in% slotNames(assay)) { + out <- slot(assay, layer) + if (!is.null(out) && length(dim(out)) == 2) { + return(out) + } + } - df <- data.frame( - i = T@i + 1L, # 1-based row index - j = group, # 1-based unique column index - x = T@x # values + stop( + "Could not retrieve layer/slot '", layer, "' from assay of class: ", + paste(class(assay), collapse = ", ") ) +} - agg <- stats::aggregate(x ~ i + j, data = df, FUN = max) +parse_peak_names_compat <- function(regions, sep = c("-", "-")) { + gr <- tryCatch( + GenomicRanges::GRanges(regions), + error = function(e) NULL + ) + if (!is.null(gr)) { + return(gr) + } - out <- Matrix::sparseMatrix( - i = agg$i, - j = agg$j, - x = agg$x, - dims = c(nrow(mat), length(u_colnames)) + gr <- tryCatch( + Signac:::StringToGRanges(regions = regions, sep = sep), + error = function(e) NULL + ) + if (!is.null(gr)) { + return(gr) + } + + stop( + "Could not parse peak names as genomic ranges. Example names: ", + paste(utils::head(regions), collapse = ", ") ) - colnames(out) <- u_colnames - rownames(out) <- rownames(mat) - out } -# -------------------------- # -# Helper: LinkPeaksCustom -# -------------------------- # +get_peak_ranges_compat <- function(assay, peak.layer = "counts") { + gr <- tryCatch( + granges(x = assay), + error = function(e) NULL + ) + if (!is.null(gr) && length(gr) > 0) { + return(gr) + } -LinkPeaksCustom <- function( - object, - peak.assay, - expression.assay, - peak.slot = "counts", - expression.slot = "data", - method = "pearson", - gene.coords = NULL, - distance = 5e+05, - min.distance = NULL, - min.cells = 10, - genes.use = NULL, - n_sample = 200, - pvalue_cutoff = 0.05, - score_cutoff = 0.05, - gene.id = FALSE, - verbose = TRUE -) { + if ("ranges" %in% slotNames(assay)) { + gr <- slot(assay, "ranges") + if (inherits(gr, "GRanges") && length(gr) > 0) { + return(gr) + } + } - if (!inherits(x = object[[peak.assay]], what = "ChromatinAssay")) { - stop("The requested assay is not a ChromatinAssay") + mat <- get_layer_compat(assay = assay, layer = peak.layer) + parse_peak_names_compat(rownames(mat)) +} + +get_feature_metadata_compat <- function(assay) { + md <- tryCatch( + assay[[]], + error = function(e) NULL + ) + if (!is.null(md) && nrow(md) > 0) { + return(md) } - if (!is.null(x = min.distance)) { - if (!is.numeric(x = min.distance)) stop("min.distance must be numeric") - if (min.distance <= 0) min.distance <- NULL + if ("meta.features" %in% slotNames(assay)) { + md <- slot(assay, "meta.features") + if (!is.null(md) && nrow(md) > 0) { + return(md) + } } - features.match <- c("GC.percent", "count", "sequence.length") - cor_method <- switch( - method, - pearson = corSparse, - spearman = SparseSpearmanCor, - stop("method must be 'pearson' or 'spearman'") + NULL +} + +get_fragments_compat <- function(assay) { + frags <- tryCatch( + Fragments(assay), + error = function(e) NULL ) + if (!is.null(frags)) { + return(frags) + } - if (is.null(x = gene.coords)) { - annot <- Annotation(object = object[[peak.assay]]) - if (is.null(x = annot)) stop("Gene annotations not found") - gene.coords <- CollapseToLongestTranscript(annot) + if ("fragments" %in% slotNames(assay)) { + frags <- slot(assay, "fragments") + if (!is.null(frags)) { + return(frags) + } } - meta.features <- GetAssayData( - object = object, assay = peak.assay, layer = "meta.features" + NULL +} + +prepare_gencode_annotation <- function(gtf, rel_levels = NULL) { + gtf_ann <- gtf + + if (!any(grepl("^chr", as.character(seqnames(gtf_ann))))) { + seqlevels(gtf_ann) <- paste0("chr", seqlevels(gtf_ann)) + } + + gtf_ann <- GenomeInfoDb::keepStandardChromosomes( + gtf_ann, + pruning.mode = "coarse" ) - if (!(all(c("GC.percent", "sequence.length") %in% colnames(meta.features)))) { - stop("Run RegionsStats before calling this function.") + + if (!is.null(rel_levels)) { + gtf_ann <- gtf_ann[as.character(seqnames(gtf_ann)) %in% rel_levels] } - if (!("count" %in% colnames(meta.features))) { - data.use <- GetAssayData(object = object[[peak.assay]], layer = "counts") - hvf.info <- FindTopFeatures(object = data.use, verbose = FALSE) - hvf.info <- hvf.info[rownames(meta.features), , drop = FALSE] - meta.features <- cbind(meta.features, hvf.info) + GenomeInfoDb::genome(gtf_ann) <- "hg38" + + if (!"gene_biotype" %in% colnames(mcols(gtf_ann)) && + "gene_type" %in% colnames(mcols(gtf_ann))) { + gtf_ann$gene_biotype <- gtf_ann$gene_type + } + + if (!"tx_id" %in% colnames(mcols(gtf_ann)) && + "transcript_id" %in% colnames(mcols(gtf_ann))) { + gtf_ann$tx_id <- gtf_ann$transcript_id } - peak.data <- GetAssayData( - object = object, assay = peak.assay, layer = peak.slot + required.cols <- c("gene_id", "gene_name", "gene_biotype", "type") + required.ok <- all(required.cols %in% colnames(mcols(gtf_ann))) && + any(c("tx_id", "transcript_id") %in% colnames(mcols(gtf_ann))) + + if (!required.ok) { + stop( + "GENCODE annotation is missing columns required by Signac. Present columns: ", + paste(colnames(mcols(gtf_ann)), collapse = ", ") + ) + } + + keep.cols <- intersect( + c("tx_id", "transcript_id", "gene_id", "gene_name", "gene_biotype", "type"), + colnames(mcols(gtf_ann)) ) + mcols(gtf_ann) <- mcols(gtf_ann)[, keep.cols, drop = FALSE] + + gtf_ann +} + +replace_peak_annotation_with_gtf <- function(object, peak.assay, annotation, peak.layer = "counts") { + object <- SeuratObject::UpdateSeuratObject(object) + + if (!peak.assay %in% names(object@assays)) { + stop("Peak assay '", peak.assay, "' not found in Seurat object.") + } + + old.assay <- object[[peak.assay]] + + if (inherits(old.assay, "GRangesAssay")) { + Annotation(old.assay) <- annotation + object[[peak.assay]] <- old.assay + } else if (inherits(old.assay, "ChromatinAssay")) { + old.counts <- get_layer_compat(old.assay, peak.layer) + old.data <- tryCatch(get_layer_compat(old.assay, "data"), error = function(e) NULL) + old.ranges <- get_peak_ranges_compat(old.assay, peak.layer = peak.layer) + old.feature.md <- get_feature_metadata_compat(old.assay) + old.fragments <- get_fragments_compat(old.assay) + + if (length(old.ranges) != nrow(old.counts)) { + stop("Number of peak ranges does not match number of peak count rows.") + } + names(old.ranges) <- rownames(old.counts) + + new.assay <- tryCatch( + as.GRangesAssay(old.assay), + error = function(e) NULL + ) + + if (is.null(new.assay)) { + if (is.null(old.data)) { + new.assay <- CreateGRangesAssay( + counts = old.counts, + ranges = old.ranges, + annotation = annotation, + fragments = old.fragments, + validate.fragments = FALSE + ) + } else { + new.assay <- CreateGRangesAssay( + counts = old.counts, + data = old.data, + ranges = old.ranges, + annotation = annotation, + fragments = old.fragments, + validate.fragments = FALSE + ) + } + } + + Annotation(new.assay) <- annotation + + if (!is.null(old.feature.md) && nrow(old.feature.md) == nrow(old.counts)) { + rownames(old.feature.md) <- rownames(old.counts) + new.assay <- tryCatch( + AddMetaData(object = new.assay, metadata = old.feature.md), + error = function(e) { + warning("Could not copy peak-level metadata to converted assay: ", conditionMessage(e)) + new.assay + } + ) + } + + object[[peak.assay]] <- new.assay + } else { + stop( + "Peak assay '", peak.assay, "' must be a Signac ChromatinAssay or GRangesAssay. Current class: ", + paste(class(old.assay), collapse = ", ") + ) + } + + if (!inherits(object[[peak.assay]], "GRangesAssay")) { + stop("Peak assay is still not a GRangesAssay after compatibility conversion.") + } + if (is.null(Annotation(object[[peak.assay]])) || length(Annotation(object[[peak.assay]])) == 0) { + stop("GENCODE annotation was not attached to peak assay.") + } + + object +} + +get_nonunique_gene_symbols <- function(gencode) { + gencode.genes <- gencode[ + gencode$type == "gene" & + as.character(seqnames(gencode)) != "chrY" + ] + + gene.dt <- data.table::data.table( + seqnames = as.character(seqnames(gencode.genes)), + start = start(gencode.genes), + end = end(gencode.genes), + strand = as.character(strand(gencode.genes)), + gene_id = gencode.genes$gene_id, + gene_name = gencode.genes$gene_name + )[!is.na(gene_name) & gene_name != ""] + + multi.chr.genes <- gene.dt[ + , .N, by = .(gene_name, seqnames) + ][, .N, by = gene_name][N > 1L, gene_name] + + multi.gene.same.chr <- gene.dt[ + , .N, by = .(gene_name, seqnames, gene_id, start, end, strand) + ][, .N, by = .(gene_name, seqnames)][N > 1L, gene_name] + + unique(c(multi.chr.genes, multi.gene.same.chr)) +} + +# Harmonised version for local use by FindCandidateLinks(). +CollapseToLongestTranscript <- function(ranges) { + range.df <- data.table::as.data.table(x = as.data.frame(x = ranges)) + + if (!("gene_id" %in% colnames(x = range.df))) { + stop("Input ranges must contain gene_id") + } + + if (!("gene_name" %in% colnames(x = range.df))) { + range.df[, gene_name := gene_id] + } - if (!(expression.slot %in% Layers(object[[expression.assay]]))) { - stop("Requested expression layer not found in assay '", expression.assay, "'") + if (!("gene_biotype" %in% colnames(x = range.df))) { + range.df[, gene_biotype := if ("gene_type" %in% colnames(x = range.df)) { + gene_type + } else { + NA_character_ + }] } - expression.data <- GetAssayData( - object = object, assay = expression.assay, layer = expression.slot + range.df$strand <- as.character(x = range.df$strand) + range.df$strand <- ifelse( + test = range.df$strand == "*", + yes = "+", + no = range.df$strand ) - peakcounts <- Matrix::rowSums(peak.data > 0) - genecounts <- Matrix::rowSums(expression.data > 0) + collapsed <- range.df[ + , list( + seqnames = unique(seqnames)[[1]], + start = min(start), + end = max(end), + strand = strand[[1]], + gene_biotype = gene_biotype[[1]], + gene_name = gene_name[[1]] + ), + by = "gene_id" + ] + + collapsed$gene_name <- make.unique(names = collapsed$gene_name) + + gene.ranges <- GenomicRanges::makeGRangesFromDataFrame( + df = collapsed, + keep.extra.columns = TRUE + ) - peaks.keep <- peakcounts > min.cells - genes.keep <- genecounts > min.cells - peak.data <- peak.data[peaks.keep, ] + return(gene.ranges) +} - if (!is.null(genes.use)) { - genes.keep <- intersect(names(genes.keep[genes.keep]), genes.use) +#' Find candidate peak-gene links +#' +#' Identify candidate peak-gene pairs within a specified distance from a gene +#' transcription start site (TSS). Candidates can be matched either to the +#' collapsed most-5-prime gene model or to the nearest transcript TSS for each +#' peak-gene pair. +#' +#' @param object A Seurat object. +#' @param peak.assay Name of assay containing peak information. +#' @param expression.assay Name of assay containing gene expression information. +#' @param peak.layer Name of layer to pull chromatin data from. +#' @param expression.layer Name of layer to pull expression data from. +#' @param gene.coords A [GenomicRanges::GRanges] object containing gene or +#' transcript coordinates. If `NULL`, gene annotations are extracted from the +#' peak assay. +#' @param distance Maximum distance from a TSS for a peak-gene pair to be +#' considered a candidate link. +#' @param min.distance Optional minimum distance from a TSS. Candidate links +#' closer than this distance are excluded. +#' @param min.cells Minimum number of cells positive for a peak or gene for it +#' to be retained. +#' @param genes.use Optional vector of genes to test. If `NULL`, genes are +#' selected from the expression assay after `min.cells` filtering. +#' @param gene.id Set to `TRUE` if genes in the expression assay are named by +#' gene ID rather than gene name. +#' @param verbose Display messages. +#' @param tss.match.strategy Strategy used to assign candidate TSSs. Use +#' `"most5"` to use one collapsed gene-level TSS, or `"transcript_nearest"` to +#' choose the nearest transcript TSS for each peak-gene pair. +#' +#' @return A list containing objects used downstream by [LinkPeaks()]: +#' `candidate.matrix`, `candidate.links`, `peaks`, `gene.coords`, +#' `candidate.link.table`, `peak.data`, and `expression.data`. +#' +#' @importFrom GenomicRanges GRanges granges seqnames start end width strand resize findOverlaps +#' @importFrom IRanges IRanges +#' @importFrom S4Vectors mcols queryHits subjectHits +#' @importFrom InteractionSet GInteractions +#' @importFrom Matrix rowSums sparseMatrix drop0 summary +#' @importFrom SeuratObject LayerData Layers +#' @importFrom stats ave +#' +#' @export +#' @concept links +FindCandidateLinks <- function( + object, + peak.assay, + expression.assay, + peak.layer = "counts", + expression.layer = "data", + gene.coords = NULL, + distance = 5e5, + min.distance = NULL, + min.cells = 10, + genes.use = NULL, + gene.id = FALSE, + verbose = TRUE, + tss.match.strategy = c("transcript_nearest", "most5") +) { + tss.match.strategy <- match.arg(tss.match.strategy) + gene.key <- if (gene.id) "gene_id" else "gene_name" + + .peak_to_tss_distance <- function(peaks, tss, peak_index, tss_index) { + peak.starts <- start(peaks)[peak_index] + peak.ends <- end(peaks)[peak_index] + tss.pos <- start(tss)[tss_index] + + ifelse( + tss.pos >= peak.starts & tss.pos <= peak.ends, + 0, + pmin(abs(tss.pos - peak.starts), abs(tss.pos - peak.ends)) + ) } - ## ---- CASE 1: no genes left for this chromosome ---- - if (length(genes.keep) == 0L) { - if (verbose) message("No genes remain after filtering; skipping.") - Links(object[[peak.assay]]) <- GenomicRanges::GRanges() - return(object) + .add_closest_gene_info <- function(df) { + df$closest_tss <- as.integer( + df$tss_distance == ave(df$tss_distance, df$peak_index, FUN = min) + ) + df$closest_gene <- df$closest_tss + df$closest_tss_id <- NA_character_ + + if ("transcript_id" %in% colnames(df)) { + by.peak <- split(seq_len(nrow(df)), df$peak_index) + + for (idx in by.peak) { + closest.idx <- idx[df$closest_tss[idx] == 1L] + df$closest_tss_id[idx] <- paste( + unique(as.character(df$transcript_id[closest.idx])), + collapse = ";" + ) + } + } + + df + } + + .add_peak_gene_overlap_info <- function(df, peaks, gene.body.coords, gene.key) { + df$peak_overlap_gene <- 0L + df$frac_peak_in_gene <- 0 + + if (nrow(df) == 0 || length(gene.body.coords) == 0) { + return(df) + } + + gene.body.values <- as.character( + mcols(gene.body.coords)[[gene.key]] + ) + + body.index <- match(as.character(df$gene), gene.body.values) + ok <- !is.na(body.index) + + if (any(ok)) { + peak.ranges <- peaks[df$peak_index[ok]] + body.ranges <- gene.body.coords[body.index[ok]] + + same.seq <- as.character(seqnames(peak.ranges)) == + as.character(seqnames(body.ranges)) + + overlap.start <- pmax( + start(peak.ranges), + start(body.ranges) + ) + overlap.end <- pmin( + end(peak.ranges), + end(body.ranges) + ) + + overlap.width <- ifelse( + same.seq, + pmax(0L, overlap.end - overlap.start + 1L), + 0L + ) + + df$peak_overlap_gene[ok] <- as.integer(overlap.width > 0L) + df$frac_peak_in_gene[ok] <- overlap.width / width(peak.ranges) + } + + df + } + + if (!inherits(object[[peak.assay]], "GRangesAssay")) { + stop("The requested peak assay is not a GRangesAssay") + } + + if (!is.null(min.distance)) { + if (!is.numeric(min.distance)) stop("min.distance should be numeric") + if (min.distance <= 0) min.distance <- NULL + if (!is.null(min.distance) && min.distance >= distance) { + stop("min.distance should be smaller than distance") + } + } + + gene.coords.input <- gene.coords + if (is.null(gene.coords.input)) { + gene.coords.input <- Annotation(object[[peak.assay]]) + if (is.null(gene.coords.input)) stop("Gene annotations not found") + } + + if (!("gene_id" %in% colnames(mcols(gene.coords.input)))) { + stop("gene.coords / annotation must contain 'gene_id'") + } + + if (!("gene_name" %in% colnames(mcols(gene.coords.input)))) { + if (gene.id) { + mcols(gene.coords.input)$gene_name <- + mcols(gene.coords.input)$gene_id + } else { + stop("gene.coords / annotation must contain 'gene_name' when gene.id = FALSE") + } + } + + if (!(peak.layer %in% Layers(object[[peak.assay]]))) { + stop("Requested peak layer not found") + } + + if (!(expression.layer %in% Layers(object[[expression.assay]]))) { + stop("Requested expression layer not found") } + peak.data <- LayerData(object[[peak.assay]], layer = peak.layer) + expression.data <- LayerData(object[[expression.assay]], layer = expression.layer) + + peaks.keep <- rowSums(peak.data > 0) > min.cells + genes.keep <- names(which(rowSums(expression.data > 0) > min.cells)) + + if (!is.null(genes.use)) { + genes.keep <- intersect(genes.keep, genes.use) + } + + peak.data <- peak.data[peaks.keep, , drop = FALSE] expression.data <- expression.data[genes.keep, , drop = FALSE] + + if (nrow(peak.data) == 0) stop("No peaks pass min.cells filtering") + if (nrow(expression.data) == 0) stop("No genes pass min.cells / genes.use filtering") + genes <- rownames(expression.data) + peaks <- granges(object[[peak.assay]])[peaks.keep] + peak.names <- rownames(peak.data) + if (is.null(peak.names)) peak.names <- as.character(peaks) + + gene.body.coords <- CollapseToLongestTranscript(gene.coords.input) if (gene.id) { - gene.coords.use <- gene.coords[gene.coords$gene_id %in% genes, ] - gene.coords.use$gene_name <- gene.coords.use$gene_id + gene.body.coords <- gene.body.coords[gene.body.coords$gene_id %in% genes] + gene.body.coords$gene_name <- gene.body.coords$gene_id } else { - gene.coords.use <- gene.coords[gene.coords$gene_name %in% genes, ] + gene.body.coords <- gene.body.coords[gene.body.coords$gene_name %in% genes] } - ## ---- CASE 2: no coordinates for these genes ---- - if (length(gene.coords.use) == 0L) { - if (verbose) message("No gene coordinates found; skipping chromosome.") - Links(object[[peak.assay]]) <- GenomicRanges::GRanges() - return(object) + gene.body.values <- as.character( + mcols(gene.body.coords)[[gene.key]] + ) + keep.body.genes <- genes[genes %in% gene.body.values] + gene.body.coords <- gene.body.coords[match(keep.body.genes, gene.body.values)] + + if (tss.match.strategy == "most5") { + gene.coords.use <- gene.body.coords + gene.values <- mcols(gene.coords.use)[[gene.key]] + + keep.genes <- genes[genes %in% gene.values] + if (length(keep.genes) == 0) { + stop("No expressed genes matched gene coordinates") + } + + gene.coords.use <- gene.coords.use[match(keep.genes, gene.values)] + gene.values <- mcols(gene.coords.use)[[gene.key]] + + candidate.matrix <- DistanceToTSS( + peaks = peaks, + genes = gene.coords.use, + distance = distance + ) + + rownames(candidate.matrix) <- peak.names + colnames(candidate.matrix) <- gene.values + + if (!is.null(min.distance)) { + inner.matrix <- DistanceToTSS( + peaks = peaks, + genes = gene.coords.use, + distance = min.distance + ) + + rownames(inner.matrix) <- rownames(candidate.matrix) + colnames(inner.matrix) <- colnames(candidate.matrix) + + candidate.matrix <- candidate.matrix - inner.matrix + candidate.matrix@x[candidate.matrix@x < 0] <- 0 + candidate.matrix <- drop0(candidate.matrix) + } + + if (sum(candidate.matrix) == 0) { + stop("No candidate peak-gene links found") + } + + candidate.tss <- resize( + gene.coords.use, + width = 1, + fix = "start" + ) + mcols(candidate.tss) <- NULL + + candidate.link.table <- as.data.frame(summary(candidate.matrix)) + colnames(candidate.link.table) <- c("peak_index", "gene_index", "x") + + candidate.link.table$peak <- rownames(candidate.matrix)[candidate.link.table$peak_index] + candidate.link.table$gene <- gene.values[candidate.link.table$gene_index] + candidate.link.table$transcript_id <- NA_character_ + + candidate.link.table$tss_distance <- .peak_to_tss_distance( + peaks = peaks, + tss = candidate.tss, + peak_index = candidate.link.table$peak_index, + tss_index = candidate.link.table$gene_index + ) + + candidate.link.table <- .add_closest_gene_info(candidate.link.table) + candidate.link.table <- .add_peak_gene_overlap_info( + df = candidate.link.table, + peaks = peaks, + gene.body.coords = gene.body.coords, + gene.key = gene.key + ) + + candidate.links <- GInteractions( + peaks[candidate.link.table$peak_index], + candidate.tss[candidate.link.table$gene_index] + ) + + candidate.links$peak <- candidate.link.table$peak + candidate.links$gene <- candidate.link.table$gene + candidate.links$transcript_id <- candidate.link.table$transcript_id + candidate.links$tss_distance <- candidate.link.table$tss_distance + candidate.links$peak_overlap_gene <- candidate.link.table$peak_overlap_gene + candidate.links$frac_peak_in_gene <- candidate.link.table$frac_peak_in_gene + candidate.links$closest_gene <- candidate.link.table$closest_gene + candidate.links$closest_tss_id <- candidate.link.table$closest_tss_id } - if (length(gene.coords.use) < nrow(expression.data)) { - message("Found coordinates for ", length(gene.coords.use), - " of ", nrow(expression.data), " genes.") + if (tss.match.strategy == "transcript_nearest") { + transcript.coords.use <- gene.coords.input + + if ("type" %in% colnames(mcols(transcript.coords.use))) { + tx.rows <- transcript.coords.use$type %in% c("transcript", "mRNA") + if (any(tx.rows, na.rm = TRUE)) { + transcript.coords.use <- transcript.coords.use[tx.rows] + } + } + + if (gene.id) { + transcript.coords.use <- transcript.coords.use[ + transcript.coords.use$gene_id %in% genes + ] + transcript.coords.use$gene_name <- transcript.coords.use$gene_id + } else { + transcript.coords.use <- transcript.coords.use[ + transcript.coords.use$gene_name %in% genes + ] + } + + if (length(transcript.coords.use) == 0) { + stop("No transcript/gene coordinates found") + } + + transcript.gene.values <- mcols(transcript.coords.use)[[gene.key]] + gene.values <- genes[genes %in% unique(transcript.gene.values)] + + if (length(gene.values) == 0) { + stop("No expressed genes matched transcript coordinates") + } + + gene.coords.use <- transcript.coords.use[match(gene.values, transcript.gene.values)] + + transcript.tss <- resize( + transcript.coords.use, + width = 1, + fix = "start" + ) + + outer.hits <- findOverlaps( + query = peaks, + subject = Extend( + transcript.tss, + upstream = distance, + downstream = distance + ), + type = "any" + ) + + if (length(outer.hits) == 0) { + stop("No candidate peak-transcript links found") + } + + outer.df <- data.frame( + peak_index = queryHits(outer.hits), + transcript_index = subjectHits(outer.hits) + ) + + outer.df$gene <- transcript.gene.values[outer.df$transcript_index] + outer.df <- outer.df[outer.df$gene %in% gene.values, , drop = FALSE] + + if (nrow(outer.df) == 0) { + stop("No candidate peak-transcript links matched expressed genes") + } + + if (!is.null(min.distance)) { + inner.hits <- findOverlaps( + query = peaks, + subject = Extend( + transcript.tss, + upstream = min.distance, + downstream = min.distance + ), + type = "any" + ) + + if (length(inner.hits) > 0) { + inner.df <- data.frame( + peak_index = queryHits(inner.hits), + transcript_index = subjectHits(inner.hits) + ) + + inner.df$gene <- transcript.gene.values[inner.df$transcript_index] + inner.df <- inner.df[inner.df$gene %in% gene.values, , drop = FALSE] + + outer.df <- outer.df[ + !(paste(outer.df$peak_index, outer.df$gene, sep = "\r") %in% + paste(inner.df$peak_index, inner.df$gene, sep = "\r")), + , + drop = FALSE + ] + } + } + + if (nrow(outer.df) == 0) { + stop("No candidate links remain after min.distance filtering") + } + + outer.df$gene_index <- match(outer.df$gene, gene.values) + + outer.df$tss_distance <- .peak_to_tss_distance( + peaks = peaks, + tss = transcript.tss, + peak_index = outer.df$peak_index, + tss_index = outer.df$transcript_index + ) + + outer.df$tss_seqname <- as.character( + seqnames(transcript.tss) + )[outer.df$transcript_index] + + outer.df$tss_position <- start(transcript.tss)[ + outer.df$transcript_index + ] + + outer.df$tss_strand <- as.character( + strand(transcript.tss) + )[outer.df$transcript_index] + + outer.df$transcript_id <- if ( + "transcript_id" %in% colnames(mcols(transcript.coords.use)) + ) { + as.character(transcript.coords.use$transcript_id[outer.df$transcript_index]) + } else { + paste0("transcript_row_", outer.df$transcript_index) + } + + outer.df <- outer.df[order( + outer.df$peak_index, + outer.df$gene_index, + outer.df$tss_distance, + outer.df$transcript_index + ), , drop = FALSE] + + outer.df <- outer.df[ + !duplicated(paste(outer.df$peak_index, outer.df$gene_index, sep = "\r")), + , + drop = FALSE + ] + + outer.df$peak <- peak.names[outer.df$peak_index] + outer.df <- .add_closest_gene_info(outer.df) + outer.df <- .add_peak_gene_overlap_info( + df = outer.df, + peaks = peaks, + gene.body.coords = gene.body.coords, + gene.key = gene.key + ) + + candidate.matrix <- sparseMatrix( + i = outer.df$peak_index, + j = outer.df$gene_index, + x = 1, + dims = c(length(peaks), length(gene.values)) + ) + + rownames(candidate.matrix) <- peak.names + colnames(candidate.matrix) <- gene.values + + candidate.tss <- GRanges( + seqnames = outer.df$tss_seqname, + ranges = IRanges( + start = outer.df$tss_position, + end = outer.df$tss_position + ), + strand = outer.df$tss_strand + ) + mcols(candidate.tss) <- NULL + + candidate.links <- GInteractions( + peaks[outer.df$peak_index], + candidate.tss + ) + + candidate.links$peak <- outer.df$peak + candidate.links$gene <- outer.df$gene + candidate.links$transcript_id <- outer.df$transcript_id + candidate.links$tss_distance <- outer.df$tss_distance + candidate.links$peak_overlap_gene <- outer.df$peak_overlap_gene + candidate.links$frac_peak_in_gene <- outer.df$frac_peak_in_gene + candidate.links$closest_gene <- outer.df$closest_gene + candidate.links$closest_tss_id <- outer.df$closest_tss_id + + candidate.link.table <- outer.df } - peaks <- granges(object[[peak.assay]])[peaks.keep] + if (verbose) { + message( + "Identified ", length(candidate.links), + " candidate peak-gene links across ", ncol(candidate.matrix), + " genes and ", sum(rowSums(candidate.matrix) > 0), + " peaks using strategy '", tss.match.strategy, "'" + ) + } + + candidate.link.table <- candidate.link.table[, c("peak", "gene"), drop = FALSE] - peak_distance_matrix <- Signac:::DistanceToTSS( + # Return only objects used downstream by LinkPeaks(). + list( + candidate.matrix = candidate.matrix, + candidate.links = candidate.links, peaks = peaks, - genes = gene.coords.use, - distance = distance + gene.coords = gene.coords.use, + candidate.link.table = candidate.link.table, + peak.data = peak.data, + expression.data = expression.data ) - peak_distance_matrix <- CollapseDuplicateColumns(peak_distance_matrix) +} - if (!is.null(min.distance)) { - peak_distance_matrix_min <- Signac:::DistanceToTSS( - peaks = peaks, - genes = gene.coords.use, - distance = min.distance +#' Link peaks to genes +#' +#' Find peaks that are correlated with the expression of nearby genes. For each +#' gene, this function computes the correlation coefficient between gene +#' expression and accessibility of candidate peaks within a specified distance +#' from the gene TSS. By default, links are retained using the observed +#' correlation coefficient only. If `calculate.zscore = TRUE`, the function also +#' computes an expected correlation coefficient for each peak using background +#' peaks matched on GC content, accessibility, and sequence length, and uses this +#' background distribution to compute a z-score and p-value. +#' +#' This function was inspired by the method originally described by SHARE-seq +#' (Sai Ma et al. 2020, Cell). Please consider citing the original SHARE-seq +#' work if using this function: +#' [doi: 10.1016/j.cell.2020.09.056](https://pubmed.ncbi.nlm.nih.gov/33098772/) +#' +#' @param object A Seurat object. +#' @param peak.assay Name of assay containing peak information. +#' @param expression.assay Name of assay containing gene expression information. +#' @param peak.layer Name of layer to pull chromatin data from. +#' @param expression.layer Name of layer to pull expression data from. +#' @param method Correlation method to use. One of `"pearson"` or `"spearman"`. +#' @param key Key to use when storing link information in the assay. +#' @param gene.coords A [GenomicRanges::GRanges] object containing gene or +#' transcript coordinates. If `NULL`, gene annotations are extracted from the +#' peak assay. +#' @param distance Maximum distance from a TSS for peaks to include as candidate +#' links. +#' @param min.distance Optional minimum distance from a TSS. Candidate links +#' closer than this distance are excluded. +#' @param min.cells Minimum number of cells positive for the peak and gene +#' needed to include them in the analysis. +#' @param genes.use Optional vector of genes to test. If `NULL`, genes are +#' selected from the expression assay after `min.cells` filtering. +#' @param n_sample Number of background peaks to sample when computing the null +#' distribution. Only used when `calculate.zscore = TRUE`. +#' @param pvalue_cutoff Maximum p-value for retaining a link. Links with a +#' p-value greater than or equal to this value are removed from the output. Only +#' used when `calculate.zscore = TRUE`. +#' @param score_cutoff Minimum absolute correlation coefficient for a link to be +#' retained. +#' @param calculate.zscore Compute background-matched z-scores and p-values. If +#' `FALSE` (default), background peak matching is skipped and links are retained +#' using `score_cutoff` only. +#' @param gene.id Set to `TRUE` if genes in the expression assay are named by +#' gene ID rather than gene name. +#' @param verbose Display messages. +#' @param tss.match.strategy Strategy used to assign candidate TSSs. Use +#' `"most5"` to use one collapsed gene-level TSS, or `"transcript_nearest"` to +#' choose the nearest transcript TSS for each peak-gene pair. +#' @param peak.slot Deprecated; use `peak.layer`. +#' @param expression.slot Deprecated; use `expression.layer`. +#' +#' @return Returns a Seurat object with results added to the links slot in the +#' peak assay, stored under `key`. The results are stored as an +#' [InteractionSet::GInteractions] object accessible via [Links()]. The metadata +#' stored on each link includes: +#' * `peak`: peak identifier +#' * `gene`: linked gene identifier +#' * `transcript_id`: selected transcript identifier, when available +#' * `tss_distance`: distance from the peak to the selected TSS +#' * `peak_overlap_gene`: binary indicator for any peak overlap with the linked +#' gene body +#' * `frac_peak_in_gene`: fraction of the peak width overlapping the linked gene +#' body +#' * `closest_gene`: binary indicator for whether this gene is closest to the +#' peak among candidate genes +#' * `closest_tss_id`: transcript IDs for the closest TSSs to the peak +#' * `score`: observed correlation coefficient +#' * `zscore`: z-score of the observed correlation coefficient; `NA` when +#' `calculate.zscore = FALSE` +#' * `pvalue`: p-value associated with the z-score; `NA` when +#' `calculate.zscore = FALSE` +#' +#' @importFrom SeuratObject LayerData Layers as.sparse +#' @importFrom Matrix sparseMatrix rowSums drop0 +#' @importMethodsFrom Matrix t +#' @importFrom GenomicRanges seqnames +#' @importFrom S4Vectors mcols DataFrame +#' @importFrom future.apply future_lapply +#' @importFrom future nbrOfWorkers +#' @importFrom pbapply pblapply +#' @importFrom lifecycle is_present deprecated deprecate_warn +#' @importFrom stats pnorm sd +#' +#' @export +#' @concept links +LinkPeaks <- function( + object, + peak.assay, + expression.assay, + peak.layer = "counts", + expression.layer = "data", + method = "pearson", + key = "linkpeaks", + gene.coords = NULL, + distance = 5e+05, + min.distance = NULL, + min.cells = 10, + genes.use = NULL, + n_sample = 200, + pvalue_cutoff = 0.05, + score_cutoff = 0.05, + calculate.zscore = FALSE, + gene.id = FALSE, + verbose = TRUE, + tss.match.strategy = c("transcript_nearest", "most5"), + peak.slot = deprecated(), + expression.slot = deprecated() +) { + tss.match.strategy <- match.arg(tss.match.strategy) + + if (!inherits(x = object[[peak.assay]], what = "GRangesAssay")) { + stop("The requested assay is not a GRangesAssay") + } + + if (!is.null(x = min.distance)) { + if (!is.numeric(x = min.distance)) { + stop("min.distance should be a numeric value") + } + if (min.distance < 0) { + warning("Requested a negative min.distance value, setting min.distance to zero") + min.distance <- NULL + } else if (min.distance == 0) { + min.distance <- NULL + } + } + + if (is_present(arg = expression.slot)) { + deprecate_warn( + when = "1.16.0", + what = "LinkPeaks(expression.slot)", + with = "LinkPeaks(expression.layer)" + ) + expression.layer <- expression.slot + } + + if (is_present(arg = peak.slot)) { + deprecate_warn( + when = "1.16.0", + what = "LinkPeaks(peak.slot)", + with = "LinkPeaks(peak.layer)" ) - peak_distance_matrix_min <- CollapseDuplicateColumns(peak_distance_matrix_min) - peak_distance_matrix <- peak_distance_matrix - peak_distance_matrix_min + peak.layer <- peak.slot + } + + features.match <- c("GC.percent", "count", "sequence.length") + + if (method == "pearson") { + cor_method <- corSparse + } else if (method == "spearman") { + cor_method <- SparseSpearmanCor + } else { + stop("method can be one of 'pearson' or 'spearman'.") + } + + meta.features <- NULL + if (isTRUE(calculate.zscore)) { + meta.features <- object[[peak.assay]][[]] + if (!(all(c("GC.percent", "sequence.length") %in% colnames(x = meta.features)))) { + stop( + "DNA sequence information for each peak has not been computed.\n", + "Run RegionStats before calling this function." + ) + } + + if (!("count" %in% colnames(x = meta.features))) { + hvf.info <- FindTopFeatures( + object = LayerData(object = object[[peak.assay]], layer = peak.layer), + verbose = FALSE + ) + hvf.info <- hvf.info[rownames(meta.features), , drop = FALSE] + meta.features <- cbind(meta.features, hvf.info) + } + } + + candidates <- tryCatch( + FindCandidateLinks( + object = object, + peak.assay = peak.assay, + expression.assay = expression.assay, + peak.layer = peak.layer, + expression.layer = expression.layer, + gene.coords = gene.coords, + distance = distance, + min.distance = min.distance, + min.cells = min.cells, + genes.use = genes.use, + gene.id = gene.id, + verbose = verbose, + tss.match.strategy = tss.match.strategy + ), + error = function(e) { + if (verbose) message("Candidate-finding skipped: ", conditionMessage(e)) + NULL + } + ) + + if (is.null(candidates)) { + Links(object = object[[peak.assay]], key = key) <- make_empty_links() + return(object) + } + + peak_distance_matrix <- candidates$candidate.matrix + candidate.links <- candidates$candidate.links + candidate.link.table <- candidates$candidate.link.table + peak.data <- candidates$peak.data + expression.data <- candidates$expression.data + gene.coords.use <- candidates$gene.coords + genes.use <- colnames(x = peak_distance_matrix) + all.peaks <- rownames(x = peak.data) + all.peak.chroms <- if (isTRUE(calculate.zscore)) { + as.character(seqnames(candidates$peaks)) + } else { + NULL } - ## ---- CASE 3: no peaks within distance window ---- if (sum(peak_distance_matrix) == 0) { - if (verbose) message("No peaks within distance; skipping chromosome.") - Links(object[[peak.assay]]) <- GenomicRanges::GRanges() + if (verbose) message("No peaks fall within distance threshold.") + Links(object = object[[peak.assay]], key = key) <- make_empty_links() return(object) } if (verbose) { - message("Testing ", nrow(expression.data), " genes and ", - sum(Matrix::rowSums(peak_distance_matrix) > 0), " peaks") + message( + "Testing ", nrow(x = expression.data), + " genes and ", sum(rowSums(x = peak_distance_matrix) > 0), + " peaks" + ) + } + + peak.data <- t(x = peak.data) + + if (nbrOfWorkers() > 1) { + mylapply <- future_lapply + } else { + mylapply <- if (verbose) pblapply else lapply } - genes.use <- colnames(peak_distance_matrix) - all.peaks <- rownames(peak.data) + res <- mylapply( + X = seq_along(along.with = genes.use), + FUN = function(i) { + peak.use <- as.logical(x = peak_distance_matrix[, genes.use[[i]]]) + gene.expression <- t(x = expression.data[genes.use[[i]], , drop = FALSE]) + + if (sum(peak.use) < 2) { + return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) + } + + peak.access <- peak.data[, peak.use, drop = FALSE] + if (inherits(x = peak.access, what = "IterableMatrix")) { + peak.access <- as.sparse(x = peak.access) + } + if (inherits(x = gene.expression, what = "IterableMatrix")) { + gene.expression <- as.sparse(x = gene.expression) + } + + coef.result <- cor_method( + X = peak.access, + Y = gene.expression + ) + rownames(x = coef.result) <- colnames(x = peak.access) + coef.result <- coef.result[ + abs(x = coef.result) > score_cutoff, + , + drop = FALSE + ] + + if (nrow(x = coef.result) == 0) { + return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) + } + + peaks.test <- rownames(x = coef.result) + + if (!isTRUE(calculate.zscore)) { + coef.vec <- as.numeric(coef.result[, 1]) + names(x = coef.vec) <- peaks.test + return(list( + "gene" = rep(i, length(x = coef.vec)), + "coef" = coef.vec, + "zscore" = rep(NA_real_, length(x = coef.vec)) + )) + } + + gene.chrom <- as.character(x = seqnames(x = gene.coords.use[i])) + + # Original LinkPeaks samples background peaks from chromosomes different + # from the gene. Use actual peak chromosomes instead of assuming a + # chr-start-end peak-name format. + trans.peaks <- all.peaks[all.peak.chroms != gene.chrom] + meta.use <- meta.features[trans.peaks, , drop = FALSE] + pk.use <- meta.features[peaks.test, , drop = FALSE] + + bg.peaks <- lapply( + X = seq_len(length.out = nrow(x = pk.use)), + FUN = function(x) { + MatchRegionStats( + meta.feature = meta.use, + query.feature = pk.use[x, , drop = FALSE], + features.match = features.match, + n = n_sample, + verbose = FALSE + ) + } + ) + + unique.bg <- unique(x = unlist(x = bg.peaks)) + if (length(unique.bg) == 0) { + return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) + } + + bg.access <- peak.data[, unique.bg, drop = FALSE] + if (inherits(x = bg.access, what = "IterableMatrix")) { + bg.access <- as.sparse(x = bg.access) + } + + bg.coef <- cor_method( + X = bg.access, + Y = gene.expression + ) + rownames(x = bg.coef) <- unique.bg + + zscores <- vector(mode = "numeric", length = length(x = peaks.test)) + for (j in seq_along(along.with = peaks.test)) { + coef.use <- bg.coef[bg.peaks[[j]], , drop = FALSE] + bg.sd <- sd(x = coef.use) + if (bg.sd == 0 || !is.finite(x = bg.sd)) { + zscores[[j]] <- 0 + } else { + zscores[[j]] <- (coef.result[j] - mean(x = coef.use)) / bg.sd + } + } + + names(x = coef.result) <- peaks.test + names(x = zscores) <- peaks.test + + pval.vec <- 2 * pnorm(q = -abs(x = zscores)) + links.keep <- pval.vec < pvalue_cutoff + + if (sum(x = links.keep) == 0) { + return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) + } + + list( + "gene" = rep(i, sum(links.keep)), + "coef" = coef.result[links.keep, , drop = TRUE], + "zscore" = zscores[links.keep] + ) + } + ) + + gene.vec <- do.call(what = c, args = lapply(X = res, FUN = `[[`, 1)) + coef.vec <- do.call(what = c, args = lapply(X = res, FUN = `[[`, 2)) + zscore.vec <- do.call(what = c, args = lapply(X = res, FUN = `[[`, 3)) + + if (length(x = coef.vec) == 0) { + if (verbose) message("No links pass score_cutoff") + Links(object = object[[peak.assay]], key = key) <- make_empty_links() + return(object) + } - peak.data <- Matrix::t(peak.data) + sig.df <- data.frame( + peak = names(x = coef.vec), + gene_index = as.integer(gene.vec), + gene = genes.use[as.integer(gene.vec)], + score = as.numeric(coef.vec), + zscore = if (isTRUE(calculate.zscore)) { + as.numeric(zscore.vec) + } else { + rep(NA_real_, length(x = coef.vec)) + }, + stringsAsFactors = FALSE + ) - coef.vec <- gene.vec <- zscore.vec <- c() + if (isTRUE(calculate.zscore)) { + sig.df$pvalue <- 2 * pnorm(q = -abs(x = sig.df$zscore)) + sig.df <- sig.df[sig.df$pvalue < pvalue_cutoff, , drop = FALSE] - if (future::nbrOfWorkers() > 1) { - mylapply <- future.apply::future_lapply + if (nrow(sig.df) == 0) { + if (verbose) message("No significant links after p-value filtering") + Links(object = object[[peak.assay]], key = key) <- make_empty_links() + return(object) + } } else { - mylapply <- if (verbose) pbapply::pblapply else lapply + sig.df$pvalue <- NA_real_ } - res <- mylapply(seq_along(genes.use), function(i) { - peak.use <- as.logical(peak_distance_matrix[, genes.use[[i]]]) - gene.expression <- Matrix::t(expression.data[genes.use[[i]], , drop = FALSE]) - gene.chrom <- as.character(seqnames(gene.coords.use[i])) + candidate.key <- paste(candidate.link.table$peak, candidate.link.table$gene, sep = "\r") + sig.key <- paste(sig.df$peak, sig.df$gene, sep = "\r") + link.idx <- match(sig.key, candidate.key) - if (verbose) message(sum(peak.use), " peaks found for ", genes.use[[i]]) - if (sum(peak.use) < 2) - return(list(gene=NULL, coef=NULL, zscore=NULL)) + if (anyNA(link.idx)) { + missing <- unique(sig.key[is.na(link.idx)]) + stop( + "Internal error: significant links could not be matched back to candidate links. ", + "First missing key: ", missing[[1]] + ) + } - peak.access <- peak.data[, peak.use, drop = FALSE] - coef.result <- cor_method(X = peak.access, Y = gene.expression) - rownames(coef.result) <- colnames(peak.access) - coef.result <- coef.result[abs(coef.result) > score_cutoff, , drop=FALSE] + links <- candidate.links[link.idx] + links$peak <- sig.df$peak + links$gene <- sig.df$gene - if (nrow(coef.result) == 0) - return(list(gene=NULL, coef=NULL, zscore=NULL)) + if (!("transcript_id" %in% colnames(mcols(links)))) { + links$transcript_id <- rep(NA_character_, length(links)) + } + if (!("closest_tss_id" %in% colnames(mcols(links)))) { + links$closest_tss_id <- rep(NA_character_, length(links)) + } + if (!("peak_overlap_gene" %in% colnames(mcols(links)))) { + links$peak_overlap_gene <- rep(NA_integer_, length(links)) + } + if (!("frac_peak_in_gene" %in% colnames(mcols(links)))) { + links$frac_peak_in_gene <- rep(NA_real_, length(links)) + } + if (!("closest_gene" %in% colnames(mcols(links)))) { + if ("closest_tss" %in% colnames(mcols(links))) { + links$closest_gene <- as.integer(links$closest_tss) + } else { + links$closest_gene <- rep(NA_integer_, length(links)) + } + } - peaks.test <- rownames(coef.result) - trans.peaks <- all.peaks[!grepl(paste0("^", gene.chrom, "-"), all.peaks)] - meta.use <- meta.features[trans.peaks, ] - pk.use <- meta.features[peaks.test, ] + links$score <- sig.df$score + links$zscore <- sig.df$zscore + links$pvalue <- sig.df$pvalue - bg.peaks <- lapply(seq_len(nrow(pk.use)), function(x) { - MatchRegionStats( - meta.feature = meta.use, - query.feature = pk.use[x, , drop = FALSE], - features.match = features.match, - n = n_sample, - verbose = FALSE + .link_mcol <- function(x, col) { + value <- mcols(x)[[col]] + + if (length(value) != length(x)) { + stop( + "Internal error: links$", col, " has length ", length(value), + " but length(links) is ", length(x), "." ) - }) + } - bg.access <- peak.data[, unlist(bg.peaks), drop = FALSE] - bg.coef <- cor_method(X = bg.access, Y = gene.expression) - rownames(bg.coef) <- colnames(bg.access) + value + } - zscores <- numeric(length(peaks.test)) - for (j in seq_along(peaks.test)) { - coef.use <- bg.coef[((j-1)*n_sample+1):(j*n_sample), ] - zscores[j] <- (coef.result[j] - mean(coef.use)) / sd(coef.use) - } + # Keep only fields required downstream for export and saved link objects. + mcols(links) <- DataFrame( + peak = as.character(.link_mcol(links, "peak")), + gene = as.character(.link_mcol(links, "gene")), + transcript_id = as.character(.link_mcol(links, "transcript_id")), + tss_distance = as.numeric(.link_mcol(links, "tss_distance")), + peak_overlap_gene = as.integer(.link_mcol(links, "peak_overlap_gene")), + frac_peak_in_gene = as.numeric(.link_mcol(links, "frac_peak_in_gene")), + closest_gene = as.integer(.link_mcol(links, "closest_gene")), + closest_tss_id = as.character(.link_mcol(links, "closest_tss_id")), + score = as.numeric(.link_mcol(links, "score")), + zscore = as.numeric(.link_mcol(links, "zscore")), + pvalue = as.numeric(.link_mcol(links, "pvalue")) + ) + + Links(object = object[[peak.assay]], key = key) <- links + object +} + + +# -------------------------- # +# Helper: link extraction and TSV output +# -------------------------- # + +is_links_object <- function(x) { + inherits(x, "GInteractions") || inherits(x, "GenomicRanges") || inherits(x, "GRanges") +} + +normalize_links_object <- function(links, key = NULL, empty.ok = TRUE) { + if (is.null(links)) { + if (empty.ok) return(make_empty_links()) + stop("Link object is NULL.") + } + + if (is_links_object(links)) { + return(links) + } - names(coef.result) <- peaks.test - names(zscores) <- peaks.test + if (is.list(links)) { + if (!is.null(key) && key %in% names(links)) { + return(normalize_links_object(links[[key]], key = NULL, empty.ok = empty.ok)) + } - list( - gene = rep(i, length(coef.result)), - coef = coef.result, - zscore = zscores + normalized <- lapply( + links, + function(x) { + tryCatch( + normalize_links_object(x, key = key, empty.ok = FALSE), + error = function(e) NULL + ) + } ) - }) + normalized <- Filter(Negate(is.null), normalized) - ## Combine outputs - gene.vec <- unlist(lapply(res, `[[`, 1)) - coef.vec <- unlist(lapply(res, `[[`, 2)) - zscore.vec <- unlist(lapply(res, `[[`, 3)) + if (length(normalized) == 0) { + if (empty.ok) return(make_empty_links()) + stop("Could not find a GInteractions/GRanges link object inside list.") + } - ## ---- CASE 4: no significant links ---- - if (length(coef.vec) == 0) { - if (verbose) message("No significant links; skipping chromosome.") - Links(object[[peak.assay]]) <- GenomicRanges::GRanges() - return(object) + nonempty <- Filter(function(x) length(x) > 0, normalized) + if (length(nonempty) == 0) { + return(make_empty_links()) + } + if (length(nonempty) == 1) { + return(nonempty[[1]]) + } + return(do.call(c, nonempty)) } - peak.key <- seq_along(unique(names(coef.vec))) - names(peak.key) <- unique(names(coef.vec)) + stop("Unsupported link object type: ", paste(class(links), collapse = ", ")) +} - coef.matrix <- Matrix::sparseMatrix( - i = gene.vec, - j = peak.key[names(coef.vec)], - x = coef.vec, - dims = c(length(genes.use), max(peak.key)) +get_links_normalized <- function(assay, key = NULL) { + links <- tryCatch( + Signac::Links(object = assay, key = key), + error = function(e) { + tryCatch( + Signac::Links(object = assay), + error = function(e2) NULL + ) + } ) - rownames(coef.matrix) <- genes.use - colnames(coef.matrix) <- names(peak.key) - links <- Signac:::LinksToGRanges(coef.matrix, gene.coords = gene.coords.use) + normalize_links_object(links, key = key) +} + +combine_links_objects <- function(link.list) { + normalized <- lapply(link.list, normalize_links_object) + nonempty <- Filter(function(x) length(x) > 0, normalized) - z.matrix <- Matrix::sparseMatrix( - i = gene.vec, - j = peak.key[names(zscore.vec)], - x = zscore.vec, - dims = c(length(genes.use), max(peak.key)) + if (length(nonempty) == 0) { + return(make_empty_links()) + } + + do.call(c, nonempty) +} + +links_to_output_df <- function(links, key = NULL) { + wanted <- c( + "peak", + "gene", + "transcript_id", + "tss_distance", + "distance", + "peak_overlap_gene", + "frac_peak_in_gene", + "closest_gene", + "closest_tss_id", + "score", + "zscore", + "pvalue" ) - rownames(z.matrix) <- genes.use - colnames(z.matrix) <- names(peak.key) - z.lnk <- Signac:::LinksToGRanges(z.matrix, gene.coords = gene.coords.use) + links <- normalize_links_object(links, key = key) + + if (length(links) == 0) { + out <- data.frame( + peak = character(), + gene = character(), + transcript_id = character(), + tss_distance = numeric(), + distance = numeric(), + peak_overlap_gene = integer(), + frac_peak_in_gene = numeric(), + closest_gene = integer(), + closest_tss_id = character(), + score = numeric(), + zscore = numeric(), + pvalue = numeric() + ) + return(out) + } + + meta <- as.data.frame(S4Vectors::mcols(links)) + if (!("distance" %in% colnames(meta)) && "tss_distance" %in% colnames(meta)) { + meta$distance <- meta$tss_distance + } + for (nm in wanted) { + if (!(nm %in% colnames(meta))) meta[[nm]] <- NA + } + meta[, wanted, drop = FALSE] +} + +write_chromosome_tsv <- function(links, chromosome, output.dir, key = NULL) { + links_df <- links_to_output_df(links, key = key) + links_df$chromosome <- chromosome + links_df <- links_df[ + , + c( + "chromosome", + "peak", + "gene", + "transcript_id", + "tss_distance", + "distance", + "peak_overlap_gene", + "frac_peak_in_gene", + "closest_gene", + "closest_tss_id", + "score", + "zscore", + "pvalue" + ), + drop = FALSE + ] - links$zscore <- z.lnk$score - links$pvalue <- pnorm(-abs(links$zscore)) - links <- links[links$pvalue < pvalue_cutoff] + data.table::fwrite( + links_df, + file = file.path(output.dir, sprintf("%s.tsv", chromosome)), + sep = "\t", + quote = FALSE + ) - Links(object[[peak.assay]]) <- links - return(object) + links_df } # -------------------------- # @@ -449,22 +1619,35 @@ data <- readRDS(seurat_rds) rel_levels <- seqlevels(REMO.v1.GRCh38) message("Loading GENCODE annotations (GTF)...") -gencode <- import(gencode_gtf) +gencode <- prepare_gencode_annotation( + gtf = import(gencode_gtf), + rel_levels = rel_levels +) + +message("Replacing peak assay annotation with supplied GENCODE GTF...") +data <- replace_peak_annotation_with_gtf( + object = data, + peak.assay = peak_assay, + annotation = gencode, + peak.layer = peak_layer +) + +nonuniq_symbols <- get_nonunique_gene_symbols(gencode) +message("Excluded non-unique gene symbols: ", length(nonuniq_symbols)) + gencode_transcripts <- gencode[ gencode$type == "transcript" & - seqnames(gencode) %in% rel_levels -] -# remove genes with the same name across multiple chr -dt_gc <- as.data.table(gencode_transcripts)[, .(gene_name, seqnames = as.character(seqnames))] -multi_chr_genes <- dt_gc[, uniqueN(chr) > 1, gene_name] -gencode_transcripts <- gencode_transcripts[ - !(gencode_transcripts$gene_name %in% multi_chr_genes) + seqnames(gencode) %in% rel_levels & + !(gencode$gene_name %in% nonuniq_symbols) ] message("Chromosomes to process: ", paste(rel_levels, collapse = ", ")) DefaultAssay(data) <- peak_assay +all.outputs <- list() +all.link.objects <- list() + # -------------------------- # # Per-chromosome linking # -------------------------- # @@ -475,43 +1658,53 @@ for (level in rel_levels) { message("Processing chromosome: ", level) gene_coords_chr <- gencode_transcripts[seqnames(gencode_transcripts) == level] - gene_set_chr <- gene_coords_chr$gene_name + gene_set_chr <- unique(as.character(gene_coords_chr$gene_name)) + gene_set_chr <- gene_set_chr[!is.na(gene_set_chr) & gene_set_chr != ""] message(" Genes on this chromosome: ", length(gene_set_chr)) - data <- LinkPeaksCustom( - object = data, - peak.assay = peak_assay, - expression.assay = expression_assay, - peak.slot = peak_layer, - expression.slot = expression_layer, - gene.coords = gene_coords_chr, - genes.use = gene_set_chr, - distance = args$distance, - min.distance = args$min_distance, - min.cells = args$min_cells, - n_sample = args$n_sample, - pvalue_cutoff = args$pvalue_cutoff, - score_cutoff = args$score_cutoff, - verbose = TRUE + data <- LinkPeaks( + object = data, + peak.assay = peak_assay, + expression.assay = expression_assay, + peak.layer = peak_layer, + expression.layer = expression_layer, + key = link_key, + gene.coords = gene_coords_chr, + genes.use = gene_set_chr, + distance = args$distance, + min.distance = args$min_distance, + min.cells = args$min_cells, + n_sample = args$n_sample, + pvalue_cutoff = args$pvalue_cutoff, + score_cutoff = args$score_cutoff, + calculate.zscore = TRUE, + verbose = TRUE, + tss.match.strategy = tss_match_strategy ) - links_gr <- Links(data[[peak_assay]]) - links_df <- if (length(links_gr) > 0) { - as.data.frame(links_gr)[, c("peak", "gene", "score", "pvalue")] - } else { - data.frame(peak = character(), gene = character(), - score = numeric(), pvalue = numeric()) - } + links_gr <- get_links_normalized(data[[peak_assay]], key = link_key) + all.link.objects[[level]] <- links_gr - data.table::fwrite( - links_df, - file = links_outfile, - sep = "\t", - quote = FALSE + links_df <- write_chromosome_tsv( + links = links_gr, + chromosome = level, + output.dir = output_dir, + key = link_key ) + all.outputs[[level]] <- links_df + message(" Wrote links to: ", links_outfile, " (", nrow(links_df), " rows)") } +combined_links_df <- data.table::rbindlist(all.outputs, use.names = TRUE, fill = TRUE) +data.table::fwrite(combined_links_df, combined_links_tsv, sep = "\t", quote = FALSE) +message("Wrote combined links TSV to: ", combined_links_tsv, " (", nrow(combined_links_df), " rows)") + +combined_links_gr <- combine_links_objects(all.link.objects) +Links(object = data[[peak_assay]], key = link_key) <- combined_links_gr +saveRDS(data, linked_seurat_rds) +message("Wrote linked Seurat RDS to: ", linked_seurat_rds) + message("All chromosomes processed.") diff --git a/linkage_analysis/scripts/02_run_aggregate_links.sh b/linkage_analysis/scripts/02_aggregate_and_validate_links.sh similarity index 54% rename from linkage_analysis/scripts/02_run_aggregate_links.sh rename to linkage_analysis/scripts/02_aggregate_and_validate_links.sh index 738a2f0..e004447 100644 --- a/linkage_analysis/scripts/02_run_aggregate_links.sh +++ b/linkage_analysis/scripts/02_aggregate_and_validate_links.sh @@ -5,7 +5,7 @@ #SBATCH --nodes=1 #SBATCH --ntasks=1 #SBATCH --cpus-per-task=1 -#SBATCH --mem=20G +#SBATCH --mem=30G #SBATCH --output=./log/%x-%j.out #SBATCH --error=./log/%x-%j.err #SBATCH -p cpu @@ -13,19 +13,24 @@ # ------------------------------------------------------------------------------ PROJECT_ROOT="/Projects/linkage_analysis" # "path/to/renv_proj" -LINKS_DIR_CSV="${PROJECT_ROOT}/scripts/config/objects_and_links.csv" # modify example table of seurat obj fpaths and their mutiome_link dirs -FILTER_CONDS_CSV="${PROJECT_ROOT}/scripts/config/link_filtering_conds.csv" # specify pvals and scores for links to keep (multiple conds possible) +LINKS_DIR_CSV="${PROJECT_ROOT}/scripts/config/objects_and_links.csv" # per-dataset multiome link directories +FILTER_CONDS_CSV="${PROJECT_ROOT}/scripts/config/link_filtering_conds.csv" # pvals/scores for links to keep # ------------------------------------------------------------------------------ -R_SCRIPT="${PROJECT_ROOT}/scripts/02_aggregate_links.R" +SCRIPT_AGGREGATE="${PROJECT_ROOT}/scripts/02_aggregate_links.R" +SCRIPT_GTEX_VALIDATE="${PROJECT_ROOT}/scripts/02a_validate_links_against_gtex.R" GENCODE_GTF="${PROJECT_ROOT}/resources/gencode.v32.basic.annotation.gtf.gz" +GTEX_FILE="${PROJECT_ROOT}/resources/GTEx_49tissues_release1.tsv.gz" COMBINED_OUTDIR="${PROJECT_ROOT}/data/combined_metadata" FILTERED_PARENT_DIR="${PROJECT_ROOT}/data/filtered_links" +GTEX_OUTPUT_DIR="${PROJECT_ROOT}/data/gtex_validation" DISTANCE_TO_TSS="500" -OVERWRITE="1" # set to 1 to pass --overwrite +PIP_CUTOFF="0.5" +GTEX_TISSUES="Whole_Blood" # use ANY for all tissues +OVERWRITE="1" # set to 1 to pass --overwrite to both steps ## Activate env with correct vers of renv installed in HPC # module load mambaforge @@ -35,8 +40,8 @@ OVERWRITE="1" # set to 1 to pass --overwrite mkdir -p ./log cd "${PROJECT_ROOT}" || exit 1 # ensures renv auto-loads via .Rprofile, plus paths are consistent -cmd=( - Rscript "${R_SCRIPT}" +cmd_aggregate=( + Rscript "${SCRIPT_AGGREGATE}" --links_dir_csv "${LINKS_DIR_CSV}" --filter_conds_csv "${FILTER_CONDS_CSV}" --combined_outdir "${COMBINED_OUTDIR}" @@ -45,12 +50,29 @@ cmd=( --distance_to_tss "${DISTANCE_TO_TSS}" ) +cmd_validate=( + Rscript "${SCRIPT_GTEX_VALIDATE}" + --links_dir_csv "${LINKS_DIR_CSV}" + --combined_outdir "${COMBINED_OUTDIR}" + --gencode_gtf "${GENCODE_GTF}" + --gtex_file "${GTEX_FILE}" + --output_dir "${GTEX_OUTPUT_DIR}" + --pip_cutoff "${PIP_CUTOFF}" + --gtex_tissues "${GTEX_TISSUES}" +) + if [[ "${OVERWRITE}" == "1" ]]; then - cmd+=(--overwrite) + cmd_aggregate+=(--overwrite) + cmd_validate+=(--overwrite) fi -echo "Running command:" -printf ' %q ' "${cmd[@]}" +echo "Running aggregate step:" +printf ' %q ' "${cmd_aggregate[@]}" +echo +"${cmd_aggregate[@]}" echo -"${cmd[@]}" +echo "Running GTEx validation step:" +printf ' %q ' "${cmd_validate[@]}" +echo +"${cmd_validate[@]}" diff --git a/linkage_analysis/scripts/02_aggregate_links.R b/linkage_analysis/scripts/02_aggregate_links.R index ff13e0d..0f710ff 100644 --- a/linkage_analysis/scripts/02_aggregate_links.R +++ b/linkage_analysis/scripts/02_aggregate_links.R @@ -1,9 +1,9 @@ +#!/usr/bin/env Rscript + suppressPackageStartupMessages({ library(argparse) library(Signac) - library(Seurat) library(GenomicRanges) - library(Matrix) library(data.table) library(rtracklayer) library(REMO.v1.GRCh38) @@ -25,7 +25,7 @@ parser$add_argument("--filtered_parent_dir", required = TRUE, parser$add_argument("--gencode_gtf", required = TRUE, help = "GENCODE GTF") parser$add_argument("--distance_to_tss", default = 500, - help = "Distance used in Signac:::DistanceToTSS (proximal genes)") + help = "Distance around TSS used to annotate proximal genes") parser$add_argument("--overwrite", action = "store_true", help = "Overwrite per-chr combined + per-condition outputs if they exist") @@ -45,39 +45,92 @@ dir.create(filtered_parent_dir, recursive = TRUE, showWarnings = FALSE) # ----------------------------- # Helpers # ----------------------------- -CollapseDuplicateColumns <- function(mat) { - if (!inherits(mat, "dgCMatrix")) stop("Input must be dgCMatrix.") - columns <- colnames(mat) - if (is.null(columns)) stop("Matrix must have column names.") - u_colnames <- unique(columns) - - T <- as(mat, "TsparseMatrix") - if (length(T@x) == 0L) { - out <- Matrix::sparseMatrix( - i = integer(0), j = integer(0), x = numeric(0), - dims = c(nrow(mat), length(u_colnames)) - ) - rownames(out) <- rownames(mat) - colnames(out) <- u_colnames - return(out) +collapse_vec <- function(x) { + x <- as.character(unlist(x, use.names = FALSE)) + x <- x[!is.na(x) & nzchar(x)] + if (length(x) == 0L) "" else paste(unique(x), collapse = ",") +} + +merge_gene_vecs <- function(...) { + x <- as.character(unlist(list(...), use.names = FALSE)) + x <- x[!is.na(x) & nzchar(x)] + unique(x) +} + +get_remo_ids <- function(gr) { + if ("REMO" %in% names(S4Vectors::mcols(gr))) { + return(as.character(gr$REMO)) + } + if ("name" %in% names(S4Vectors::mcols(gr))) { + return(as.character(gr$name)) } + stop("REMO.v1.GRCh38 must contain either a 'REMO' or 'name' metadata column") +} + +get_unique_symbol_transcripts <- function(gencode, rel_levels) { + gencode_genes <- gencode[ + gencode$type == "gene" & + GenomicRanges::seqnames(gencode) %in% rel_levels + ] + + gene_dt <- data.table( + seqnames = as.character(GenomicRanges::seqnames(gencode_genes)), + start = GenomicRanges::start(gencode_genes), + end = GenomicRanges::end(gencode_genes), + strand = as.character(GenomicRanges::strand(gencode_genes)), + gene_id = as.character(gencode_genes$gene_id), + gene_name = as.character(gencode_genes$gene_name) + )[!is.na(gene_name) & gene_name != ""] + + multi_chr_genes <- gene_dt[ + , .N, by = .(gene_name, seqnames) + ][, .N, by = gene_name][N > 1L, gene_name] + + multi_gene_same_chr <- gene_dt[ + , .N, by = .(gene_name, seqnames, gene_id, start, end, strand) + ][, .N, by = .(gene_name, seqnames)][N > 1L, gene_name] + + nonuniq_symbols <- unique(c(multi_chr_genes, multi_gene_same_chr)) + + gencode[ + gencode$type %in% c("transcript", "mRNA") & + GenomicRanges::seqnames(gencode) %in% rel_levels & + !(as.character(gencode$gene_name) %in% nonuniq_symbols) + ] +} - group <- match(columns[T@j + 1L], u_colnames) - df <- data.frame(i = T@i + 1L, j = group, x = T@x) - agg <- stats::aggregate(x ~ i + j, data = df, FUN = max) +get_proximal_genes_by_peak <- function(peaks, transcripts, distance_to_tss) { + peak_names <- Signac::GRangesToString(peaks) - out <- Matrix::sparseMatrix( - i = agg$i, j = agg$j, x = agg$x, - dims = c(nrow(mat), length(u_colnames)) + if (length(transcripts) == 0L) { + out <- data.table(peak = character(), proximal_genes = list()) + setkey(out, peak) + return(out) + } + + tss <- GenomicRanges::resize(transcripts, width = 1, fix = "start") + tss_window <- GenomicRanges::resize( + tss, + width = 2L * as.integer(distance_to_tss) + 1L, + fix = "center" ) - colnames(out) <- u_colnames - rownames(out) <- rownames(mat) - out -} + GenomicRanges::start(tss_window) <- pmax(1L, GenomicRanges::start(tss_window)) -collapse_vec <- function(x) { - x <- as.character(unlist(x)) - if (length(x) == 0L) "" else paste(x, collapse = ",") + hits <- GenomicRanges::findOverlaps(peaks, tss_window, type = "any") + if (length(hits) == 0L) { + out <- data.table(peak = character(), proximal_genes = list()) + setkey(out, peak) + return(out) + } + + hit_dt <- data.table( + peak = peak_names[S4Vectors::queryHits(hits)], + gene = as.character(transcripts$gene_name[S4Vectors::subjectHits(hits)]) + )[!is.na(gene) & gene != ""] + + out <- hit_dt[, .(proximal_genes = list(unique(gene))), by = peak] + setkey(out, peak) + out } # ----------------------------- @@ -99,56 +152,46 @@ if (overwrite && file.exists(summary_fpath)) file.remove(summary_fpath) rel_levels <- seqlevels(REMO.v1.GRCh38) -peak_strs <- GRangesToString(REMO.v1.GRCh38) -peak_grps <- REMO.v1.GRCh38$REMO +peak_strs <- Signac::GRangesToString(REMO.v1.GRCh38) +peak_grps <- get_remo_ids(REMO.v1.GRCh38) names(peak_grps) <- peak_strs -gencode <- import(gencode_gtf) -gencode_transcripts <- gencode[gencode$type == "transcript" & seqnames(gencode) %in% rel_levels] - -dt_gc <- as.data.table(gencode_transcripts)[, .(gene_name, seqnames = as.character(seqnames))] -multi_chr_genes <- dt_gc[, .(n_chr = uniqueN(seqnames)), by = gene_name][n_chr > 1] - -pdm <- Signac:::DistanceToTSS(peaks = REMO.v1.GRCh38, genes = gencode_transcripts, distance = distance_to_tss) -pdm <- CollapseDuplicateColumns(pdm) - -T <- as(pdm, "TsparseMatrix") -dt_long <- data.table( - peak = rownames(pdm)[T@i + 1L], - gene = colnames(pdm)[T@j + 1L] +gencode <- rtracklayer::import(gencode_gtf) +gencode_transcripts <- get_unique_symbol_transcripts(gencode, rel_levels) +pcre_by_peak <- get_proximal_genes_by_peak( + peaks = REMO.v1.GRCh38, + transcripts = gencode_transcripts, + distance_to_tss = distance_to_tss ) -dt_hits <- dt_long[, .(genes = list(unique(gene))), by = peak] - -pcre_by_peak <- dt_hits[, .(proximal_genes = list(unique(genes))), by = peak] -setkey(pcre_by_peak, peak) # ----------------------------- # Combine PGLs per chr across datasets + add metadata # ----------------------------- +required_link_cols <- c("peak", "gene", "score", "pvalue") chr_links_list <- list() for (chr in rel_levels) { links_list <- vector("list", length(link_dirs)) + for (i in seq_along(link_dirs)) { link_dir <- link_dirs[[i]] link_fpath <- file.path(link_dir, paste0(chr, ".tsv")) links <- fread(link_fpath) - # Ensure that genes with the same name across multiple chr - # are filtered out of downstream steps - links <- links[!(gene %in% multi_chr_genes$gene_name), ] - n_before <- nrow(links) - is_multi <- links$gene %in% multi_chr_genes$gene_name - links <- links[!is_multi] - - cat(sprintf( - "Removed %d multi_chr_genes; %d/%d links removed\n", - length(unique(links$gene[is_multi])), - sum(is_multi), - n_before - )) - links[, remo_name := peak_grps[peak]] - links[, dataset := basename(link_dir)] + missing_cols <- setdiff(required_link_cols, names(links)) + if (length(missing_cols) > 0L) { + stop( + "Link file is missing required columns: ", link_fpath, + "\nMissing: ", paste(missing_cols, collapse = ", ") + ) + } + + links[, `:=`( + score = as.numeric(score), + pvalue = as.numeric(pvalue), + remo_name = peak_grps[peak], + dataset = basename(link_dir) + )] links_list[[i]] <- links } @@ -156,17 +199,19 @@ for (chr in rel_levels) { chr_links <- rbindlist(links_list, use.names = TRUE, fill = TRUE)[order(peak, gene)] chr_links_list[[chr]] <- chr_links - chr_all_peaks <- REMO.v1.GRCh38[seqnames(REMO.v1.GRCh38) == chr] - cat(sprintf("%s: %d linked peaks out of %d peaks\n", - chr, length(unique(chr_links$peak)), length(chr_all_peaks)), - file = summary_fpath, append = TRUE) - cat(sprintf("%s: %d linked modules out of %d modules\n\n", - chr, length(unique(chr_links$remo_name)), length(unique(chr_all_peaks$name))), - file = summary_fpath, append = TRUE) + chr_all_peaks <- REMO.v1.GRCh38[GenomicRanges::seqnames(REMO.v1.GRCh38) == chr] + cat(sprintf( + "%s: %d linked peaks out of %d peaks\n", + chr, length(unique(chr_links$peak)), length(chr_all_peaks) + ), file = summary_fpath, append = TRUE) + cat(sprintf( + "%s: %d linked modules out of %d modules\n\n", + chr, length(unique(chr_links$remo_name)), length(unique(get_remo_ids(chr_all_peaks))) + ), file = summary_fpath, append = TRUE) outlink_fpath <- file.path(combined_outdir, paste0(chr, ".tsv")) if (overwrite || !file.exists(outlink_fpath)) { - write.table(chr_links, outlink_fpath, sep = "\t", quote = FALSE, row.names = FALSE) + fwrite(chr_links, outlink_fpath, sep = "\t", quote = FALSE) } } @@ -177,8 +222,8 @@ for (cond in conds) { cnd <- config[condition == cond] p_value <- cnd$pvalue corr_coeff <- cnd$score - if (is.null(p_value) || is.na(p_value)) p_value <- 1 - if (is.null(corr_coeff) || is.na(corr_coeff)) corr_coeff <- 0 + if (length(p_value) == 0L || is.na(p_value)) p_value <- 1 + if (length(corr_coeff) == 0L || is.na(corr_coeff)) corr_coeff <- 0 output_dir <- file.path(filtered_parent_dir, cond) pgl_metadata_dir <- file.path(output_dir, "pgl_metadata") @@ -195,13 +240,13 @@ for (cond in conds) { chr_pgl_fpath <- file.path(pgl_metadata_dir, paste0(chr, ".tsv")) if (overwrite || !file.exists(chr_pgl_fpath)) { - write.table(sig_links, chr_pgl_fpath, sep = "\t", quote = FALSE, row.names = FALSE) + fwrite(sig_links, chr_pgl_fpath, sep = "\t", quote = FALSE) } all_pgl_metadata_list[[chr]] <- sig_links - chr_peak_strs <- GRangesToString(REMO.v1.GRCh38[seqnames(REMO.v1.GRCh38) == chr]) - chr_remos <- REMO.v1.GRCh38[seqnames(REMO.v1.GRCh38) == chr]$name - peaks_dt <- data.table(peak = chr_peak_strs, remo_name = chr_remos) + chr_gr <- REMO.v1.GRCh38[GenomicRanges::seqnames(REMO.v1.GRCh38) == chr] + chr_peak_strs <- Signac::GRangesToString(chr_gr) + peaks_dt <- data.table(peak = chr_peak_strs, remo_name = get_remo_ids(chr_gr)) chr_genesPerPeak <- pcre_by_peak[peaks_dt, on = "peak"] @@ -213,17 +258,18 @@ for (cond in conds) { setkey(dcre_by_peak, peak) chr_genesPerPeak <- dcre_by_peak[chr_genesPerPeak, on = "peak"] - chr_genesPerPeak[, all_genes := Map(function(p, d) unique(c(p, d)), proximal_genes, distal_genes)] + chr_genesPerPeak[, all_genes := Map(merge_gene_vecs, proximal_genes, distal_genes)] cols <- c("proximal_genes", "distal_genes", "all_genes") - chr_genesPerPeak[, (cols) := lapply(.SD, function(col) vapply(col, collapse_vec, character(1L))), - .SDcols = cols] + chr_genesPerPeak[, (cols) := lapply(.SD, function(col) { + vapply(col, collapse_vec, character(1L)) + }), .SDcols = cols] setcolorder(chr_genesPerPeak, c("peak", "remo_name", "proximal_genes", "distal_genes", "all_genes")) chr_gpp_fpath <- file.path(genes_per_peak_dir, paste0(chr, ".tsv")) if (overwrite || !file.exists(chr_gpp_fpath)) { - write.table(chr_genesPerPeak, chr_gpp_fpath, sep = "\t", quote = FALSE, row.names = FALSE) + fwrite(chr_genesPerPeak, chr_gpp_fpath, sep = "\t", quote = FALSE) } all_genesPerPeak_list[[chr]] <- chr_genesPerPeak } @@ -231,12 +277,12 @@ for (cond in conds) { all_flinks <- rbindlist(all_pgl_metadata_list, use.names = TRUE, fill = TRUE) all_flinks_fpath <- file.path(pgl_metadata_dir, "pgls_metadata_all_chr.tsv") if (overwrite || !file.exists(all_flinks_fpath)) { - write.table(all_flinks, all_flinks_fpath, sep = "\t", quote = FALSE, row.names = FALSE) + fwrite(all_flinks, all_flinks_fpath, sep = "\t", quote = FALSE) } all_gpp <- rbindlist(all_genesPerPeak_list, use.names = TRUE, fill = TRUE) all_gpp_fpath <- file.path(genes_per_peak_dir, "genes_per_peak_all_chr.tsv") if (overwrite || !file.exists(all_gpp_fpath)) { - write.table(all_gpp, all_gpp_fpath, sep = "\t", quote = FALSE, row.names = FALSE) + fwrite(all_gpp, all_gpp_fpath, sep = "\t", quote = FALSE) } } diff --git a/linkage_analysis/scripts/02a_validate_links_against_gtex.R b/linkage_analysis/scripts/02a_validate_links_against_gtex.R new file mode 100644 index 0000000..13825dc --- /dev/null +++ b/linkage_analysis/scripts/02a_validate_links_against_gtex.R @@ -0,0 +1,580 @@ +#!/usr/bin/env Rscript + +suppressPackageStartupMessages({ + library(argparse) + library(data.table) + library(ggplot2) + library(GenomicRanges) + library(IRanges) + library(S4Vectors) + library(rtracklayer) +}) + +# ----------------------------- +# Args +# ----------------------------- +parser <- ArgumentParser( + description = "Validate per-dataset and combined peak-gene links against GTEx/SGL fine-mapping positives" +) + +parser$add_argument("--links_dir_csv", required = TRUE, + help = "CSV used by 02_aggregate_links.R. Required column: links_dir; optional column: dataset") +parser$add_argument("--combined_outdir", required = TRUE, + help = "Directory written by 02_aggregate_links.R containing combined per-chr link TSVs") +parser$add_argument("--gencode_gtf", required = TRUE, + help = "GENCODE GTF used for mapping GTEx ENSG IDs to unique gene symbols") +parser$add_argument("--gtex_file", required = TRUE, + help = "GTEx/SGL fine-mapping TSV/TSV.GZ with columns gene,chromosome,end,tissue,pip") +parser$add_argument("--output_dir", required = TRUE, + help = "Output directory for GTEx validation plots and tables") +parser$add_argument("--pip_cutoff", type = "double", default = 0.5, + help = "Minimum GTEx/SGL PIP used to define positive variants [default: 0.5]") +parser$add_argument("--gtex_tissues", default = "Whole_Blood", + help = "Comma-separated GTEx tissue filter; use ANY to include all tissues [default: Whole_Blood]") +parser$add_argument("--use_signed_score_metrics", action = "store_true", + help = "Use signed score and zscore instead of absolute values [default: use absolute values]") +parser$add_argument("--plot_distributions", action = "store_true", + help = "Also plot metric distributions by GTEx support") +parser$add_argument("--overwrite", action = "store_true", + help = "Overwrite existing validation outputs") + +args <- parser$parse_args() + +links_dir_csv <- args$links_dir_csv +combined_outdir <- args$combined_outdir +gencode_gtf <- args$gencode_gtf +gtex_file <- args$gtex_file +output_dir <- args$output_dir +pip_cutoff <- args$pip_cutoff +gtex_tissues <- if (toupper(trimws(args$gtex_tissues)) %in% c("ANY", "ALL", "")) { + NULL +} else { + trimws(strsplit(args$gtex_tissues, ",", fixed = TRUE)[[1]]) +} +use_abs_score_metrics <- !isTRUE(args$use_signed_score_metrics) +plot_distributions <- isTRUE(args$plot_distributions) +overwrite <- isTRUE(args$overwrite) + +dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + +# ----------------------------- +# Helpers +# ----------------------------- +check_cols <- function(dt, cols, name) { + missing <- setdiff(cols, names(dt)) + if (length(missing) > 0L) { + stop(name, " is missing columns: ", paste(missing, collapse = ", ")) + } +} + +strip_ensg_version <- function(x) sub("\\..*$", "", x) + +as_ucsc_chr <- function(x) { + x <- as.character(x) + x <- sub("^CHR", "chr", x, ignore.case = TRUE) + x <- ifelse(grepl("^chr", x), x, paste0("chr", x)) + sub("chrMT$", "chrM", x) +} + +make_peak_key <- function(seqnames, start, end) { + paste(as_ucsc_chr(seqnames), as.integer(start), as.integer(end), sep = "-") +} + +normalize_peak_key <- function(x) { + x <- as.character(x) + x <- gsub(":", "-", x, fixed = TRUE) + x <- gsub("_", "-", x, fixed = TRUE) + parts <- tstrsplit(x, "-", fixed = TRUE) + if (length(parts) < 3L) { + stop("Could not parse peak IDs like chr-start-end from peak column.") + } + make_peak_key(parts[[1]], parts[[2]], parts[[3]]) +} + +peak_key_to_granges <- function(peak_key) { + peak_key <- unique(as.character(peak_key)) + parts <- tstrsplit(peak_key, "-", fixed = TRUE) + if (length(parts) < 3L) { + stop("Could not parse normalized peak keys like chr-start-end.") + } + gr <- GRanges( + seqnames = parts[[1]], + ranges = IRanges(start = as.integer(parts[[2]]), end = as.integer(parts[[3]])) + ) + gr$peak_key <- peak_key + gr +} + +safe_label <- function(x) { + x <- gsub("[^A-Za-z0-9_.-]+", "_", x) + x <- gsub("_+", "_", x) + sub("_$", "", sub("^_", "", x)) +} + +make_gtf_maps <- function(gtf_path) { + message("Importing GTF and building unique ENSG-symbol map...") + gtf <- rtracklayer::import(gtf_path) + gene_gr <- gtf[S4Vectors::mcols(gtf)$type == "gene"] + + gene_dt <- data.table( + seqnames = as.character(GenomicRanges::seqnames(gene_gr)), + start = GenomicRanges::start(gene_gr), + end = GenomicRanges::end(gene_gr), + strand = as.character(GenomicRanges::strand(gene_gr)), + gene_id = as.character(S4Vectors::mcols(gene_gr)$gene_id), + gene_name = as.character(S4Vectors::mcols(gene_gr)$gene_name) + )[!is.na(gene_name) & gene_name != ""] + + gene_dt[, gene_id_stripped := strip_ensg_version(gene_id)] + + multi_chr <- gene_dt[ + , .N, by = .(gene_name, seqnames) + ][, .N, by = gene_name][N > 1L, gene_name] + + multi_gene_same_chr <- gene_dt[ + , .N, by = .(gene_name, seqnames, gene_id_stripped, start, end, strand) + ][, .N, by = .(gene_name, seqnames)][N > 1L, gene_name] + + nonuniq_symbols <- unique(c(multi_chr, multi_gene_same_chr)) + gene_dt <- gene_dt[!(gene_name %in% nonuniq_symbols)] + + ensg_to_symbol <- gene_dt$gene_name + names(ensg_to_symbol) <- gene_dt$gene_id_stripped + + list( + ensg_to_symbol = ensg_to_symbol, + nonuniq_symbols = nonuniq_symbols + ) +} + +read_link_file <- function(path, dataset = NULL) { + links <- fread(path) + if (nrow(links) == 0L) return(links) + + check_cols(links, c("peak", "gene", "score"), path) + + if (!("zscore" %in% names(links))) { + if (all(c("pvalue", "score") %in% names(links))) { + message("zscore absent in ", path, "; deriving signed zscore from pvalue and score.") + links[, zscore := sign(as.numeric(score)) * abs(qnorm(as.numeric(pvalue) / 2, lower.tail = FALSE))] + max_z <- max(abs(links[is.finite(zscore), zscore]), na.rm = TRUE) + if (is.finite(max_z)) { + links[is.infinite(zscore), zscore := sign(as.numeric(score)) * max_z] + } + } else { + stop("Link file must contain zscore, or pvalue + score so zscore can be derived: ", path) + } + } + + distance_cols <- intersect(c("tss_distance", "distance", "peak_tss_distance"), names(links)) + if (length(distance_cols) == 0L) { + stop("Link file must contain tss_distance, distance, or peak_tss_distance: ", path) + } + + links[, `:=`( + peak = as.character(peak), + gene = as.character(gene), + score = as.numeric(score), + zscore = as.numeric(zscore), + distance_bp = abs(as.numeric(get(distance_cols[[1L]]))), + peak_key = normalize_peak_key(peak), + source_file = basename(path) + )] + links[, chr := tstrsplit(peak_key, "-", fixed = TRUE)[[1L]]] + + if (!is.null(dataset) && !("dataset" %in% names(links))) { + links[, dataset := dataset] + } + + unique(links) +} + +load_links_dir <- function(link_dir, dataset) { + files <- list.files(link_dir, pattern = "^chr.*\\.tsv(\\.gz)?$", full.names = TRUE) + files <- files[!grepl("all_chromosomes\\.tsv(\\.gz)?$", basename(files))] + if (length(files) == 0L) { + stop("No per-chromosome chr*.tsv files found in ", link_dir) + } + rbindlist(lapply(files, read_link_file, dataset = dataset), use.names = TRUE, fill = TRUE) +} + +load_combined_outdir <- function(combined_dir) { + files <- list.files(combined_dir, pattern = "^chr.*\\.tsv(\\.gz)?$", full.names = TRUE) + files <- files[!grepl("all_chromosomes\\.tsv(\\.gz)?$", basename(files))] + if (length(files) == 0L) return(NULL) + rbindlist(lapply(files, read_link_file), use.names = TRUE, fill = TRUE) +} + +make_gtex_ground_truth <- function( + gtex_path, + link_peaks, + gene_map, + link_genes, + tissues = NULL, + pip.cutoff = 0.5 +) { + link_chromosomes <- sort(unique(as.character(GenomicRanges::seqnames(link_peaks)))) + + gtex <- fread( + gtex_path, + select = c("gene", "chromosome", "end", "tissue", "pip") + ) + check_cols(gtex, c("gene", "chromosome", "end", "tissue", "pip"), "GTEx/SGL file") + + gtex[, `:=`( + gene_ensg = strip_ensg_version(gene), + chromosome = as_ucsc_chr(chromosome), + pip = as.numeric(pip) + )] + gtex[, gene_symbol := gene_map$ensg_to_symbol[gene_ensg]] + + if (!is.null(tissues)) { + gtex <- gtex[tissue %in% tissues] + } + + gtex <- gtex[ + !is.na(gene_symbol) & + gene_symbol %in% link_genes & + chromosome %in% link_chromosomes & + chromosome != "chrY" & + !is.na(end) & + is.finite(pip) & + pip >= pip.cutoff + ] + + if (nrow(gtex) == 0L) { + warning("No GTEx/SGL positives remain after filtering. All evaluated links will be negatives.") + return(data.table( + peak_key = character(), gene = character(), max_pip = numeric(), + tissues = character(), n_gtex_variants = integer() + )) + } + + gtex_gr <- GRanges( + seqnames = gtex$chromosome, + ranges = IRanges(start = as.integer(gtex$end), end = as.integer(gtex$end)), + gene_symbol = gtex$gene_symbol, + tissue = gtex$tissue, + pip = gtex$pip + ) + + hits <- findOverlaps(gtex_gr, link_peaks, ignore.strand = TRUE) + if (length(hits) == 0L) { + warning("No GTEx/SGL positive variants overlap evaluated link peaks.") + return(data.table( + peak_key = character(), gene = character(), max_pip = numeric(), + tissues = character(), n_gtex_variants = integer() + )) + } + + hit_dt <- data.table( + peak_key = link_peaks$peak_key[S4Vectors::subjectHits(hits)], + gene = S4Vectors::mcols(gtex_gr)$gene_symbol[S4Vectors::queryHits(hits)], + tissue = S4Vectors::mcols(gtex_gr)$tissue[S4Vectors::queryHits(hits)], + pip = S4Vectors::mcols(gtex_gr)$pip[S4Vectors::queryHits(hits)] + ) + + hit_dt[ + , .( + max_pip = max(pip, na.rm = TRUE), + tissues = paste(sort(unique(tissue)), collapse = ";"), + n_gtex_variants = .N + ), + by = .(peak_key, gene) + ] +} + +fast_auc <- function(predictor, response) { + dt <- data.table(pred = predictor, resp = as.logical(response))[!is.na(pred) & !is.na(resp)] + if (nrow(dt) == 0L) return(NA_real_) + + n_pos <- as.numeric(sum(dt$resp)) + n_neg <- as.numeric(sum(!dt$resp)) + if (n_pos == 0 || n_neg == 0) return(NA_real_) + + r <- rank(dt$pred, ties.method = "average") + (sum(r[dt$resp]) - n_pos * (n_pos + 1) / 2) / (n_pos * n_neg) +} + +calc_roc <- function(predictor, response) { + dt <- data.table(pred = predictor, resp = as.logical(response))[!is.na(pred) & !is.na(resp)] + if (nrow(dt) == 0L || sum(dt$resp) == 0L || sum(!dt$resp) == 0L) { + return(data.table(pred = numeric(), tpr = numeric(), fpr = numeric())) + } + + setorder(dt, -pred) + dt_agg <- dt[, .(pos = sum(resp), neg = sum(!resp)), by = pred] + dt_agg[, `:=`( + tpr = cumsum(pos) / sum(pos), + fpr = cumsum(neg) / sum(neg) + )] + + rbind( + data.table(pred = Inf, tpr = 0, fpr = 0), + dt_agg[, .(pred, tpr, fpr)] + ) +} + +write_if_needed <- function(dt, path) { + if (overwrite || !file.exists(path)) { + fwrite(dt, path, sep = "\t", quote = FALSE) + } +} + +plot_roc_with_auc <- function(eval_dt, metric_spec, out_dir, file_label, display_label, plot_subtitle) { + auc_dt <- rbindlist(lapply(names(metric_spec), function(metric_name) { + col <- metric_spec[[metric_name]] + data.table( + metric = metric_name, + auc = fast_auc(eval_dt[[col]], eval_dt$is_gtex_valid) + ) + })) + write_if_needed(auc_dt, file.path(out_dir, paste0("roc_auc_metrics_", file_label, ".tsv"))) + + n_pos <- sum(eval_dt$is_gtex_valid, na.rm = TRUE) + n_neg <- sum(!eval_dt$is_gtex_valid, na.rm = TRUE) + if (n_pos == 0L || n_neg == 0L) { + warning("ROC cannot be plotted for ", display_label, ": ", n_pos, " positives and ", n_neg, " negatives.") + return(auc_dt) + } + + roc_dt <- rbindlist(lapply(names(metric_spec), function(metric_name) { + col <- metric_spec[[metric_name]] + auc <- auc_dt[metric == metric_name, auc] + roc <- calc_roc(eval_dt[[col]], eval_dt$is_gtex_valid) + if (nrow(roc) == 0L || !is.finite(auc)) return(NULL) + roc[, metric := sprintf("%s (AUC = %.3f)", metric_name, auc)] + roc + }), fill = TRUE) + + if (nrow(roc_dt) == 0L) return(auc_dt) + + p_roc <- ggplot(roc_dt, aes(x = fpr, y = tpr, color = metric)) + + geom_line(linewidth = 0.9) + + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray60") + + coord_equal(xlim = c(0, 1), ylim = c(0, 1), expand = FALSE) + + labs( + title = display_label, + subtitle = plot_subtitle, + x = "False positive rate", + y = "True positive rate", + color = NULL + ) + + theme_bw(base_size = 12) + + theme(legend.position = "bottom", legend.direction = "vertical") + + ggsave(file.path(out_dir, paste0("roc_metrics_", file_label, ".png")), p_roc, width = 7.5, height = 7, dpi = 300, bg = "white") + ggsave(file.path(out_dir, paste0("roc_metrics_", file_label, ".pdf")), p_roc, width = 7.5, height = 7, bg = "white") + + auc_dt +} + +plot_metric_distributions <- function(eval_dt, metric_spec, out_dir, file_label, display_label, plot_subtitle) { + long_dt <- melt( + eval_dt, + id.vars = "is_gtex_valid", + measure.vars = unname(metric_spec), + variable.name = "metric_col", + value.name = "value" + ) + inv_map <- setNames(names(metric_spec), unname(metric_spec)) + long_dt[, metric := inv_map[as.character(metric_col)]] + long_dt <- long_dt[is.finite(value)] + + p <- ggplot(long_dt, aes(x = is_gtex_valid, y = value, fill = is_gtex_valid)) + + geom_boxplot(outlier.shape = NA) + + facet_wrap(~ metric, scales = "free_y") + + labs( + title = display_label, + subtitle = paste(plot_subtitle, "metric distributions", sep = " — "), + x = "GTEx/SGL-supported link?", + y = "Metric value" + ) + + theme_bw(base_size = 12) + + theme(legend.position = "none") + + ggsave(file.path(out_dir, paste0("metric_distributions_", file_label, ".png")), p, width = 9, height = 6, dpi = 300, bg = "white") + ggsave(file.path(out_dir, paste0("metric_distributions_", file_label, ".pdf")), p, width = 9, height = 6, bg = "white") +} + +write_validation_diagnostics <- function(links_eval, truth_dt, out_dir, file_label) { + link_peak_genes <- unique(links_eval[, .(peak_key, gene)]) + truth_peak_genes <- unique(truth_dt[, .(peak_key, gene)]) + + diag_dt <- data.table( + diagnostic = c( + "evaluated_links", + "evaluated_unique_peak_gene_pairs", + "gtex_truth_peak_gene_pairs", + "gtex_supported_links_exact_peak_gene", + "links_not_gtex_supported_exact_peak_gene", + "link_peaks_with_any_gtex_variant", + "link_genes_with_any_gtex_variant" + ), + value = c( + nrow(links_eval), + nrow(link_peak_genes), + nrow(truth_peak_genes), + sum(links_eval$is_gtex_valid, na.rm = TRUE), + sum(!links_eval$is_gtex_valid, na.rm = TRUE), + length(intersect(unique(link_peak_genes$peak_key), unique(truth_peak_genes$peak_key))), + length(intersect(unique(link_peak_genes$gene), unique(truth_peak_genes$gene))) + ) + ) + + write_if_needed(diag_dt, file.path(out_dir, paste0("validation_diagnostics_", file_label, ".tsv"))) + diag_dt +} + +run_validation <- function(links_dt, out_dir, file_label, display_label, plot_subtitle, gtf_maps) { + dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + + if (nrow(links_dt) == 0L) { + warning("No links available for ", display_label) + return(NULL) + } + + links_dt <- copy(links_dt) + links_dt <- links_dt[!is.na(peak_key) & !is.na(gene) & gene != ""] + links_dt[, gene := as.character(gene)] + + link_peaks <- peak_key_to_granges(links_dt$peak_key) + link_genes <- sort(unique(links_dt$gene)) + link_chromosomes <- sort(unique(links_dt$chr)) + + truth_dt <- make_gtex_ground_truth( + gtex_path = gtex_file, + link_peaks = link_peaks, + gene_map = gtf_maps, + link_genes = link_genes, + tissues = gtex_tissues, + pip.cutoff = pip_cutoff + ) + + write_if_needed(truth_dt, file.path(out_dir, paste0("gtex_peak_gene_truth_", file_label, ".tsv"))) + + links_dt[, peak_gene := paste(peak_key, gene, sep = "__")] + truth_dt[, peak_gene := paste(peak_key, gene, sep = "__")] + links_dt[, is_gtex_valid := peak_gene %in% truth_dt$peak_gene] + + truth_annot <- truth_dt[, .(peak_gene, max_pip, tissues, n_gtex_variants)] + links_eval <- merge(links_dt, truth_annot, by = "peak_gene", all.x = TRUE, sort = FALSE) + links_eval[is.na(max_pip), `:=`(max_pip = 0, tissues = NA_character_, n_gtex_variants = 0L)] + + links_eval[, score_metric := if (use_abs_score_metrics) abs(score) else score] + links_eval[, zscore_metric := if (use_abs_score_metrics) abs(zscore) else zscore] + links_eval[, exp_distance_200kb_metric := exp(-distance_bp / 200000)] + links_eval[, score_exp_distance_200kb_metric := score_metric * exp_distance_200kb_metric] + + metric_spec <- c( + "score" = "score_metric", + "zscore" = "zscore_metric", + "exp(-distance / 200 kb)" = "exp_distance_200kb_metric", + "score * exp(-distance / 200 kb)" = "score_exp_distance_200kb_metric" + ) + + summary_dt <- data.table( + label = display_label, + chromosomes = paste(link_chromosomes, collapse = ","), + gtex_tissues = if (is.null(gtex_tissues)) "any" else paste(gtex_tissues, collapse = ","), + pip_cutoff = pip_cutoff, + use_abs_score_metrics = use_abs_score_metrics, + n_links = nrow(links_eval), + n_gtex_supported_links = sum(links_eval$is_gtex_valid), + n_gtex_unsupported_links = sum(!links_eval$is_gtex_valid), + n_link_peaks = length(unique(links_eval$peak_key)), + n_link_genes = length(unique(links_eval$gene)), + n_gtex_peak_gene_truth_pairs = nrow(truth_dt) + ) + + write_if_needed(summary_dt, file.path(out_dir, paste0("validation_summary_", file_label, ".tsv"))) + write_if_needed(links_eval, file.path(out_dir, paste0("links_with_gtex_labels_and_metrics_", file_label, ".tsv"))) + + diag_dt <- write_validation_diagnostics(links_eval, truth_dt, out_dir, file_label) + auc_dt <- plot_roc_with_auc( + eval_dt = links_eval, + metric_spec = metric_spec, + out_dir = out_dir, + file_label = file_label, + display_label = display_label, + plot_subtitle = plot_subtitle + ) + if (plot_distributions) { + plot_metric_distributions( + eval_dt = links_eval, + metric_spec = metric_spec, + out_dir = out_dir, + file_label = file_label, + display_label = display_label, + plot_subtitle = plot_subtitle + ) + } + + message(display_label, ": ", sum(links_eval$is_gtex_valid), " GTEx/SGL-supported links out of ", nrow(links_eval)) + invisible(list(summary = summary_dt, diagnostics = diag_dt, auc = auc_dt)) +} + +# ----------------------------- +# Main workflow +# ----------------------------- +message("Links config: ", links_dir_csv) +message("Combined outdir: ", combined_outdir) +message("GENCODE GTF: ", gencode_gtf) +message("GTEx file: ", gtex_file) +message("Output directory: ", output_dir) +message("PIP cutoff: ", pip_cutoff) +message("GTEx tissues: ", if (is.null(gtex_tissues)) "any" else paste(gtex_tissues, collapse = ", ")) +message("Score metrics: ", if (use_abs_score_metrics) "absolute" else "signed") + +links_cfg <- fread(links_dir_csv) +check_cols(links_cfg, "links_dir", "links_dir_csv") + +if ("dataset" %in% names(links_cfg)) { + dataset_names <- as.character(links_cfg$dataset) +} else { + dataset_names <- basename(normalizePath(links_cfg$links_dir, mustWork = FALSE)) +} +dataset_names[is.na(dataset_names) | dataset_names == ""] <- paste0("dataset_", which(is.na(dataset_names) | dataset_names == "")) +dataset_names <- make.unique(dataset_names) + +links_cfg[, dataset := dataset_names] + +gtf_maps <- make_gtf_maps(gencode_gtf) + +per_dataset_links <- vector("list", nrow(links_cfg)) +names(per_dataset_links) <- links_cfg$dataset + +for (i in seq_len(nrow(links_cfg))) { + dataset <- links_cfg$dataset[[i]] + link_dir <- links_cfg$links_dir[[i]] + message("\nLoading links for dataset: ", dataset) + links_dt <- load_links_dir(link_dir, dataset = dataset) + per_dataset_links[[dataset]] <- links_dt + + run_validation( + links_dt = links_dt, + out_dir = file.path(output_dir, "per_dataset", safe_label(dataset)), + file_label = safe_label(dataset), + display_label = dataset, + plot_subtitle = "GTEx/SGL ROC across all chromosomes", + gtf_maps = gtf_maps + ) +} + +message("\nLoading combined links across datasets...") +combined_links <- load_combined_outdir(combined_outdir) +if (is.null(combined_links)) { + message("No per-chr combined files found in combined_outdir; using the row-bound per-dataset links.") + combined_links <- rbindlist(per_dataset_links, use.names = TRUE, fill = TRUE) +} + +run_validation( + links_dt = combined_links, + out_dir = file.path(output_dir, "all_datasets"), + file_label = "all_datasets", + display_label = "All datasets", + plot_subtitle = "GTEx/SGL ROC across all aggregated links", + gtf_maps = gtf_maps +) + +message("\nDone. Validation outputs written to: ", output_dir) diff --git a/linkage_analysis/scripts/03_get_conditional_probs.R b/linkage_analysis/scripts/03_get_conditional_probs.R index b55b12b..29bdeb5 100644 --- a/linkage_analysis/scripts/03_get_conditional_probs.R +++ b/linkage_analysis/scripts/03_get_conditional_probs.R @@ -1,26 +1,24 @@ #!/usr/bin/env Rscript suppressPackageStartupMessages({ - library(argparse) - library(data.table) - library(rtracklayer) - library(GenomicRanges) - library(Signac) - library(REMO.v1.GRCh38) + library(argparse) + library(data.table) + library(Signac) + library(REMO.v1.GRCh38) }) parser <- ArgumentParser() parser$add_argument("--renv_dir", required = TRUE, - help = "Path to renv project root") + help = "Path to renv project root") parser$add_argument("--sig_pgl_fpath", required = TRUE, - help = "Significant PGL metadata TSV path (must include column 'peak')") + help = "Significant PGL metadata TSV path (must include column 'peak')") parser$add_argument("--output_dir", required = TRUE, - help = "Output directory") + help = "Output directory") parser$add_argument("--output_name", default = "conditional_prob_unlinked_singleCRE.tsv", - help = "Output TSV filename") + help = "Output TSV filename") parser$add_argument("--overwrite", action = "store_true", - help = "Overwrite output file if exists") + help = "Overwrite output file if exists") args <- parser$parse_args() @@ -34,18 +32,18 @@ overwrite <- isTRUE(args$overwrite) dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) if (file.exists(output_fpath) && !overwrite) { - message(sprintf("Output exists and --overwrite not set, skipping: %s", output_fpath)) - quit(save = "no", status = 0) + message(sprintf("Output exists and --overwrite not set, skipping: %s", output_fpath)) + quit(save = "no", status = 0) } -remo <- import(remo_bed, format = "BED") -start(remo) <- start(remo) - 1 - sig_links <- fread(sig_pgl_fpath) +if (!("peak" %in% names(sig_links))) { + stop("sig_pgl_fpath must contain a 'peak' column", call. = FALSE) +} gr_dt <- as.data.table(REMO.v1.GRCh38) -gr_dt[, peak := GRangesToString(REMO.v1.GRCh38)] -gr_dt[, remo_name := REMO.v1.GRCh38$REMO] +gr_dt[, peak := Signac::GRangesToString(REMO.v1.GRCh38)] +gr_dt[, remo_name := as.character(REMO.v1.GRCh38$REMO)] linked_peaks <- unique(sig_links$peak) gr_dt[, linked := peak %in% linked_peaks] @@ -58,46 +56,40 @@ gr_dt[, singleCRE := (n_peaks == 1L)] freq <- gr_dt[, .N, by = .(linked, singleCRE)] -N_ls <- freq[linked == TRUE & singleCRE == TRUE, N] -N_lm <- freq[linked == TRUE & singleCRE == FALSE, N] -N_us <- freq[linked == FALSE & singleCRE == TRUE, N] -N_um <- freq[linked == FALSE & singleCRE == FALSE, N] +get_n <- function(linked_value, single_value) { + n <- freq[linked == linked_value & singleCRE == single_value, N] + if (length(n) == 0L) 0L else n +} -N_ls <- ifelse(length(N_ls) == 0L, 0L, N_ls) -N_lm <- ifelse(length(N_lm) == 0L, 0L, N_lm) -N_us <- ifelse(length(N_us) == 0L, 0L, N_us) -N_um <- ifelse(length(N_um) == 0L, 0L, N_um) +N_ls <- get_n(TRUE, TRUE) +N_lm <- get_n(TRUE, FALSE) +N_us <- get_n(FALSE, TRUE) +N_um <- get_n(FALSE, FALSE) N_total <- N_ls + N_lm + N_us + N_um N_linked <- N_ls + N_lm N_unlinked <- N_us + N_um N_singleCRE <- N_ls + N_us -N_multiCRE <- N_lm + N_um -P_unlinked <- N_unlinked / N_total -P_unlinked_given_singleCRE <- N_us / N_singleCRE -P_singleCRE_given_unlinked <- N_us / N_unlinked -P_linked_given_singleCRE <- N_ls / N_singleCRE -P_singleCRE_given_linked <- N_ls / N_linked -P_singleCRE <- N_singleCRE / N_total +safe_div <- function(num, den) ifelse(den == 0, NA_real_, num / den) cond_summary <- data.table( - metric = c( - "P(unlinked)", - "P(unlinked | singleCRE)", - "P(singleCRE | unlinked)", - "P(linked | singleCRE)", - "P(singleCRE | linked)", - "P(singleCRE)" - ), - value = c( - P_unlinked, - P_unlinked_given_singleCRE, - P_singleCRE_given_unlinked, - P_linked_given_singleCRE, - P_singleCRE_given_linked, - P_singleCRE - ) + metric = c( + "P(unlinked)", + "P(unlinked | singleCRE)", + "P(singleCRE | unlinked)", + "P(linked | singleCRE)", + "P(singleCRE | linked)", + "P(singleCRE)" + ), + value = c( + safe_div(N_unlinked, N_total), + safe_div(N_us, N_singleCRE), + safe_div(N_us, N_unlinked), + safe_div(N_ls, N_singleCRE), + safe_div(N_ls, N_linked), + safe_div(N_singleCRE, N_total) + ) ) print(cond_summary) diff --git a/linkage_analysis/scripts/03a_get_glm_coeffs.R b/linkage_analysis/scripts/03a_get_glm_coeffs.R index 3a8775a..0512de8 100644 --- a/linkage_analysis/scripts/03a_get_glm_coeffs.R +++ b/linkage_analysis/scripts/03a_get_glm_coeffs.R @@ -3,8 +3,6 @@ suppressPackageStartupMessages({ library(argparse) library(data.table) - library(Matrix) - library(rtracklayer) library(GenomicRanges) library(Signac) library(REMO.v1.GRCh38) @@ -75,37 +73,10 @@ get_linkedPeaks_dt <- function(gr, pgl_dt) { setkey(genes_per_peak, peak_name) # ---- Build peak metadata dt ---- - mcols_names <- names(mcols(gr)) - if ("REMO" %in% mcols_names) { - names_vec <- as.character(gr$REMO) - - } else if (length(mcols_names) > 0L) { - id_col <- mcols_names[1L] - names_vec <- as.character(mcols(gr)[[id_col]]) - - cat(sprintf( - "[WARN] GRanges has no 'REMO' metadata column; using first metadata column '%s' for 'names'.\n", - id_col - )) - - prefix <- "REMOv1.GRCh38-" - ok <- startsWith(names_vec, prefix) - ok[is.na(ok)] <- FALSE - - if (!all(ok)) { - bad_idx <- which(!ok) - ex_bad <- head(unique(names_vec[bad_idx]), 5) - cat(sprintf( - "[WARN] %d/%d values in '%s' do not start with '%s'. Examples: %s\n", - length(bad_idx), length(names_vec), - id_col, prefix, - paste(ex_bad, collapse = " | ") - )) - } - - } else { - stop("GRanges has no metadata columns (mcols is empty); cannot derive REMO identifier.") + if (!("REMO" %in% names(mcols(gr)))) { + stop("REMO.v1.GRCh38 subset is expected to contain a 'REMO' metadata column.") } + names_vec <- as.character(gr$REMO) dt <- data.table( seqname = as.character(seqnames(gr)), @@ -352,9 +323,6 @@ only_chrs <- parse_chr_list(args$only_chrs) sig_links <- fread(sig_pgl_fpath) -if ("pvalue" %in% names(sig_links)) { - sig_links[, pvalue := formatC(as.numeric(pvalue), format = "e", digits = 3)] -} rel_levels <- seqlevels(REMO.v1.GRCh38) if (length(only_chrs) > 0L) { diff --git a/linkage_analysis/scripts/03b_build_glm_summary.R b/linkage_analysis/scripts/03b_build_glm_summary.R index c512259..5043a17 100644 --- a/linkage_analysis/scripts/03b_build_glm_summary.R +++ b/linkage_analysis/scripts/03b_build_glm_summary.R @@ -39,11 +39,6 @@ if (file.exists(out_fpath) && !overwrite) { quit(save = "no", status = 0) } -tsv_files <- list.files(glm_dir, pattern = "\\.tsv$", full.names = TRUE, recursive = FALSE) -if (length(tsv_files) == 0L) { - stop(sprintf("No .tsv files found in: %s", glm_dir), call. = FALSE) -} - # ----------------------------- # Find TSV files # ----------------------------- @@ -136,7 +131,6 @@ out_dt[, chr_num := NULL] # ----------------------------- # Write output # ----------------------------- -out_fpath <- file.path(glm_dir, "glm_results_all.tsv") fwrite(out_dt, out_fpath, sep = "\t", quote = FALSE, na = "NA") message(sprintf("Wrote: %s", out_fpath))