---
title: "GRIN2: Genomic Random Interval Analysis"
author:
  - "Abdelrahman Elsayed, PhD"
  - "Stanley Pounds, PhD"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{GRIN2: Genomic Random Interval Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
library(GRIN2)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE,
  error = FALSE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  fig.width = 8,
  fig.height = 6
)

options(width = 80, digits = 3)

data(
  list = c(
    "clin_data",
    "lesion_data",
    "expr_data",
    "hg38_gene_annotation",
    "hg38_chrom_size",
    "hg38_cytoband",
    "pathways",
    "grin.results",
    "example_exon_annotation",
    "hg38_exon_chrom_size"
  ),
  package = "GRIN2"
)
```


```{r table-display-helper, include=FALSE}
show.results.table <- function(data, columns, order.by = NULL, n = 6,
                               caption = NULL, digits = 2) {
  # Optionally order the table and exclude missing ordering values
  if (!is.null(order.by) && order.by %in% names(data)) {
    order.values <- as.numeric(as.character(data[[order.by]]))
    data <- data[
      order(order.values, na.last = NA),
      ,
      drop = FALSE
    ]
  }

  # Retain the requested columns
  columns <- intersect(columns, names(data))
  display <- head(data[, columns, drop = FALSE], n)

  # Identify p- and q-value columns
  significance.columns <- grep(
    "(^[pq][0-9]*\\.)|(_[pq]val(\\.adj)?$)",
    names(display),
    value = TRUE
  )

  # Use scientific notation without changing the original results
  display[significance.columns] <- lapply(
    display[significance.columns],
    function(x) {
      formatC(
        as.numeric(as.character(x)),
        format = "e",
        digits = digits
      )
    }
  )

  knitr::kable(
    display,
    row.names = FALSE,
    caption = caption
  )
}
```

## Overview

The **GRIN2** package implements the Genomic Random Interval (GRIN)
framework for identifying genes and genomic loci affected by genomic lesions
more often than expected by chance. The package also provides tools for:

- standard gene- and chromosome-level GRIN analysis;
- exon-level target-size modeling for selected lesion types;
- lesion constellation analysis;
- genome-wide and regional lesion visualization;
- association analysis between genomic lesions and gene expression;
- association analysis between genomic lesions and clinical outcomes; and
- association analysis between gene expression and clinical outcomes.

This vignette presents an end-to-end workflow using only datasets included with GRIN2. Every code chunk is evaluated when the vignette is built. Internet-dependent operations, including downloading annotation bundles and retrieving transcript databases, are intentionally not run in this vignette, allowing it to be built reproducibly and offline.

## Example data

The bundled example data are derived from a previously published T-cell acute
lymphoblastic leukemia (T-ALL) cohort containing RNA-sequencing, whole-exome
sequencing, genomic lesion, and clinical outcome data. Additional information
about the study is available in Liu et al. (2017),
[The genomic landscape of pediatric and young adult T-lineage acute lymphoblastic leukemia](https://doi.org/10.1038/ng.3909).

To keep the package compact and the vignette quick to build, the bundled
gene-annotation and expression datasets are reduced example datasets containing
417 genes. They are intended for demonstrating the GRIN2 workflow and should
not be interpreted as complete genome-wide datasets.

The datasets used in this vignette are summarized below.

```{r bundled-data}
data.summary <- data.frame(
  object = c(
    "lesion_data",
    "expr_data",
    "clin_data",
    "hg38_gene_annotation",
    "hg38_chrom_size",
    "hg38_cytoband",
    "pathways",
    "grin.results",
    "example_exon_annotation",
    "hg38_exon_chrom_size"
  ),
  purpose = c(
    "GRCh38 genomic lesions by subject and lesion type",
    "Gene-by-subject expression matrix",
    "Subject-level clinical and outcome data",
    "Example GRCh38 gene annotation",
    "GRCh38 chromosome lengths",
    "GRCh38 cytoband annotation",
    "Example pathway definitions",
    "Precomputed GRIN result object for plotting examples",
    "Example exon coordinates for exon-level GRIN analysis",
    "Chromosome-level exonic target sizes for GRCh38"
  )
)

knitr::kable(data.summary, row.names = FALSE)
```

All genomic coordinates supplied to GRIN2 must use the same genome assembly.
The current GRIN2 annotation workflow supports the human GRCh38 assembly.
Coordinates from another assembly should be converted to GRCh38 before the
analysis.

### Required input structure

GRIN2 uses specific column names to connect lesion, expression, annotation,
and clinical data. The lesion dataset must contain the five required columns
`ID`, `chrom`, `loc.start`, `loc.end`, and `hit.type`. The `ID` column identifies
the subject, `chrom` specifies the chromosome, `loc.start` and `loc.end` define
the genomic interval, and `hit.type` identifies the lesion type.

In the expression dataset, the first column must be named `gene` and contain
unique, unversioned Ensembl gene identifiers. The remaining columns contain
expression measurements for individual subjects. The clinical dataset must
contain an `ID` column identifying each subject. These identifiers are used to
match subjects across the lesion, expression, and clinical datasets.

The examples below show the required structure of each input dataset.

```{r inspect-inputs}
# Selected lesion records representing different lesion types
lesion.rows <- c(1, 2, 4849, 5066, 6239, 6854)

knitr::kable(
  lesion_data[lesion.rows, , drop = FALSE],
  row.names = FALSE,
  caption = "Example lesion data"
)

knitr::kable(
  head(expr_data[, seq_len(min(6, ncol(expr_data))), drop = FALSE]),
  row.names = FALSE,
  caption = "Example expression data"
)

knitr::kable(
  head(clin_data),
  row.names = FALSE,
  caption = "Example clinical data"
)
```

## Genome annotations

The bundled `hg38_gene_annotation` and `hg38_chrom_size` objects provide a
compact, reproducible dataset for the examples in this vignette.

```{r inspect-annotations}
selected.genes <- c("RPL5", "NRAS", "CDKN2A", "IKZF5", "WT1", "EZH2")

annotation.example <- hg38_gene_annotation[
  match(selected.genes, hg38_gene_annotation$gene.name),
  ,
  drop = FALSE
]

knitr::kable(
  annotation.example,
  row.names = FALSE,
  caption = "Example gene annotation data"
)

knitr::kable(head(hg38_chrom_size), row.names = FALSE)
```

For a complete analysis, `get.ensembl.annotation()` can retrieve versioned
GRCh38 Ensembl gene, exon, and regulatory-element annotation bundles. Valid
downloads are cached and reused, and MD5 checksums are used to verify cached
files. Because annotation retrieval requires network access, it is not run in
this vignette; see `?get.ensembl.annotation` for the current interface.

## GRIN analysis

### Standard gene-level analysis

The standard GRIN model uses the genomic length of each gene together with the
length of its chromosome to calculate lesion probabilities.

```{r standard-grin}
standard.grin.results <- grin.stats(
  lsn.data = lesion_data,
  gene.data = hg38_gene_annotation,
  chr.size = hg38_chrom_size
)
```

### Exon-level target-size analysis

For protein-altering variants that can occur only in coding regions, a whole
gene and whole chromosome may be inappropriate target spaces. Examples include
missense, nonsense, and coding frameshift variants, particularly in studies
based on whole-exome sequencing.

GRIN2 can use the annotated exon length of each gene and the total annotated
exonic target size of the corresponding chromosome for selected lesion types.
Lesion types not listed in `exon_level` continue to use conventional gene and
chromosome target sizes in the same analysis.

The following analysis uses exon-level target sizes for mutations while
retaining the standard model for all other lesion types.

```{r exon-level-grin}
grin.exon.results <- grin.stats(
  lsn.data = lesion_data,
  gene.data = hg38_gene_annotation,
  chr.size = hg38_chrom_size,
  exons.annotation = example_exon_annotation,
  exon.chrom.size = hg38_exon_chrom_size,
  exon_level = "mutation"
)

grin.exon.results$exon_level
head(grin.exon.results$gene.exon.size)
knitr::kable(head(grin.exon.results$exon.chrom.size), row.names = FALSE)
```

The values supplied to `exon_level` must match the lesion-type labels in
`lesion_data`. Versioned exon annotations containing one selected transcript
per gene can be retrieved using
`get.ensembl.annotation(annotation.type = "exon")`. When supplying custom exon
annotations, use one selected transcript per gene so that alternative or
overlapping exons are not counted repeatedly.

### Inspecting GRIN results

The primary statistical results are stored in the `gene.hits` component.

```{r inspect-grin-results}
grin.table <- standard.grin.results$gene.hits

show.grin.columns <- function(pattern, max.columns = 8) {
  annotation.columns <- intersect(
    c("gene", "gene.name", "chrom", "loc.start", "loc.end"),
    names(grin.table)
  )

  result.columns <- grep(
    pattern,
    names(grin.table),
    value = TRUE
  )

  selected.columns <- unique(c(
    annotation.columns,
    head(result.columns, max.columns)
  ))

  order.column <- if ("p2.nsubj" %in% names(grin.table)) {
    "p2.nsubj"
  } else {
    NULL
  }

  show.results.table(
    data = grin.table,
    columns = selected.columns,
    order.by = order.column
  )
}
```

The subject-count columns report how many unique subjects are affected by each
lesion type.

```{r subject-count-results}
show.grin.columns("^nsubj\\.")
```

Lesion-specific probability and false-discovery rate columns evaluate whether
each gene is affected more often than expected for an individual lesion type.

```{r lesion-probability-results}
show.grin.columns("^[pq]\\.nsubj\\.")
```

Constellation columns evaluate whether a gene is affected by at least one,
two, three, or more distinct lesion types.

```{r constellation-results}
show.grin.columns("^[pq][0-9]+\\.nsubj$")
```

Corresponding hit-level columns count all lesions, whereas subject-level
columns count an affected subject once for the relevant lesion category.

## Exporting results

`write.grin.xlsx()` exports the GRIN analysis results and supporting information
to a multi-sheet Excel workbook. The workbook contains the following sheets:

- `gene.hits`: the main GRIN results, including lesion counts, numbers of
  affected subjects, enrichment p-values, and FDR-adjusted q-values.
- `lsn.data`: the lesion data used in the analysis.
- `gene.data`: the gene annotation data used in the analysis.
- `chr.size`: the chromosome sizes used to calculate lesion-hit probabilities.
- `interpretation`: descriptions of the results and output columns.
- `method.paragraph`: a summary of the GRIN methodology and relevant references.

Beginning with GRIN2 version 2.1.0, the potentially very large `gene.lsn.data`
table is excluded from the workbook to reduce file size and improve export
performance.

```{r export-results}
output.file <- file.path(tempdir(), "GRIN2_example_results.xlsx")

write.grin.xlsx(
  grin.result = standard.grin.results,
  output.file = output.file
)

stopifnot(file.exists(output.file))
unlink(output.file)
```

## Visualizing genomic lesions

### Genome-wide lesion plot

```{r genomewide-lesion-plot, fig.height=7, fig.width=8, fig.cap="Genome-wide distribution and statistical significance of genomic lesions."}
genomewide.lsn.plot(
  standard.grin.results,
  max.log10q = 50
)
```

### Stacked lesion-count plot

```{r stacked-barplot, fig.height=8, fig.width=9, fig.cap="Numbers of subjects affected by different lesion types in selected genes."}
genes.of.interest <- c(
  "CDKN2A", "NOTCH1", "CDKN2B", "TAL1", "FBXW7", "PTEN", "IRF8",
  "NRAS", "BCL11B", "MYB", "LEF1", "RB1", "MLLT3", "EZH2", "ETV6",
  "CTCF", "JAK1", "KRAS", "RUNX1", "IKZF1", "KMT2A", "RPL11",
  "TCF7", "WT1", "JAK2", "JAK3", "FLT3"
)

grin.barplt(
  standard.grin.results,
  genes.of.interest
)
```

### Regional lesion plot

The bundled `grin.results` object allows regional plotting without rerunning
the analysis. The example below uses only bundled data and disables the
transcript and ideogram tracks, avoiding any external annotation lookup.

```{r chr9-lesion-plot, fig.keep="last", fig.height=8, fig.width=9, fig.cap="Distribution of all lesion types across chromosome 9 without transcript annotations or a chromosome ideogram."}
lsn.transcripts.plot(
  grin.res = grin.results,
  chrom = 9,
  plot.start = 1,
  plot.end = 138394717,
  transTrack = FALSE,
  show.ideogram = FALSE,
  point.size.mm = 1
)
```

Gene-centered transcript plots and transcript selection are supported by
`lsn.transcripts.plot()`, but require an external `EnsDb` transcript annotation
object. Those network-dependent examples are documented on the function's help
page rather than included in this offline vignette.

### Preparing an OncoPrint matrix

`grin.oncoprint.mtx()` converts GRIN results into a gene-by-subject lesion
matrix compatible with `ComplexHeatmap::oncoPrint()`. Matrix rows follow the
gene order requested by the user. Gene symbols are used as row labels when
available; otherwise, Ensembl gene IDs are used. The function also ensures
that all output gene labels are unique.

```{r prepare-oncoprint}
oncoprint.genes <- c(
  "ENSG00000101307", "ENSG00000171862", "ENSG00000138795",
  "ENSG00000139083", "ENSG00000162434", "ENSG00000134371",
  "ENSG00000118058", "ENSG00000171843", "ENSG00000139687",
  "ENSG00000184674"
)

oncoprint.mtx <- grin.oncoprint.mtx(
  grin.results,
  oncoprint.genes
)

dim(oncoprint.mtx)
oncoprint.mtx[, seq_len(min(6, ncol(oncoprint.mtx))), drop = FALSE]

onco.props <- onco.print.props(
  lesion_data,
  hgt = c(gain = 5, loss = 4, mutation = 2, fusion = 1)
)

str(onco.props, max.level = 1)
```

## Gene-by-subject lesion matrices

GRIN2 provides two complementary lesion-matrix formats. Both begin by
preparing gene and lesion coordinates and identifying gene-lesion overlaps.

```{r prepare-gene-lesion-overlaps}
gene.lsn <- prep.gene.lsn.data(
  lsn.data = lesion_data,
  gene.data = hg38_gene_annotation
)

gene.lsn.overlap <- find.gene.lsn.overlaps(gene.lsn)
```

### Lesion-group matrix

`prep.lsn.type.matrix()` returns one row per gene and one column per subject.
Each entry identifies the subject's lesion group for that gene. Subjects with
more than one lesion type affecting the same gene are assigned to the
`"multiple"` group. This categorical representation is used when subjects are
compared across gene-specific lesion groups, including comparisons of survival
distributions with `grin.logRank()` and comparisons of gene expression across
lesion groups in the lesion-expression workflow.

```{r lesion-group-matrix}
gene.lsn.type.mtx <- prep.lsn.type.matrix(
  gene.lsn.overlap,
  min.ngrp = 5
)

dim(gene.lsn.type.mtx)
gene.lsn.type.mtx[
  seq_len(min(6, nrow(gene.lsn.type.mtx))),
  seq_len(min(6, ncol(gene.lsn.type.mtx))),
  drop = FALSE
]
```

### Binary lesion matrix

`prep.binary.lsn.mtx()` creates a separate binary row for each gene-lesion
combination, such as `NOTCH1_mutation`. An entry of 1 indicates that the subject
has the specified lesion, whereas 0 indicates that the lesion is absent. This
representation supports presence-versus-absence analyses of individual lesion
types and is used by `grin.assoc.lsn.outcome()` to evaluate associations between
specific gene-lesion combinations and clinical outcomes.

```{r binary-lesion-matrix}
lsn.binary.mtx <- prep.binary.lsn.mtx(
  gene.lsn.overlap,
  min.ngrp = 5
)

dim(lsn.binary.mtx)
lsn.binary.mtx[
  seq_len(min(6, nrow(lsn.binary.mtx))),
  seq_len(min(6, ncol(lsn.binary.mtx))),
  drop = FALSE
]
```

## Association between lesions and gene expression

### Preparing matched lesion and expression data

`alex.prep.lsn.expr()` matches the lesion and expression data by gene and
subject, then arranges both datasets so that corresponding rows represent the
same genes and corresponding columns represent the same subjects. The returned
object is used as input to `KW.hit.express()`, which tests whether gene
expression differs among genomic lesion groups.

```{r prepare-alex-data}
alex.data <- alex.prep.lsn.expr(
  expr_data,
  lesion_data,
  hg38_gene_annotation,
  min.expr = 1,
  min.pts.lsn = 5
)

dim(alex.data$alex.lsn)
dim(alex.data$alex.expr)
```

### Kruskal-Wallis association analysis

`KW.hit.express()` uses the Kruskal–Wallis test to determine whether gene
expression differs among the observed lesion groups. The `min.grp.size`
argument specifies the minimum number of subjects required in each group for
that group to be included in the test, helping to avoid unstable comparisons
based on sparsely represented lesion groups.

```{r alex-kruskal-wallis}
alex.kw.results <- KW.hit.express(
  alex.data,
  hg38_gene_annotation,
  min.grp.size = 5
)

# select columns to display in the result table
show.results.table(
  data = alex.kw.results,
  columns = c(
    "gene",
    "gene.name",
    "p.KW",
    "q.KW"
  ),
  order.by = "q.KW",
  caption = "Genes with the smallest Kruskal–Wallis q-values"
)
```

### Waterfall plot

Waterfall plots provide a side-by-side representation of genomic lesions and
gene expression across subjects for a selected gene. Subjects are ordered first
by lesion group and then, within each group, by the expression level of the
selected gene. This arrangement facilitates visual comparison of expression
patterns across the different lesion groups.

```{r wt1-waterfall-plot, fig.height=6, fig.width=7, fig.cap="WT1 expression and lesion groups across subjects."}
WT1.waterfall.data <- alex.waterfall.prep(
  alex.data,
  alex.kw.results,
  "WT1",
  lesion_data
)

alex.waterfall.plot(
  WT1.waterfall.data,
  lesion_data
)
```

## Association with clinical outcomes

Create survival outcomes once and reuse the resulting clinical data for the
three outcome-association workflows.

```{r prepare-clinical-outcomes}
clinical <- clin_data

clinical$EFS <- survival::Surv(
  clinical$efs.time,
  clinical$efs.censor
)

clinical$OS <- survival::Surv(
  clinical$os.time,
  clinical$os.censor
)
```

### Associations between individual lesions and clinical outcomes

`grin.assoc.lsn.outcome()` evaluates whether the presence or absence of each
gene-lesion combination is associated with a clinical outcome. For a binary
outcome, such as response versus no response, the function uses logistic
regression to compare the odds of the outcome between subjects with and without
the lesion. For a time-to-event outcome, such as event-free survival, it uses a
Cox proportional hazards model to evaluate whether the lesion is associated
with the event hazard while accounting for follow-up time.

```{r lesion-outcome-association}
lesion.outcome.results <- grin.assoc.lsn.outcome(
  lsn.mtx = lsn.binary.mtx,
  clin.data = clinical,
  annotation.data = hg38_gene_annotation,
  clinvars = c("MRD.binary", "EFS")
)

# select columns to display in the result table
show.results.table(
  data = lesion.outcome.results,
  columns = c(
    "Gene_lsn",
    "gene.name",
    "cox_EFS_pval",
    "cox_EFS_qval",
    "logistic_MRD.binary_pval",
    "logistic_MRD.binary_qval"
  ),
  order.by = "cox_EFS_pval",
  caption = paste(
    "Gene-lesion associations ordered by EFS Cox-model p-value,",
    "with corresponding MRD logistic-regression results"
  )
)
```

### Lesion-group log-rank analysis

`grin.logRank()` uses the log-rank test to determine whether time-to-event
distributions differ among the lesion groups observed for each gene. The
`min.grp.size` argument specifies the minimum total number of evaluable subjects
required in each retained lesion group; both subjects with events and censored
subjects contribute to this number. In addition to the log-rank p-values and
FDR-adjusted q-values, the returned results report the numbers of subjects with
and without events for each lesion group.

```{r lesion-group-logrank}
logrank.results <- grin.logRank(
  lsn.mtx = gene.lsn.type.mtx,
  clin.data = clinical,
  annotation.data = hg38_gene_annotation,
  clinvars = "EFS",
  min.grp.size = 4
)

# select columns to display in the result table 
show.results.table(
  data = logrank.results,
  columns = c(
    "gene",
    "gene.name",
    "logRank_EFS_pval",
    "logRank_EFS_qval"
  ),
  order.by = "logRank_EFS_pval",
  caption = "Genes with the smallest EFS log-rank p-values"
)
```

### Gene-expression outcome associations

`grin.assoc.expr.outcome()` performs gene-level outcome association analyses
using a gene-by-subject expression matrix. Cox proportional hazards models are
used for time-to-event outcomes, and logistic regression models are used for
binary outcomes. Optional covariates are included in each gene-level model.

The first column of `expr.mtx` must be named `gene` and contain unique,
unversioned Ensembl gene IDs. The remaining columns must contain normalized
expression values, with subject IDs as column names.

```{r expression-outcome-association}
expression.outcome.results <- grin.assoc.expr.outcome(
  expr.mtx = expr_data,
  clin.data = clinical,
  annotation.data = hg38_gene_annotation,
  clinvars = c("MRD.binary", "EFS"),
  covariate = "WBC"
)

# Extract columns to display in the results table 
show.results.table(
  data = expression.outcome.results,
  columns = c(
    "gene",
    "gene.name",
    "logistic_MRD.binary_pval.adj",
    "logistic_MRD.binary_qval.adj",
    "cox_EFS_pval_adj",
    "cox_EFS_qval_adj"
  ),
  order.by = "cox_EFS_pval_adj",
  caption = paste(
    "Gene-expression associations ordered by the adjusted EFS Cox-model",
    "p-value, with corresponding adjusted MRD logistic-regression results"
  )
)

```

## Summary

This vignette demonstrated a complete GRIN2 workflow using bundled data:

1. running conventional and exon-level GRIN analyses;
2. inspecting, exporting, and visualizing GRIN results;
3. preparing categorical and binary gene-lesion matrices;
4. associating genomic lesions with gene expression;
5. testing lesion-group associations with survival outcomes using
   `grin.logRank()`;
6. testing binary lesion associations with binary and survival outcomes using
   `grin.assoc.lsn.outcome()`; and
7. testing gene-expression associations with binary and survival outcomes
   using `grin.assoc.expr.outcome()`.

For complete argument descriptions and additional function-specific details,
consult the individual GRIN2 help pages.

## References

- Pounds S, et al. (2013). A genomic random interval model for statistical
  analysis of genomic lesion data. *Bioinformatics*, 29(17), 2088–2095.
  <https://doi.org/10.1093/bioinformatics/btt372>
- Cao X, Elsayed AH, Pounds SB. (2023). Statistical methods inspired by
  challenges in pediatric cancer multi-omics. In Fridley B and Wang X (Eds.),
  *Statistical Genomics*. Methods in Molecular Biology, vol. 2629,
  pp. 349–373. Humana, New York, NY.
  <https://doi.org/10.1007/978-1-0716-2986-4_16>
- Pounds S, Cheng C. (2006). Robust estimation of the false discovery rate.
  *Bioinformatics*, 22(16), 1979–1987.
  <https://doi.org/10.1093/bioinformatics/btl328>
- Liu Y, et al. (2017). The genomic landscape of pediatric and young adult
  T-lineage acute lymphoblastic leukemia. *Nature Genetics*, 49(8), 1211–1218.
  <https://doi.org/10.1038/ng.3909>