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

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

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

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

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

## ----standard-grin------------------------------------------------------------
standard.grin.results <- grin.stats(
  lsn.data = lesion_data,
  gene.data = hg38_gene_annotation,
  chr.size = hg38_chrom_size
)

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

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

## ----subject-count-results----------------------------------------------------
show.grin.columns("^nsubj\\.")

## ----lesion-probability-results-----------------------------------------------
show.grin.columns("^[pq]\\.nsubj\\.")

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

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

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

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

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

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

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

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

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

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

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

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