GPlinksR constructs gene-peak regulatory links by combining three sources of
evidence: enhancer-based links, promoter overlaps, and nearest-gene mappings.
The package is designed for ATAC-RNA integration workflows where users already
have a set of peaks and a set of genes, but it also includes a wrapper that can
extract these inputs from common Bioconductor containers such as
SingleCellExperiment (Lun and Risso 2024) and
MultiAssayExperiment (Ramos, Morgan, and Carey 2024).
This vignette uses a small built-in example selected from a KPMP-derived kidney single-cell analysis. It demonstrates how to prepare real peak and gene inputs, run the main function, and use the wrapper when data are already stored in Bioconductor containers.
If you are running these chunks manually from the package source directory
rather than knitting the vignette, load the current source package first with
devtools::load_all(".") so that you use the local development version rather
than an older installed copy.
Install GPlinksR from Bioconductor with BiocManager:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("GPlinksR")The core function in GPlinksR is build_gp_links(). It expects:
pk: peak coordinates in "chr:start-end" format, or a data.frame with
peak coordinate columns.gn: a character vector of gene symbols.For users working with single-cell analysis objects, build_gp_links_wrapper()
can extract these inputs automatically and pass them to build_gp_links().
We begin with the package example dataset, which contains 300 peak coordinates and 100 gene symbols selected from a KPMP-derived analysis.
data("gp_example_inputs", package = "GPlinksR")
pk <- gp_example_inputs$pk
gn <- gp_example_inputs$gn
length(pk)## [1] 300length(gn)## [1] 100For the executable examples below, we use the first three peaks and four genes. Two small files installed with the package contain the corresponding real PEREGRINE version 19 enhancer links and hg38 enhancer coordinates. A local HGNC mapping avoids a BioMart network request while the vignette is built.
pk_demo <- pk[seq_len(3)]
gn_demo <- gn[seq_len(4)]
pk_demo## [1] "chr1:819770-822338" "chr1:876277-877834" "chr1:923554-925121"gn_demo## [1] "TTLL10" "FAM87B" "SAMD11" "FAM41C"The main function is called directly below. This chunk is evaluated during the
vignette build and prints the result produced by build_gp_links().
example_enh_file <- system.file(
"extdata", "gp_example_enhancer_links.tsv",
package = "GPlinksR", mustWork = TRUE
)
example_coords_file <- system.file(
"extdata", "gp_example_enhancer_coordinates.tsv",
package = "GPlinksR", mustWork = TRUE
)
example_gene_map <- data.frame(
hgnc_symbol = c("TTLL10", "SAMD11"),
hgnc_id = c("HGNC:26693", "HGNC:28706")
)
gp <- build_gp_links(
pk = pk_demo,
gn = gn_demo,
enh_file = example_enh_file,
enh_coords_file = example_coords_file,
gene_map = example_gene_map
)
gp## Peak Gene Src
## 1 chr1:819770-822338 TTLL10 enh
## 2 chr1:876277-877834 SAMD11 enh
## 3 chr1:876277-877834 FAM41C prom
## 4 chr1:923554-925121 SAMD11 prom
## 5 chr1:819770-822338 FAM87B clo
## 6 chr1:876277-877834 FAM41C clo
## 7 chr1:923554-925121 SAMD11 cloIf you want to use a previously cached PEREGRINE file, you may first download or retrieve it with:
enh_file <- get_peregrine_file(19)
gp <- build_gp_links(pk = pk_demo, gn = gn_demo, enh_file = enh_file)Some users store peaks as separate chromosome, start, and end columns. This is also supported.
peak_parts <- do.call(rbind, strsplit(pk_demo, "[:-]"))
peak_df <- data.frame(
chr = peak_parts[, 1],
start = as.integer(peak_parts[, 2]),
end = as.integer(peak_parts[, 3])
)
peak_df## chr start end
## 1 chr1 819770 822338
## 2 chr1 876277 877834
## 3 chr1 923554 925121This can be supplied directly to build_gp_links():
gp_from_df <- build_gp_links(
pk = peak_df,
gn = gn_demo,
enh_file = example_enh_file,
enh_coords_file = example_coords_file,
gene_map = example_gene_map
)
identical(gp_from_df, gp)## [1] TRUEFor many Bioconductor workflows, peaks and genes are already stored in a
MultiAssayExperiment. In that setting, the wrapper lets users specify which
experiment contains peaks and which contains genes, and then handles the input
extraction automatically.
The following code shows a simple simulated setup. Here, the RNA experiment is
represented by a SummarizedExperiment (Morgan et al. 2024), while the ATAC
experiment stores peak ranges.
if (requireNamespace("MultiAssayExperiment", quietly = TRUE) &&
requireNamespace("SummarizedExperiment", quietly = TRUE) &&
requireNamespace("IRanges", quietly = TRUE)) {
library(MultiAssayExperiment)
library(SummarizedExperiment)
library(GenomicRanges)
library(IRanges)
library(S4Vectors)
sample_ids <- paste0("cell", seq_len(3))
peak_parts_mae <- do.call(
rbind,
strsplit(pk_demo[seq_len(3)], "[:-]")
)
atac_counts <- matrix(sample.int(10, 9, TRUE), nrow = 3)
colnames(atac_counts) <- sample_ids
atac_se <- SummarizedExperiment(
assays = list(counts = atac_counts),
rowRanges = GRanges(
seqnames = peak_parts_mae[, 1],
ranges = IRanges(
start = as.integer(peak_parts_mae[, 2]),
end = as.integer(peak_parts_mae[, 3])
)
)
)
rna_counts <- matrix(sample.int(10, 12, TRUE), nrow = 4)
colnames(rna_counts) <- sample_ids
rna_se <- SummarizedExperiment(
assays = list(counts = rna_counts),
rowData = DataFrame(
symbol = gn_demo[seq_len(4)]
)
)
mae <- MultiAssayExperiment(
experiments = list(ATAC = atac_se, RNA = rna_se),
colData = S4Vectors::DataFrame(row.names = sample_ids)
)
}Once the object has been created, the wrapper call is concise:
gp_from_mae <- build_gp_links_wrapper(
x = mae,
peak_experiment = "ATAC",
gene_experiment = "RNA",
gene_col = "symbol",
enh_file = example_enh_file,
enh_coords_file = example_coords_file,
gene_map = example_gene_map
)
head(gp_from_mae)## Peak Gene Src
## 1 chr1:819770-822338 TTLL10 enh
## 2 chr1:876277-877834 SAMD11 enh
## 3 chr1:876277-877834 FAM41C prom
## 4 chr1:923554-925121 SAMD11 prom
## 5 chr1:819770-822338 FAM87B clo
## 6 chr1:876277-877834 FAM41C cloThe wrapper also supports a SingleCellExperiment (Lun and Risso 2024)
workflow in which gene-level measurements are stored in the main object and peak
features are stored in an altExp().
if (requireNamespace("SingleCellExperiment", quietly = TRUE) &&
requireNamespace("SummarizedExperiment", quietly = TRUE)) {
library(SingleCellExperiment)
library(SummarizedExperiment)
library(S4Vectors)
sample_ids <- paste0("cell", seq_len(4))
rna_counts <- matrix(sample.int(10, 16, TRUE), nrow = 4)
rownames(rna_counts) <- gn_demo[seq_len(4)]
colnames(rna_counts) <- sample_ids
peak_counts <- matrix(sample.int(10, 12, TRUE), nrow = 3)
colnames(peak_counts) <- sample_ids
peak_rowdata <- DataFrame(
PeakRegion = pk_demo[seq_len(3)]
)
sce <- SingleCellExperiment(
assays = list(counts = rna_counts),
rowData = DataFrame(symbol = rownames(rna_counts))
)
altExp(sce, "ATAC") <- SummarizedExperiment(
assays = list(counts = peak_counts),
rowData = peak_rowdata
)
}The wrapper then extracts genes from the main object and peaks from the named alternative experiment:
gp_from_sce <- build_gp_links_wrapper(
x = sce,
peak_experiment = "ATAC",
gene_col = "symbol",
enh_file = example_enh_file,
enh_coords_file = example_coords_file,
gene_map = example_gene_map
)
head(gp_from_sce)## Peak Gene Src
## 1 chr1:819770-822338 TTLL10 enh
## 2 chr1:876277-877834 SAMD11 enh
## 3 chr1:876277-877834 FAM41C prom
## 4 chr1:923554-925121 SAMD11 prom
## 5 chr1:819770-822338 FAM87B clo
## 6 chr1:876277-877834 FAM41C cloThe wrapper is intended to reduce the amount of manual preprocessing required by new users, but a few object conventions are still important:
"chr:start-end" format.MultiAssayExperiment input, peak and gene experiments should be named
explicitly.SingleCellExperiment input, the main object is assumed to contain genes
and altExp() is assumed to contain peaks.If your object uses different column names, you can point the wrapper to them
with peak_col and gene_col.
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SingleCellExperiment_1.35.2 MultiAssayExperiment_1.39.0
## [3] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [5] GenomicRanges_1.65.1 Seqinfo_1.3.0
## [7] IRanges_2.47.2 S4Vectors_0.51.6
## [9] BiocGenerics_0.59.10 generics_0.1.4
## [11] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [13] GPlinksR_0.99.2 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] KEGGREST_1.53.6 rjson_0.2.23
## [3] xfun_0.60 bslib_0.11.0
## [5] lattice_0.22-9 vctrs_0.7.3
## [7] tools_4.6.1 bitops_1.1-0
## [9] curl_7.1.0 parallel_4.6.1
## [11] AnnotationDbi_1.75.2 tibble_3.3.1
## [13] RSQLite_3.53.3 blob_1.3.0
## [15] pkgconfig_2.0.3 BiocBaseUtils_1.15.1
## [17] Matrix_1.7-6 data.table_1.18.4
## [19] cigarillo_1.3.1 lifecycle_1.0.5
## [21] compiler_4.6.1 Rsamtools_2.29.0
## [23] Biostrings_2.81.6 codetools_0.2-20
## [25] EnsDb.Hsapiens.v86_2.99.0 GenomeInfoDb_1.49.1
## [27] htmltools_0.5.9 sass_0.4.10
## [29] lazyeval_0.2.3 RCurl_1.98-1.19
## [31] yaml_2.3.12 pillar_1.11.1
## [33] crayon_1.5.3 jquerylib_0.1.4
## [35] BiocParallel_1.47.0 DelayedArray_0.39.3
## [37] cachem_1.1.0 abind_1.4-8
## [39] tidyselect_1.2.1 digest_0.6.39
## [41] dplyr_1.2.1 restfulr_0.0.17
## [43] bookdown_0.47 grid_4.6.1
## [45] fastmap_1.2.0 SparseArray_1.13.2
## [47] cli_3.6.6 magrittr_2.0.5
## [49] GenomicFeatures_1.65.0 S4Arrays_1.13.0
## [51] XML_3.99-0.23 UCSC.utils_1.9.0
## [53] bit64_4.8.2 rmarkdown_2.31
## [55] XVector_0.53.0 httr_1.4.8
## [57] bit_4.6.0 otel_0.2.0
## [59] png_0.1-9 memoise_2.0.1
## [61] evaluate_1.0.5 knitr_1.51
## [63] BiocIO_1.23.3 rtracklayer_1.73.0
## [65] rlang_1.3.0 glue_1.8.1
## [67] DBI_1.3.0 ensembldb_2.37.3
## [69] BiocManager_1.30.27 jsonlite_2.0.0
## [71] AnnotationFilter_1.37.0 R6_2.6.1
## [73] ProtGenerics_1.45.0 GenomicAlignments_1.49.1Lun, Aaron, and Davide Risso. 2024. SingleCellExperiment: S4 Classes for Single Cell Data.
Morgan, Martin, Valerie Obenchain, Jim Hester, and Hervé Pagès. 2024. SummarizedExperiment: SummarizedExperiment Container.
Ramos, Marcel, Mike Morgan, and Vince Carey. 2024. MultiAssayExperiment: Software for Integrative Analysis of Multiomics Experiments in Bioconductor.