From 7068a36065f66bef1063105905a268d94f6e10fb Mon Sep 17 00:00:00 2001 From: JavenTyr Date: Fri, 3 Jul 2026 04:52:21 +0800 Subject: [PATCH] add FindCandidateLinks; mod LinkPeaks; mod CollapseToLongestTranscript - changed default peak-TSS matching strategy - added link attributes - LinkPeaks no longer runs zscore calcs by default - removed expectation for gene_biotype when using gencode annotations --- R/links.R | 1105 +++++++++++++++++++++++++++++++++++++------------ R/utilities.R | 40 +- 2 files changed, 874 insertions(+), 271 deletions(-) diff --git a/R/links.R b/R/links.R index 1c5be282..a484483a 100644 --- a/R/links.R +++ b/R/links.R @@ -189,71 +189,601 @@ ConnectionsToLinks <- function( return(gi) } + + +#' 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)) + ) + } + + .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.body.coords <- gene.body.coords[gene.body.coords$gene_id %in% genes] + gene.body.coords$gene_name <- gene.body.coords$gene_id + } else { + gene.body.coords <- gene.body.coords[gene.body.coords$gene_name %in% genes] + } + + 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 (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 + } + + 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] + + # Return only objects used downstream by LinkPeaks(). + list( + candidate.matrix = candidate.matrix, + candidate.links = candidate.links, + peaks = peaks, + gene.coords = gene.coords.use, + candidate.link.table = candidate.link.table, + peak.data = peak.data, + expression.data = expression.data + ) +} + #' 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 -#' the gene expression and accessibility of each peak within a given distance -#' from the gene TSS, and computes an expected correlation coefficient for each -#' peak given the GC content, accessibility, and length of the peak. The -#' expected coefficient values for the peak are then used to compute a z-score -#' and p-value. +#' 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 peak.layer Name of layer to pull chromatin data from -#' @param expression.assay Name of assay containing gene expression information -#' @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 GRanges object containing coordinates of genes in the -#' expression assay. If NULL, extract from gene annotations stored in the assay. -#' @param distance Distance threshold for peaks to include in regression model -#' @param min.distance Minimum distance between peak and TSS to include in -#' regression model. If NULL (default), no minimum distance is used. +#' @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 in the results. -#' @param genes.use Genes to test. If NULL, determine from expression assay. -#' @param n_sample Number of peaks to sample at random when computing the null -#' distribution. -#' @param pvalue_cutoff Minimum p-value required to retain a link. Links with a -#' p-value equal or greater than this value will be removed from the output. -#' @param score_cutoff Minimum absolute value correlation coefficient for a link -#' to be retained -#' @param gene.id Set to TRUE if genes in the expression assay are named -#' using gene IDs rather than gene names. -#' @param verbose Display messages -#' @param peak.slot Deprecated (use `peak.layer`) -#' @param expression.slot Deprecated (used `expression.layer`) +#' 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 stats pnorm sd -#' @importFrom Matrix sparseMatrix rowSums +#' @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 -#' @importMethodsFrom Matrix t -#' -#' @return Returns a Seurat object with results added to the `links` slot in the -#' assay, stored under the key specified in the function. The results are stored -#' as an [InteractionSet::GInteractions] object accessible via the [Links()] -#' function. This contains the [GenomicRanges::GRanges] for the pair of linked -#' regions (peak and gene), with `anchor1` corresponding to the peak region and -#' `anchor2` corresponding to the gene region linked to the peak. The following -#' metadata is also stored in the `GInteractions` object: -#' * `anchor2.gene_id`: The gene ID for the linked gene -#' * `anchor2.gene_name`: The name of the linked gene -#' * `score`: the correlation coefficient between the accessibility of the -#' peak and expression of the gene -#' * `zscore`: the z-score of the correlation coefficient, computed based on -#' the distribution of correlation coefficients from a set of background peaks -#' * `pvalue`: the p-value associated with the z-score for the link +#' @importFrom stats pnorm sd #' #' @export #' @concept links @@ -273,27 +803,31 @@ LinkPeaks <- function( 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" - ) + 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", @@ -302,6 +836,7 @@ LinkPeaks <- function( ) expression.layer <- expression.slot } + if (is_present(arg = peak.slot)) { deprecate_warn( when = "1.16.0", @@ -310,8 +845,9 @@ LinkPeaks <- function( ) peak.layer <- peak.slot } + features.match <- c("GC.percent", "count", "sequence.length") - + if (method == "pearson") { cor_method <- corSparse } else if (method == "spearman") { @@ -320,258 +856,307 @@ LinkPeaks <- function( stop("method can be one of 'pearson' or 'spearman'.") } - if (is.null(x = gene.coords)) { - annot <- Annotation(object = object[[peak.assay]]) - if (is.null(x = annot)) { - stop("Gene annotations not found") + 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) } - gene.coords <- CollapseToLongestTranscript( - ranges = annot - ) - } - 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) - } - if (!(peak.layer %in% Layers(object = object[[peak.assay]]))) { - stop("Requested peak layer not found") - } - peak.data <- LayerData(object = object[[peak.assay]], layer = peak.layer) - if (!(expression.layer %in% Layers(object = object[[expression.assay]]))) { - stop("Requested expression layer not found") } - expression.data <- LayerData( - object = object[[expression.assay]], layer = expression.layer + + 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 + } ) - peakcounts <- rowSums(x = peak.data > 0) - genecounts <- rowSums(x = expression.data > 0) - peaks.keep <- peakcounts > min.cells - genes.keep <- genecounts > min.cells - peak.data <- peak.data[peaks.keep, ] - if (!is.null(x = genes.use)) { - genes.keep <- intersect( - x = names(x = genes.keep[genes.keep]), y = genes.use - ) + + if (is.null(candidates)) { + Links(object = object[[peak.assay]], key = key) <- make_empty_links() + return(object) } - expression.data <- expression.data[genes.keep, , drop = FALSE] - genes <- rownames(x = expression.data) - if (gene.id) { - gene.coords.use <- gene.coords[gene.coords$gene_id %in% genes, ] - gene.coords.use$gene_name <- gene.coords.use$gene_id + + 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 { - gene.coords.use <- gene.coords[gene.coords$gene_name %in% genes, ] - } - if (length(x = gene.coords.use) == 0) { - stop("Could not find gene coordinates for requested genes") - } - if (length(x = gene.coords.use) < nrow(x = expression.data)) { - message( - "Found gene coordinates for ", - length(x = gene.coords.use), - " genes" - ) - } - peaks <- granges(x = object[[peak.assay]]) - peaks <- peaks[peaks.keep] - peak_distance_matrix <- DistanceToTSS( - peaks = peaks, - genes = gene.coords.use, - distance = distance - ) - if (!is.null(x = min.distance)) { - peak_distance_matrix_min <- DistanceToTSS( - peaks = peaks, - genes = gene.coords.use, - distance = min.distance - ) - peak_distance_matrix <- peak_distance_matrix - peak_distance_matrix_min + NULL } + if (sum(peak_distance_matrix) == 0) { - stop( - "No peaks fall within distance threshold\n", - "Have you set the proper genome and seqlevelsStyle for ", - peak.assay, - " assay?" - ) + 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(x = expression.data), - " genes and ", - sum(rowSums(x = peak_distance_matrix) > 0), + "Testing ", nrow(x = expression.data), + " genes and ", sum(rowSums(x = peak_distance_matrix) > 0), " peaks" ) } - genes.use <- colnames(x = peak_distance_matrix) - all.peaks <- rownames(x = peak.data) peak.data <- t(x = peak.data) - coef.vec <- c() - gene.vec <- c() - zscore.vec <- c() if (nbrOfWorkers() > 1) { mylapply <- future_lapply } else { - mylapply <- ifelse(test = verbose, yes = pblapply, no = lapply) + mylapply <- if (verbose) pblapply else lapply } - - # run in parallel across genes + 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]) - gene.chrom <- as.character(x = seqnames(x = gene.coords.use[i])) if (sum(peak.use) < 2) { - # no peaks close to gene return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) - } else { - 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)) - } else { - # select peaks at random with matching GC content and accessibility - # sample from peaks on a different chromosome to the gene - peaks.test <- rownames(x = coef.result) - trans.peaks <- all.peaks[ - !grepl(pattern = paste0("^", gene.chrom, "-"), x = all.peaks) - ] - meta.use <- meta.features[trans.peaks, ] - pk.use <- meta.features[peaks.test, ] - 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 - ) - } - ) - # run background correlations - unique.bg <- unique(x = unlist(x = bg.peaks)) - 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 + 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 ) - 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]], ] - 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 - zscore.vec <- c(zscore.vec, zscores) - gene.vec <- c(gene.vec, rep(i, length(x = coef.result))) - coef.vec <- c(coef.vec, coef.result) } - gc(verbose = FALSE) - pval.vec <- 2 * pnorm(q = -abs(x = zscore.vec)) - links.keep <- pval.vec < pvalue_cutoff - if (sum(x = links.keep) == 0) { - return(list("gene" = NULL, "coef" = NULL, "zscore" = NULL)) + ) + + 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 { - gene.vec <- gene.vec[links.keep] - coef.vec <- coef.vec[links.keep] - zscore.vec <- zscore.vec[links.keep] - return(list( - "gene" = gene.vec, - "coef" = coef.vec, - "zscore" = zscore.vec - )) + 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] + ) } ) - # combine results + 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 significant links found") - } + if (verbose) message("No links pass score_cutoff") + Links(object = object[[peak.assay]], key = key) <- make_empty_links() return(object) } - peak.key <- seq_along( - along.with = unique(x = names(x = coef.vec)) - ) - names(x = peak.key) <- unique(x = names(x = coef.vec)) - coef.matrix <- sparseMatrix( - i = gene.vec, - j = peak.key[names(x = coef.vec)], - x = coef.vec, - dims = c(length(x = genes.use), max(peak.key)) - ) - rownames(x = coef.matrix) <- genes.use - colnames(x = coef.matrix) <- names(x = peak.key) - links <- LinksToGInteractions( - linkmat = coef.matrix, - gene.coords = gene.coords.use - ) - # add zscores - z.matrix <- sparseMatrix( - i = gene.vec, - j = peak.key[names(x = zscore.vec)], - x = zscore.vec, - dims = c(length(x = genes.use), max(peak.key)) + + 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 ) - rownames(x = z.matrix) <- genes.use - colnames(x = z.matrix) <- names(x = peak.key) - z.lnk <- LinksToGInteractions( - linkmat = z.matrix, - gene.coords = gene.coords.use + + 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 (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 { + sig.df$pvalue <- NA_real_ + } + + 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 (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]] + ) + } + + links <- candidate.links[link.idx] + links$peak <- sig.df$peak + links$gene <- sig.df$gene + + 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)) + } + } + + links$score <- sig.df$score + links$zscore <- sig.df$zscore + links$pvalue <- sig.df$pvalue + + .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), "." + ) + } + + value + } + + # 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$zscore <- z.lnk$score - links$pvalue <- 2 * pnorm(q = -abs(x = links$zscore)) - links <- links[links$pvalue < pvalue_cutoff] + Links(object = object[[peak.assay]], key = key) <- links - return(object) + object } ### Not exported ### diff --git a/R/utilities.R b/R/utilities.R index b16a0e7d..b402b5f0 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -907,32 +907,50 @@ globalVariables( #' @importFrom GenomicRanges makeGRangesFromDataFrame #' @importFrom data.table as.data.table CollapseToLongestTranscript <- function(ranges) { - range.df <- as.data.table(x = ranges) + range.df <- 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 (!"gene_biotype" %in% colnames(x = range.df)) { + range.df[, gene_biotype := if ("gene_type" %in% colnames(x = range.df)) { + gene_type + } else { + NA_character_ + }] + } + range.df$strand <- as.character(x = range.df$strand) range.df$strand <- ifelse( test = range.df$strand == "*", yes = "+", no = range.df$strand ) + collapsed <- range.df[ , list( - unique(seqnames), - min(start), - max(end), - strand[[1]], - gene_biotype[[1]], - gene_name[[1]] + seqnames = unique(seqnames)[1], + start = min(start), + end = max(end), + strand = strand[[1]], + gene_biotype = gene_biotype[[1]], + gene_name = gene_name[[1]] ), - "gene_id" + by = "gene_id" ] - colnames(x = collapsed) <- c( - "gene_id", "seqnames", "start", "end", "strand", "gene_biotype", "gene_name" - ) + collapsed$gene_name <- make.unique(names = collapsed$gene_name) + gene.ranges <- makeGRangesFromDataFrame( df = collapsed, keep.extra.columns = TRUE ) + return(gene.ranges) }