---
title: "GPlinksR: Building Gene-Peak Links from Example Single-Cell Inputs"
author: "Xinran Wang, Kelly Street, Huaiyu Mi, and Bryan Queme"
date: "`r Sys.Date()`"
bibliography: bibFile.bib
output:
    BiocStyle::html_document:
        toc: true
vignette: >
    %\VignetteIndexEntry{GPlinksR: Gene-Peak Links from Example Inputs}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    echo = TRUE,
    eval = TRUE,
    message = FALSE,
    warning = FALSE
)
set.seed(123)
library(GPlinksR)
```

# Introduction

`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` [@SingleCellExperiment] and
`MultiAssayExperiment` [@MultiAssayExperiment].

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.

# Installation

Install `GPlinksR` from Bioconductor with `BiocManager`:

```{r installation, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("GPlinksR")
```

# Overview

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()`.

# Minimal Example with Real Peak and Gene Vectors

We begin with the package example dataset, which contains 300 peak coordinates
and 100 gene symbols selected from a KPMP-derived analysis.

```{r example-inputs}
data("gp_example_inputs", package = "GPlinksR")

pk <- gp_example_inputs$pk
gn <- gp_example_inputs$gn

length(pk)
length(gn)
```

For 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.

```{r demo-subset}
pk_demo <- pk[seq_len(3)]
gn_demo <- gn[seq_len(4)]

pk_demo
gn_demo
```

The main function is called directly below. This chunk is evaluated during the
vignette build and prints the result produced by `build_gp_links()`.

```{r direct-run}
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
```

If you want to use a previously cached PEREGRINE file, you may first download
or retrieve it with:

```{r cached-enhancer-file, eval=FALSE}
enh_file <- get_peregrine_file(19)
gp <- build_gp_links(pk = pk_demo, gn = gn_demo, enh_file = enh_file)
```

# Using a Peak Data Frame

Some users store peaks as separate chromosome, start, and end columns. This is
also supported.

```{r peak-data-frame}
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
```

This can be supplied directly to `build_gp_links()`:

```{r data-frame-run}
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)
```

# Wrapper Workflow for MultiAssayExperiment

For 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` [@SummarizedExperiment], while the ATAC
experiment stores peak ranges.

```{r mae-setup}
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:

```{r mae-run}
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)
```

# Wrapper Workflow for SingleCellExperiment

The wrapper also supports a `SingleCellExperiment` [@SingleCellExperiment]
workflow in which gene-level measurements are stored in the main object and peak
features are stored in an `altExp()`.

```{r sce-setup}
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:

```{r sce-run}
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)
```

# Input Expectations

The wrapper is intended to reduce the amount of manual preprocessing required by
new users, but a few object conventions are still important:

* Peak features should ultimately resolve to coordinates in
    `"chr:start-end"` format.
* Gene features should ultimately resolve to gene symbols.
* For `MultiAssayExperiment` input, peak and gene experiments should be named
    explicitly.
* For `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`.

# Session Info

```{r session}
sessionInfo()
```