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:

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:

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.

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.

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.

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.

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