# GRIN2: Genomic Random Interval Analysis

#### Abdelrahman Elsayed, PhD

#### Stanley Pounds, PhD

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

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.

    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)

object | purpose
--- | ---
lesion_data | GRCh38 genomic lesions by subject and lesion type
expr_data | Gene-by-subject expression matrix
clin_data | Subject-level clinical and outcome data
hg38_gene_annotation | Example GRCh38 gene annotation
hg38_chrom_size | GRCh38 chromosome lengths
hg38_cytoband | GRCh38 cytoband annotation
pathways | Example pathway definitions
grin.results | Precomputed GRIN result object for plotting examples
example_exon_annotation | Example exon coordinates for exon-level GRIN analysis
hg38_exon_chrom_size | Chromosome-level exonic target sizes for GRCh38

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.

    # 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"
    )

Example lesion data
ID | chrom | loc.start | loc.end | lsn.type
--- | --- | --- | --- | ---
PARFIH | 16 | 67616879 | 67616879 | mutation
PARFIH | 3 | 49722238 | 49722238 | mutation
PARHBT | 8 | 54461421 | 70043021 | gain
PASZJW | 9 | 21993125 | 22103925 | loss
PATWIJ | 21 | 32859097 | 32886028 | loss
PATHGY | 11 | 47360041 | 47360041 | fusion

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

Example expression data
gene | PARFIH | PARFPJ | PARFXJ | PARKLK | PARLJA
--- | --- | --- | --- | --- | ---
ENSG00000177076 | 1.80 | 0.473 | 0.928 | 1.215 | 1.83
ENSG00000122729 | 2.57 | 1.181 | 1.467 | 1.202 | 1.88
ENSG00000167107 | 1.89 | 0.925 | 2.036 | 0.700 | 2.26
ENSG00000130402 | 3.29 | 3.392 | 3.242 | 1.975 | 4.85
ENSG00000173137 | 1.91 | 1.194 | 0.998 | 0.839 | 1.35
ENSG00000035687 | 4.79 | 2.926 | 5.095 | 4.054 | 4.27

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

Example clinical data
ID | Sex | Race | Age_Days | WBC | MRD29 | MRD.binary | os.time | os.censor | efs.time | efs.censor
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
PARASZ | Female | White | 1610 | 27 | 0.000 | 0 | 8.45 | 0 | 8.45 | 0
PARAYM | Male | White | 5569 | 400 | 0.000 | 0 | 9.31 | 0 | 9.31 | 0
PARCVM | Male | White | 1459 | 255 | 0.000 | 0 | 6.64 | 0 | 6.64 | 0
PAREGZ | Female | White | 3128 | 204 | 0.490 | 1 | 8.45 | 0 | 8.45 | 0
PARFDL | Male | White | 7491 | 27 | 0.000 | 0 | 8.42 | 0 | 8.42 | 0
PARFIH | Male | White | 1925 | 102 | 0.028 | 0 | 3.01 | 0 | 3.01 | 0

## Genome annotations

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

    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"
    )

Example gene annotation data
gene | chrom | loc.start | loc.end | description | gene.name | biotype | chrom.strand | chrom.band
--- | --- | --- | --- | --- | --- | --- | --- | ---
ENSG00000122406 | 1 | 9.28e+07 | 9.28e+07 | ribosomal protein L5 [Source:HGNC Symbol;Acc:HGNC:10360] | RPL5 | protein_coding | 1 | p22.1
ENSG00000213281 | 1 | 1.15e+08 | 1.15e+08 | NRAS proto-oncogene, GTPase [Source:HGNC Symbol;Acc:HGNC:7989] | NRAS | protein_coding | -1 | p13.2
ENSG00000147889 | 9 | 2.20e+07 | 2.20e+07 | cyclin dependent kinase inhibitor 2A [Source:HGNC Symbol;Acc:HGNC:1787] | CDKN2A | protein_coding | -1 | p21.3
ENSG00000095574 | 10 | 1.23e+08 | 1.23e+08 | IKAROS family zinc finger 5 [Source:HGNC Symbol;Acc:HGNC:14283] | IKZF5 | protein_coding | -1 | q26.13
ENSG00000184937 | 11 | 3.24e+07 | 3.24e+07 | WT1 transcription factor [Source:HGNC Symbol;Acc:HGNC:12796] | WT1 | protein_coding | -1 | p13
ENSG00000106462 | 7 | 1.49e+08 | 1.49e+08 | enhancer of zeste 2 polycomb repressive complex 2 subunit [Source:HGNC Symbol;Acc:HGNC:3527] | EZH2 | protein_coding | -1 | q36.1

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

chrom | size
--- | ---
1 | 2.49e+08
2 | 2.42e+08
3 | 1.98e+08
4 | 1.90e+08
5 | 1.82e+08
6 | 1.71e+08

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.

    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.

    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
    #> [1] "mutation"
    head(grin.exon.results$gene.exon.size)
    #> [1] 2114 4955 2061 8314 5412 8721
    knitr::kable(head(grin.exon.results$exon.chrom.size), row.names = FALSE)

chrom | size
--- | ---
1 | 10535792
10 | 4368198
11 | 5896028
12 | 6076059
13 | 2476848
14 | 3687507

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.

    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.

    show.grin.columns("^nsubj\\.")

gene | gene.name | chrom | loc.start | loc.end | nsubj.fusion | nsubj.gain | nsubj.loss | nsubj.mutation
--- | --- | --- | --- | --- | --- | --- | --- | ---
ENSG00000171862 | PTEN | 10 | 8.79e+07 | 8.80e+07 | 0 | 2 | 23 | 37
ENSG00000118513 | MYB | 6 | 1.35e+08 | 1.35e+08 | 4 | 26 | 2 | 13
ENSG00000139083 | ETV6 | 12 | 1.16e+07 | 1.19e+07 | 2 | 0 | 15 | 7
ENSG00000111276 | CDKN1B | 12 | 1.27e+07 | 1.27e+07 | 0 | 0 | 20 | 4
ENSG00000184937 | WT1 | 11 | 3.24e+07 | 3.24e+07 | 0 | 0 | 9 | 24
ENSG00000215301 | DDX3X | X | 4.13e+07 | 4.14e+07 | 2 | 1 | 2 | 4

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

    show.grin.columns("^[pq]\\.nsubj\\.")

gene | gene.name | chrom | loc.start | loc.end | p.nsubj.fusion | p.nsubj.gain | p.nsubj.loss | p.nsubj.mutation | q.nsubj.fusion | q.nsubj.gain | q.nsubj.loss | q.nsubj.mutation
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
ENSG00000171862 | PTEN | 10 | 8.79e+07 | 8.80e+07 | 1.00e+00 | 6.32e-01 | 9.15e-37 | 1.48e-78 | 1.00e+00 | 1.00e+00 | 6.96e-35 | 1.54e-76
ENSG00000118513 | MYB | 6 | 1.35e+08 | 1.35e+08 | 3.02e-11 | 7.20e-53 | 1.00e+00 | 2.14e-28 | 3.15e-09 | 3.00e-50 | 9.12e-01 | 7.45e-27
ENSG00000139083 | ETV6 | 12 | 1.16e+07 | 1.19e+07 | 1.15e-04 | 1.00e+00 | 1.33e-13 | 1.07e-07 | 5.97e-03 | 1.00e+00 | 1.04e-12 | 1.49e-06
ENSG00000111276 | CDKN1B | 12 | 1.27e+07 | 1.27e+07 | 1.00e+00 | 1.00e+00 | 2.58e-23 | 3.46e-07 | 1.00e+00 | 1.00e+00 | 1.09e-21 | 4.38e-06
ENSG00000184937 | WT1 | 11 | 3.24e+07 | 3.24e+07 | 1.00e+00 | 1.00e+00 | 4.24e-07 | 1.06e-51 | 1.00e+00 | 1.00e+00 | 1.64e-06 | 6.33e-50
ENSG00000215301 | DDX3X | X | 4.13e+07 | 4.14e+07 | 7.95e-07 | 9.96e-01 | 9.97e-01 | 1.16e-07 | 5.53e-05 | 1.00e+00 | 9.12e-01 | 1.55e-06

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

    show.grin.columns("^[pq][0-9]+\\.nsubj$")

gene | gene.name | chrom | loc.start | loc.end | p1.nsubj | p2.nsubj | p3.nsubj | p4.nsubj | q1.nsubj | q2.nsubj | q3.nsubj | q4.nsubj
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
ENSG00000171862 | PTEN | 10 | 8.79e+07 | 8.80e+07 | 5.93e-78 | 5.03e-72 | 5.31e-01 | 1.00e+00 | 1.37e-76 | 2.10e-69 | 1.00e+00 | 1.00e+00
ENSG00000118513 | MYB | 6 | 1.35e+08 | 1.35e+08 | 2.88e-52 | 2.76e-55 | 1.10e-31 | 1.00e+00 | 4.98e-51 | 5.75e-53 | 4.59e-29 | 1.00e+00
ENSG00000139083 | ETV6 | 12 | 1.16e+07 | 1.19e+07 | 5.32e-13 | 6.86e-14 | 6.02e-12 | 1.00e+00 | 1.53e-12 | 9.53e-12 | 1.26e-09 | 1.00e+00
ENSG00000111276 | CDKN1B | 12 | 1.27e+07 | 1.27e+07 | 1.03e-22 | 7.20e-13 | 1.00e+00 | 1.00e+00 | 8.93e-22 | 7.50e-11 | 1.00e+00 | 1.00e+00
ENSG00000184937 | WT1 | 11 | 3.24e+07 | 3.24e+07 | 4.25e-51 | 1.08e-12 | 1.00e+00 | 1.00e+00 | 6.79e-50 | 8.98e-11 | 1.00e+00 | 1.00e+00
ENSG00000215301 | DDX3X | X | 4.13e+07 | 4.14e+07 | 4.62e-07 | 3.79e-12 | 1.00e+00 | 9.90e-01 | 7.50e-07 | 2.64e-10 | 1.00e+00 | 1.00e+00

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.

    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

    genomewide.lsn.plot(
      standard.grin.results,
      max.log10q = 50
    )

Genome-wide distribution and statistical significance of genomic lesions.

### Stacked lesion-count plot

    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
    )

Numbers of subjects affected by different lesion types in selected genes.

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

    lsn.transcripts.plot(
      grin.res = grin.results,
      chrom = 9,
      plot.start = 1,
      plot.end = 138394717,
      transTrack = FALSE,
      show.ideogram = FALSE,
      point.size.mm = 1
    )

Distribution of all lesion types across chromosome 9 without transcript annotations or a chromosome ideogram.

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.

    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)
    #> [1]   6 132
    oncoprint.mtx[, seq_len(min(6, ncol(oncoprint.mtx))), drop = FALSE]
    #>       PARCVM      PARJPL  PARMIH  PARMKM      PARMRF      PARNXJ
    #> PTEN  "mutation;" "loss;" "loss;" "mutation;" "mutation;" "loss;"
    #> LEF1  ""          ""      ""      ""          ""          ""
    #> ETV6  ""          ""      ""      ""          ""          ""
    #> JAK1  ""          ""      ""      ""          ""          ""
    #> MLLT3 ""          "loss;" ""      "loss;"     ""          ""
    #> RB1   ""          ""      ""      ""          ""          ""

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

    str(onco.props, max.level = 1)
    #> List of 3
    #>  $ alter_func          :List of 5
    #>  $ col                 : Named chr [1:4] "black" "blue" "olivedrab" "red"
    #>   ..- attr(*, "names")= chr [1:4] "fusion" "loss" "mutation" "gain"
    #>  $ heatmap_legend_param:List of 3

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

    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.

    gene.lsn.type.mtx <- prep.lsn.type.matrix(
      gene.lsn.overlap,
      min.ngrp = 5
    )

    dim(gene.lsn.type.mtx)
    #> [1] 302 264
    gene.lsn.type.mtx[
      seq_len(min(6, nrow(gene.lsn.type.mtx))),
      seq_len(min(6, ncol(gene.lsn.type.mtx))),
      drop = FALSE
    ]
    #>                 PARASZ PARAYM PARCVM PAREGZ PARFDL PARFIH
    #> ENSG00000005700 "loss" "none" "loss" "none" "none" "loss"
    #> ENSG00000010810 "loss" "none" "loss" "none" "none" "loss"
    #> ENSG00000014123 "loss" "none" "loss" "none" "none" "loss"
    #> ENSG00000056972 "loss" "none" "loss" "none" "none" "loss"
    #> ENSG00000057663 "loss" "none" "loss" "none" "none" "loss"
    #> ENSG00000065615 "loss" "none" "loss" "none" "none" "loss"

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

    lsn.binary.mtx <- prep.binary.lsn.mtx(
      gene.lsn.overlap,
      min.ngrp = 5
    )

    dim(lsn.binary.mtx)
    #> [1] 287 264
    lsn.binary.mtx[
      seq_len(min(6, nrow(lsn.binary.mtx))),
      seq_len(min(6, ncol(lsn.binary.mtx))),
      drop = FALSE
    ]
    #>                          PATXKW PASHNK PARMUC PATKWU PARXMV PARMIH
    #> ENSG00000005339_mutation      0      1      1      1      1      1
    #> ENSG00000005700_loss          0      0      0      0      0      0
    #> ENSG00000006283_gain          0      0      0      0      0      0
    #> ENSG00000010438_loss          0      0      0      0      0      0
    #> ENSG00000010810_loss          0      0      0      0      0      0
    #> ENSG00000014123_loss          0      0      0      0      0      0

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

    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)
    #> [1] 302 263
    dim(alex.data$alex.expr)
    #> [1] 302 263

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

    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"
    )

Genes with the smallest Kruskal–Wallis q-values
gene | gene.name | p.KW | q.KW
--- | --- | --- | ---
ENSG00000198642 | KLHL9 | 7.14e-21 | 1.96e-18
ENSG00000120159 | CAAP1 | 2.42e-20 | 3.32e-18
ENSG00000162367 | TAL1 | 1.15e-18 | 1.06e-16
ENSG00000137073 | UBAP2 | 4.61e-18 | 3.17e-16
ENSG00000107185 | RGP1 | 1.47e-16 | 8.06e-15
ENSG00000147889 | CDKN2A | 5.65e-16 | 2.59e-14

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

    WT1.waterfall.data <- alex.waterfall.prep(
      alex.data,
      alex.kw.results,
      "WT1",
      lesion_data
    )

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

WT1 expression and lesion groups across subjects.

## Association with clinical outcomes

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

    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.

    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"
      )
    )

Gene-lesion associations ordered by EFS Cox-model p-value, with corresponding MRD logistic-regression results
Gene_lsn | gene.name | cox_EFS_pval | cox_EFS_qval | logistic_MRD.binary_pval | logistic_MRD.binary_qval
--- | --- | --- | --- | --- | ---
ENSG00000102974_mutation | CTCF | 5.65e-04 | 8.75e-02 | 7.22e-01 | 4.93e-01
ENSG00000118513_gain | MYB | 6.01e-03 | 4.45e-01 | 8.11e-02 | 9.90e-02
ENSG00000171862_loss | PTEN | 1.37e-02 | 4.45e-01 | 3.05e-01 | 2.76e-01
ENSG00000088038_mutation | CNOT3 | 1.37e-02 | 4.45e-01 | 7.67e-01 | 5.07e-01
ENSG00000142208_mutation | AKT1 | 2.64e-02 | 4.45e-01 | 2.06e-01 | 1.95e-01
ENSG00000177076_loss | ACER2 | 7.97e-02 | 4.45e-01 | 5.25e-03 | 4.02e-02

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

    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"
    )

Genes with the smallest EFS log-rank p-values
gene | gene.name | logRank_EFS_pval | logRank_EFS_qval
--- | --- | --- | ---
ENSG00000102974 | CTCF | 1.18e-04 | 2.88e-02
ENSG00000088038 | CNOT3 | 1.31e-03 | 1.59e-01
ENSG00000118513 | MYB | 6.92e-03 | 5.54e-01
ENSG00000171862 | PTEN | 1.12e-02 | 5.54e-01
ENSG00000142208 | AKT1 | 3.63e-02 | 5.54e-01
ENSG00000177076 | ACER2 | 4.61e-02 | 5.54e-01

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

    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"
      )
    )

Gene-expression associations ordered by the adjusted EFS Cox-model p-value, with corresponding adjusted MRD logistic-regression results
gene | gene.name | logistic_MRD.binary_pval.adj | logistic_MRD.binary_qval.adj | cox_EFS_pval_adj | cox_EFS_qval_adj
--- | --- | --- | --- | --- | ---
ENSG00000145916 | RMND5B | 2.46e-01 | 2.33e-01 | 0.026 | 0.997
ENSG00000056972 | TRAF3IP2 | 7.15e-01 | 4.82e-01 | 0.031 | 0.997
ENSG00000221946 | FXYD7 | 1.15e-02 | 2.88e-02 | 0.037 | 0.997
ENSG00000070610 | GBA2 | 1.31e-03 | 4.61e-03 | 0.046 | 0.997
ENSG00000137103 | TMEM8B | 6.12e-01 | 4.32e-01 | 0.051 | 0.997
ENSG00000127152 | BCL11B | 9.79e-05 | 5.94e-04 | 0.062 | 0.997

## 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