---
title: "perturbMatch reference signatures via ExperimentHub"
author:
- name: Jamie Soul
  affiliation: University of Liverpool
  email: jamie.soul@liverpool.ac.uk
package: perturbMatchData
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    number_sections: true
abstract: >
  perturbMatchData is an ExperimentHub package that distributes the reference
  perturbation expression signatures used by the perturbMatch package for
  upstream regulator analysis. This vignette describes the available datasets
  and shows how to discover, download and load them.
vignette: >
  %\VignetteIndexEntry{perturbMatch reference signatures via ExperimentHub}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```

# Introduction

`perturbMatchData` is an
[ExperimentHub](https://bioconductor.org/packages/ExperimentHub/) package that
provides a collection of human reference perturbation expression signatures as
HDF5-backed `SummarizedExperiment` objects. Each signature summarises the
differential expression observed after a single-gene perturbation (knockdown,
knockout or overexpression) in an RNA-seq experiment.

The signatures were assembled with an automated Nextflow pipeline and
summarised in two ways: using the **characteristic direction** method and
using **limma** (log2 fold-change together with a *pi* score,
`log2FC * -log10(padj)`). Both are provided with and without SVA batch
correction, giving four datasets in total.

These reference databases are used by the
[perturbMatch](https://github.com/soulj/perturbMatch) package to
identify the upstream regulators most likely to drive a query gene expression
signature, and can also be used directly to ask what happens to gene expression
when a given gene is perturbed.

# Installation

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

```{r load-pkg}
library(perturbMatchData)
library(SummarizedExperiment)
```

```{r hub-ready, include = FALSE}
## The records are registered with ExperimentHub once the package is accepted,
## so the chunks that retrieve them do nothing while the hub has none
hub_ready <- tryCatch(
    length(perturbMatch_reference()) > 0L,
    error = function(e) FALSE
)
```

# Available datasets

The four datasets are:

| Accessor                    | Method                    | Batch correction |
|-----------------------------|---------------------------|------------------|
| `perturbMatch_chrdir()`     | Characteristic direction  | No               |
| `perturbMatch_chrdir_sva()` | Characteristic direction  | SVA              |
| `perturbMatch_limma()`      | limma (log2FC / pi score) | No               |
| `perturbMatch_limma_sva()`  | limma (log2FC / pi score) | SVA              |

# Discovering the data on ExperimentHub

Most work starts from the accessors in the previous section, which handle
retrieval for you. The underlying records can also be browsed like any other
ExperimentHub resource, which is useful for checking what is hosted and what is
already in the local cache. `perturbMatch_reference()` returns the records
prepared by this package:

```{r discover, eval = hub_ready}
library(ExperimentHub)

perturbMatch_reference()

## querying the hub directly also matches on the other metadata fields
eh <- ExperimentHub()
query(eh, "perturbMatchData")
```

Each dataset is stored as a record of signature matrices and a record of
signature annotation, alongside one gene annotation record shared by all four.
Every record carries a stable ExperimentHub accession, so a dataset can also be
composed by hand from those accessions.

# Loading a dataset

Use the named accessors to assemble a dataset as a `SummarizedExperiment`.
Each resource is downloaded once and subsequent calls read from the local
ExperimentHub cache. The signature matrices are large, 790 MB for `chrdir` and
673 MB for `chrdir_sva`, and 3.1 GB and 2.6 GB for `limma` and `limma_sva`,
which carry three assays each. The following chunk is therefore shown rather
than run when this vignette is built:

```{r load, eval = FALSE}
se <- perturbMatch_chrdir()
se
```

To see which perturbations are available without retrieving those matrices, ask
for the annotation alone. This is what the rest of this vignette uses:

```{r load-skeleton, eval = hub_ready}
se <- perturbMatch_chrdir(withAssays = FALSE)
se
```

The accessors wrap the standard ExperimentHub interface, so the resources can
also be retrieved by accession and composed by hand. Look the accessions up by
title rather than writing them out, so that they follow the deposit:

```{r load-by-id, eval = hub_ready}
recs <- perturbMatch_reference()
ids <- setNames(names(recs), recs$title)
ids
```

What distinguishes the four datasets, including which HDF5 datasets the assay
resource holds and the scoring metric, is attached to the `colData` resource,
so composing an object by hand means carrying it across:

```{r compose-by-hand, eval = FALSE}
eh <- ExperimentHub()
h5 <- eh[[ids[["perturbMatch_chrdir_assays"]]]]
colData <- eh[[ids[["perturbMatch_chrdir_colData"]]]]

se <- SummarizedExperiment(
    assays = list(chrdir = HDF5Array::HDF5Array(h5, "chrdir")),
    rowData = eh[[ids[["perturbMatch_rowData"]]]],
    colData = colData
)
metadata(se) <- metadata(colData)
```

# Exploring the data

Each object is a `SummarizedExperiment` with genes in the rows and
perturbation signatures in the columns. The assays are HDF5-backed delayed
arrays, so the data is read from disk on demand rather than held in memory.

```{r explore-annotation, eval = hub_ready}
## Signature-level metadata (one row per perturbation)
head(colData(se))

## Gene annotation data from NCBI
head(rowData(se))

## How the scoring matrix should be derived from the assay(s)
metadata(se)$ref_metric
```

The signature matrix itself is only present when the object was built with
`withAssays = TRUE`, so this chunk is again shown rather than run:

```{r explore-assay, eval = FALSE}
assay(se)[1:10, 1:3]
```

The `limma` datasets carry three assays (`log2fc`, `padj` and the combined
`pi` score) and the `chrdir` datasets a single `chrdir` assay. The
`metadata(se)$ref_metric` field records which assays downstream tools such as
`perturbMatch::calcEnrichment()` should use.

# Data provenance

The scripts used to build these objects from the signature results are
included with the package.

```{r provenance}
list.files(system.file("scripts", package = "perturbMatchData"))
```

`run_nextflow.sh` records the
[LLM-GeoSig](https://github.com/soulj/LLM-GeoSig) run that annotated the GEO
series and computed the signatures from the GEO RNA-seq counts matrices,
`make-data.R` assembles each `SummarizedExperiment` and `make-metadata.R`
writes the ExperimentHub metadata.

The case and control groupings, the perturbed gene and the perturbation type
were assigned by large language models rather than curated by hand, using the
agreement of two models and a curated subset for calibration. Treat the
`colData` annotation as automated annotation: it is accurate enough to rank
regulators across thousands of signatures, but any individual signature should
be checked against its GEO series before it is relied on.

The assembled resources are deposited on Zenodo under the concept DOI
[10.5281/zenodo.21792179](https://doi.org/10.5281/zenodo.21792179), which
resolves to the current version of the deposit. The ExperimentHub records point
at a fixed versioned record, so the files retrieved by a given release of this
package do not change.

# Citing the data

Please cite the Zenodo deposit above alongside the methods the signatures were
computed with.

The source counts are the NCBI-generated RNA-seq counts matrices for each GEO
series, described at
<https://www.ncbi.nlm.nih.gov/geo/info/rnaseqcounts.html>.

Clark NR, Hu KS, Feldmann AS, Kou Y, Chen EY, Duan Q, Ma'ayan A. The
characteristic direction: a geometrical approach to identify differentially
expressed genes. *BMC Bioinformatics*. 2014;15:79.
doi:[10.1186/1471-2105-15-79](https://doi.org/10.1186/1471-2105-15-79)

Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers
differential expression analyses for RNA-sequencing and microarray studies.
*Nucleic Acids Research*. 2015;43(7):e47.
doi:[10.1093/nar/gkv007](https://doi.org/10.1093/nar/gkv007)

Xiao Y, Hsiao TH, Suresh U, Chen HIH, Wu X, Wolf SE, Chen Y. A novel
significance score for gene selection and ranking. *Bioinformatics*.
2014;30(6):801-807. doi:
[10.1093/bioinformatics/btr671](https://doi.org/10.1093/bioinformatics/btr671)

Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by
surrogate variable analysis. *PLoS Genetics*. 2007;3(9):e161.
doi:[10.1371/journal.pgen.0030161](https://doi.org/10.1371/journal.pgen.0030161)

# Session info

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