NB-FOXR2 project [source]


1 Configuration

Configuration of project directory & analysis outputs:

Show full config

source(here("rr_helpers.R"))

# Set up outputs
message("Document index: ", doc_id)
## Document index: 06.1
# Specify where to save outputs
out        <- here("output", doc_id); dir.create(out, recursive = TRUE)
figout     <- here("figures", doc_id, "/")
cache      <- paste0("/project/kleinman/bhavyaa.chandarana/cache/NB-FOXR2/public/", doc_id, "/")

The root directory of this project is:

## /project/kleinman/bhavyaa.chandarana/from_hydra/2023-05-NB-FOXR2/public

Outputs and figures will be saved at these paths, relative to project root:

## public/output/06.1
## public/figures/06.1//

Setting a random seed:

set.seed(100)

2 Overview

This document explores projections/mapping of tumors to normal using both bulk and single-cell RNAseq data.

With bulk RNAseq data, the strategy is to evaluate enrichment of cell type-specific gene signatures (derived in 02) in each bulk transcriptome. We can then evaluate the highest scoring signature(s) in each sample, or compare enrichment scores of gene signatures between tumor types.

3 Libraries & functions

library(here)
library(magrittr)
library(tidyr)
library(dplyr)
library(readr)
library(readxl)
library(stringr)
library(glue)
library(purrr)
library(ggplot2)
library(ggrepel)
library(tibble)
library(broom)
library(cowplot)
library(Seurat)
library(ComplexHeatmap)
library(data.table)
library(ggExtra)
library(Rtsne)

source(here("include/style.R"))
source(here("code/functions/RNAseq.R"))
source(here("code/functions/scRNAseq.R"))
source(here("code/functions/ssGSEA.R"))

ggplot2::theme_set(theme_min())

square_theme <- theme(aspect.ratio = 1)
large_text <- theme(text = element_text(size = 15))

conflicted::conflicts_prefer(dplyr::filter)
# Function to load an RData file and return it
loadRData <- function(fileName){
  
  load(fileName)
  get(ls()[ls() != "fileName"])
  
}

4 Load metadata

meta <- read_tsv(here("output/00/metadata_patients_NGS.tsv"))
## Rows: 176 Columns: 20
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (20): Sample, ID_Patient, ID_Sample, ID, FOXR2_positive, Group, Source, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
meta_bulk <- meta %>% filter(RNAseq == "Y")

5 Load signatures

Load cell type specific gene signatures, compiled in document 02, to use as input for ssGSEA:

load(here("output/02/signatures.Rda"))

6 Load counts

6.1 Pediatric brain tumors (PBT)

Load the raw counts from the bulk RNAseq pipeline, for the test set of pediatric brain tumors:

# load counts aggregated with in-house pipeline
test_raw_counts <- file.path(here("data/RNAseq/pipeline_l3/2023-05-test_pbt/counts/Ensembl.ensGene.exon.raw.tsv.gz")) %>%
    read.table(header = T, row.names = 1, check.names = F, sep = "\t")

# rename genes with ENSEMBL ID only
head(rownames(test_raw_counts) %<>% strsplit(":") %>% sapply(getElement, 1))
## [1] "ENSG00000118473" "ENSG00000162426" "ENSG00000157191" "ENSG00000169504"
## [5] "ENSG00000142920" "ENSG00000186094"
test_raw_counts <- as.matrix(test_raw_counts)

test_raw_counts[1:5, 1:5]
##                 00-009P 00-011R 00-020P 00-027P 00-029P
## ENSG00000118473     118    1400    1873    1978    1274
## ENSG00000162426      16     345     482     839     558
## ENSG00000157191     947    2035    3668    6433    6711
## ENSG00000169504    1947   14701   31704   19765   29772
## ENSG00000142920     142     394     930    1465    1241

Load raw counts for the background set of brain tumors/normal brain:

# load counts aggregated with in-house pipeline
background_raw_counts <- file.path(here("data/RNAseq/pipeline_l3/2023-04-background_samples/counts/Ensembl.ensGene.exon.raw.tsv.gz")) %>%
    read.table(header = T, row.names = 1, check.names = F, sep = "\t")

# rename genes with ENSEMBL ID only
head(rownames(background_raw_counts) %<>% strsplit(":") %>% sapply(getElement, 1))
## [1] "ENSG00000118473" "ENSG00000162426" "ENSG00000157191" "ENSG00000169504"
## [5] "ENSG00000142920" "ENSG00000186094"
background_raw_counts <- as.matrix(background_raw_counts)

background_raw_counts[1:5, 1:5]
##                 EBT_1 EBT_2 EBT_3 EBT_4 EBT_5
## ENSG00000118473   384   557    54   160   515
## ENSG00000162426   154   128   199   176   118
## ENSG00000157191  2205  1570  1133   992  1065
## ENSG00000169504  7780  3180  9300  5670  4520
## ENSG00000142920   303   246   512   230   182

6.2 Extracranial neuroblastoma (EC-NB)

For extra-cranial neuroblastomas from Gartlguber et al, Nature Cancer, 2021 (n=579) processed data was obtained from the Shiny App associated with the study.

# load raw counts provided by the authors
ecnb_counts <- readRDS(here("data/RNAseq/external_data/Gartlgruber_NatCancer_2021/tumor_RNAseq_Counts_Matrix.rds"))

head(rownames(ecnb_counts))
## [1] "ENSG00000223972.4|DDX11L1"    "ENSG00000227232.4|WASH7P"    
## [3] "ENSG00000243485.2|MIR1302-11" "ENSG00000237613.2|FAM138A"   
## [5] "ENSG00000268020.2|OR4G4P"     "ENSG00000240361.1|OR4G11P"
# extract ENSEMBL id only from gene names
head(rownames(ecnb_counts) %<>% strsplit("\\.") %>% sapply(getElement, 1))
## [1] "ENSG00000223972" "ENSG00000227232" "ENSG00000243485" "ENSG00000237613"
## [5] "ENSG00000268020" "ENSG00000240361"
ecnb_counts <- as.matrix(ecnb_counts)

ecnb_counts[1:5, 1:5]
##                 NSP001-PT01 NSP002-PT01 NSP003-RM01 NSP005-RT01 NSP006-RM01
## ENSG00000223972           4           0           0           3           1
## ENSG00000227232         431        1348        1242        1117        1391
## ENSG00000243485           0           0           0           0           0
## ENSG00000237613           0           0           0           0           0
## ENSG00000268020           0           0           0           0           0

7 Combine metadata

7.1 Load EC-NB metadata

Load the sample metadata for Gartlgruber EC-NB dataset:

meta_ecnb <- readRDS(here("data/RNAseq/external_data/Gartlgruber_NatCancer_2021/annotation_tumor_RNAseq.rds"))

# load annotation of ChIPseq data to look for overlapping samples
meta_ecnb_chip <- readRDS(here("data/ChIPseq/external_data/Gartlgruber_NatCancer_2021/annotation_tumor_ChIPseq.rds"))

meta_ecnb <- meta_ecnb %>%
    left_join(meta_ecnb_chip, by = c("ProjectID", "MYCN", "Stage", "Age", "Risk", "Relapse")) %>% 
    dplyr::rename(ID = ProjectID) %>%
    mutate(Group = case_when(
        MYCN == "Amp" ~ "EC-NB (MYCN-amp)",
        TRUE ~ "EC-NB"))

Add metadata column for FOXR2 expression status in EC-NB:

# Load files and remove unused ones to save space
load(here("output/03/Gartlgruber_et_al_counts.Rda"))
rm(ecnb_counts_norm)
rm(ecnb_counts_vst)

# DESeq norm expression threshold used for FOXR2 +/-
threshold <- 2

foxr2_pos_samples <- ecnb_counts_tidy %>% 
    filter(gene_symbol == "FOXR2") %>% 
    filter(gene_expression > threshold) %>% 
    .$sample

meta_ecnb <- meta_ecnb %>% mutate(FOXR2_positive = case_when(ID %in% foxr2_pos_samples ~ "Y",
                                      T ~ "N")) %>% 
    mutate(Group = "EC-NB") %>% 
    mutate(Source = "Gartlgruber et al. 2021")

rm(ecnb_counts_tidy)

7.2 Combine with PBT metadata

Combine with brain tumor dataset:

meta_all <- bind_rows(meta_bulk, meta_ecnb)

8 Run ssGSEA

Run ssGSEA for each dataset:

# test set
bulk_ssgsea <- compute_scores_bin(expr_mat  = test_raw_counts,
                                  gene_sets = signatures_ens_all,
                                  save_le   = FALSE,
                                  alpha     = 0.75,
                                  normalize = FALSE,
                                  n_cores   = 1)

saveRDS(bulk_ssgsea, file = glue("{out}/ssgsea_test.Rds"))
# background set
bulk_ssgsea_bg <- compute_scores_bin(expr_mat  = background_raw_counts,
                                     gene_sets = signatures_ens_all,
                                     save_le   = FALSE,
                                     alpha     = 0.75,
                                     normalize = FALSE,
                                     n_cores   = 1)

saveRDS(bulk_ssgsea_bg, file = glue("{out}/ssgsea_background.Rds"))
# extra-cranial NB
bulk_ssgsea_ecnb <- compute_scores_bin(expr_mat  = ecnb_counts,
                                       gene_sets = signatures_ens_all,
                                       save_le   = FALSE,
                                       alpha     = 0.75,
                                       normalize = FALSE,
                                       n_cores   = 1)

saveRDS(bulk_ssgsea_ecnb, file = glue("{out}/ssgsea_ecnb.Rds"))

9 Identify indiscriminant signatures

First, we will evaluate the ssGSEA scores for reference signatures in our background set of tumors and normal brain. To do this, we count the number of samples where each signature appears in the top 10. We'll consider the top decile as indiscriminate.

# control set
bulk_ssgsea_bg <- readRDS(glue("{out}/ssgsea_background.Rds"))
dim(bulk_ssgsea_bg)
## [1] 412 211
# filter to retained signatures
bulk_ssgsea_bg_keep <- bulk_ssgsea_bg %>%
    filter(Signature %in% signatures_keep)

dim(bulk_ssgsea_bg_keep)
## [1] 390 211
# tidy dataframe
bulk_ssgsea_bg_top <- bulk_ssgsea_bg_keep %>%
    gather(Sample, Score, 2:ncol(.)) %>% 
    group_by(Sample) %>% 
    top_n(10, Score)

# get n signatures - this many signatures appear in the top 10 of any sample
n_sigs_top <- length(unique(bulk_ssgsea_bg_top$Signature))
n_sigs_top
## [1] 135
# count the number of times each signature appears in the top 10 signatures for any sample
counts_top <- bulk_ssgsea_bg_top %>% 
    group_by(Signature) %>% 
    dplyr::count() 

# show the values
counts_top %>% arrange(desc(n)) %>% DT::datatable()
# get threshold at top decile
(q <- quantile(counts_top$n, 0.9)[[1]])
## [1] 41

Then, we can identify indiscriminant signatures by calculating a threshold \(q\) where the probability that \(X\) is less than or equal to \(q\) is 90%.

counts_top %>% 
    ggplot(aes(x = n)) +
    geom_histogram() +
    geom_vline(xintercept = q, color = "red") +
    xlab("Number of samples where the signatures is in the top 10") +
    ylab("Number of signatures")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

In other words, we identify indiscriminant signatures as those which are in the top 10 of more samples than 90% of other signatures. As a sanity check, we have 135 signatures appearing in the top 10 of any sample, and then are taking the top 13.5 signatures as indiscriminant.

# indiscriminant signatures
(sigs_indiscriminant <- counts_top %>% filter(n >= q) %>% pull(Signature))
##  [1] "CGE-ANK3+,ENC1+"                  "CGE-GRIA1+,SST+"                 
##  [3] "F-e15_CEXN2"                      "HF OPC"                          
##  [5] "HP AST-FB"                        "HP Neu-NRGN-I"                   
##  [7] "HP Neu-NRGN-II"                   "HP OPC"                          
##  [9] "HP Oligodendrocytes"              "Human fetal GE OPC"              
## [11] "Human fetal GE P3"                "Human fetal cerebellum 09-Ast"   
## [13] "Str/SVZ 33-Neurogenic progenitor" "VSVZ-1 Neuron"                   
## [15] "VSVZ-2 Neuron"
length(sigs_indiscriminant)
## [1] 15

In addition, we also filter out a MGE signature which we determined to be not specific to its associated cell type, based on expression of the genes in this signature across multiple cell types. (See following section.)

sigs_indiscriminant <- c("MGE-NR2F1+,MEIS2+", sigs_indiscriminant)
length(sigs_indiscriminant)
## [1] 16
save(sigs_indiscriminant, file = glue("{out}/sigs_indiscriminant.Rda"))

9.1 Non-specific MGE signature

In previous iterations of this analysis, MGE signature MGE-NR2F1+,MEIS2+ appeared in top ssGSEA matches in the majority of DIPG-H3K27M. This was unexpected because DIPG-H3K27M appear in the pons, not in the forebrain, where the MGE structure appears during development.

Let us check the detection rate of each gene of this signature in cells from a normal developmental atlas, namely MGE neurons vs. cell types that give rise to DIPG (OPC and oligodendrocytes) in both the pons and forebrain. We will also compare this with non-neural cell types (vascular, immune, endothelial) as a negative control.

Define a function to plot detection rate per gene in a given signature:

# Load full cell type ontology table for mouse developing forebrain and pons
ontology_jessa_full <- fread("/project/kleinman/selin.jessa/from_narval/HGG-oncohistones/public/output/05/TABLE_mouse_cluster_info.tsv") %>%
  select(Label, Level3_type) %>%
  rename(Label_broad = Level3_type)

# Function producing box plot of detection rate per gene in a signature
# Within cell populations: MGE neurons, Forebrain OPC & oligodendrocytes, Pons OPC & oligodendrocytes, and non-neural cells
plot_mge_neur_opc_box_plot <- function(signature,
                                    signature_id,
                                    conv_hg_to_mm){
  print(signature_id)

  ontology <- ontology_jessa_full

  if(conv_hg_to_mm){

    # Convert genes from human to mouse using cached biomaRt reference

    # Produces file "genes_lds"
    load(here("data/singlecell/references_genome/biomaRt_mm_to_hg_lds.Rda"))

    # Store converted signature
    mm_signature <- genes_lds %>% filter(HGNC.symbol %in% signature)

    # Print any human genes that could not be converted
    if(length(mm_signature) != length(signature)){

      not_conv <- setdiff(signature, mm_signature$HGNC.symbol)
      print(glue("The following genes of the provided signature could not be converted from human to mouse using biomaRt:"))
      print(list(not_conv))

    }

    # Replace user-provided human signature with signature in mouse
    signature <- mm_signature$MGI.symbol
    rm(genes_lds)

  }

  # Load developing mouse atlas objects
  atlas_dir <- "/project/kleinman/selin.jessa/from_hydra/atlas/data"
  seurat_ct <- loadRData(glue("{atlas_dir}/joint_cortex_extended/joint_cortex_extended.seurat_v3.Rda"))
  suppressMessages({
    seurat_ct <- UpdateSeuratObject(seurat_ct)
  })
  seurat_po <- loadRData(glue("{atlas_dir}/joint_pons_extended/joint_pons_extended.seurat_v3.Rda"))
  suppressMessages({
    seurat_po <- UpdateSeuratObject(seurat_po)
  })

  # Store forebrain and pons objects in list and remove the originals
  seurat_list <- c(seurat_ct, seurat_po)
  names(seurat_list) <- c("cortex", "pons")
  rm(seurat_ct)
  rm(seurat_po)

  # Add broad labels and categories for plot to Seurat objects
  for(seurat_id in names(seurat_list)){

    seurat <- seurat_list[[seurat_id]]
    # Get metadata from Seurat object
    labels <- seurat@meta.data %>%
      select(ID_20210710) %>%
      rownames_to_column(var = "cellname")
    rownames(labels) <- NULL

    # Add broad labels
    labels <- labels %>%
      left_join(ontology %>% select(Label, Label_broad),
                by = c('ID_20210710' = 'Label')) %>%
      column_to_rownames(var = "cellname")
    seurat <- AddMetaData(seurat, labels)

    # Group cells of selected cell type(s) for different rows of plot
    if(seurat_id == "cortex"){

      labels_ct <- labels %>%
    mutate(plot_row = case_when(Label_broad %in% c("MGE inhibitory neurons") ~ "MGE inhibitory",
                                Label_broad %in% c("OPC", "Oligodendrocytes") ~ "OPC/Oligo (cortex)",
                                Label_broad %in% c("Meninges", "Endothelial", "Immune", "Pericytes") ~ "Endothelial/Immune (cortex)",
                                TRUE ~ "Other"))
      seurat <- AddMetaData(seurat, labels_ct)

    } else if(seurat_id == "pons"){

      labels_po <- labels %>%
    mutate(plot_row = case_when(Label_broad %in% c("OPC", "Oligodendrocytes") ~ "OPC/Oligo (pons)",
                                TRUE ~ "Other"))
      seurat <- AddMetaData(seurat, labels_po)

    }

    seurat_list[[seurat_id]] <- seurat
  }

  # Retrieve detection rate for each gene in the signature

  seurat_list[["cortex"]] <- seurat_list[["cortex"]][signature, ]
  seurat_list[["pons"]] <- seurat_list[["pons"]][signature, ]
  mat_list <- c()
  mat_list[["cortex"]] <- as.matrix(seurat_list[["cortex"]]@assays$RNA@data)
  mat_list[["pons"]] <- as.matrix(seurat_list[["pons"]]@assays$RNA@data)
  pct_list <- c()

  for(name in names(mat_list)){

    x <- mat_list[[name]]

    # Binarize
    x[x > 0] <- 1

    # Transpose preserving cell names
    cellname <- colnames(x)
    x <- as.data.table(t(x))
    rownames(x) <- cellname

    # Add in_group info for cells
    x[, plot_row := as.character(seurat_list[[name]]@meta.data[["plot_row"]])]

    # Get prop of cells expressing a gene, within group and out of group
    pct <- x[, lapply(.SD, prop), by = plot_row]

    # Convert to tidy format for plotting
    pct <- pct %>%
      pivot_longer(cols = colnames(pct)[-1], names_to = "Gene")

    pct_list[[name]] <- pct
  }

  # Clean up
  rm(seurat_list)
  rm(mat_list)

  # Join cortex and pons detection rates, removing "Other" cell type category from each
  pct_list[["cortex"]] <- pct_list[["cortex"]] %>%
    as.data.frame() %>%
    filter(plot_row != "Other") %>%
    rename(Cell_class = plot_row)
  pct_list[["pons"]] <- pct_list[["pons"]] %>%
    as.data.frame() %>%
    filter(plot_row != "Other") %>%
    rename(Cell_class = plot_row)
  pct <- rbind(pct_list[["cortex"]], pct_list[["pons"]])

  suppressMessages({
  # Produce horizontal box plot
  pct %>%
    mutate(Cell_class = factor(Cell_class,
                    levels = rev(c("MGE inhibitory",
                               "OPC/Oligo (cortex)",
                               "OPC/Oligo (pons)",
                               "Endothelial/Immune (cortex)")))) %>%
    ggplot(aes(x = value, y = Cell_class)) +
    ggtitle(glue("{signature_id}")) +
    geom_boxplot() +
    xlab("% of cells expressing gene") +
    ylab("") +
    scale_x_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     limits = c(-0.05, 1.05)) +
    theme(axis.title.x = element_text(hjust = 0.5)) 
  })

}

Plot detection rate in the signature of interest MGE-NR2F1+,MEIS2+, along with two other MGE signatures used in the study:

sigs <- c("MGE-NR2F1+,MEIS2+",
          "F-e12_MGINH",
          "Human fetal GE MGE1")

plot_list <- map(sigs,
    ~ plot_mge_neur_opc_box_plot(signature = signatures_sym_all[[.x]],
                              signature_id = .x,
                              conv_hg_to_mm = T) + large_text)
## [1] "MGE-NR2F1+,MEIS2+"
## The following genes of the provided signature could not be converted from human to mouse using biomaRt:
## [[1]]
## [1] "NNAT"    "BEX3"    "HMGN2"   "TMSB15A" "PPIA"    "VDAC1"   "TCEAL2" 
## [8] "LDHA"    "RACK1"  
## 
## [1] "F-e12_MGINH"
## The following genes of the provided signature could not be converted from human to mouse using biomaRt:
## [[1]]
##  [1] "NNAT"      "H3F3B"     "HIST3H2BB" "ATPIF1"    "AES"       "USMG5"    
##  [7] "AHSA2"     "UBL5"      "DNAJA1"    "MAP1LC3B2"
## 
## [1] "Human fetal GE MGE1"
## The following genes of the provided signature could not be converted from human to mouse using biomaRt:
## [[1]]
##  [1] "DLX6-AS1"     "CADPS"        "CEMIP2"       "RIPOR2"       "TMEM161B-AS1"
##  [6] "C1orf61"      "MIR181A1HG"   "AES"          "HIST1H2BG"    "MALAT1"      
## [11] "MEG3"         "AKR1C1"       "KCNQ1OT1"
plot_grid(plotlist = plot_list, nrow = 1)

The MGE-NR2F1+,MEIS2+ gene signature actually has a lower detection rate in MGE inhibitory neurons than the OPC of the forebrain and pons. Therefore, we remove this signature from our analysis for non-specific expression to the MGE.

9.2 TABLE: Reference signatures

Here, we'll prepare a supplementary table containing all the reference gene signatures and their associated cell types/datasets from which they originate.

# collapse signatures into strings
sigs_sym_df <- imap_dfr(signatures_sym_all,
                        ~ data.frame(Cluster = .y, Signature_sym = glue_collapse(.x, sep = ",")))

sigs_ens_df <- imap_dfr(signatures_ens_all,
                        ~ data.frame(Cluster = .y, Signature_ENS = glue_collapse(.x, sep = ",")))

# join with cell type annotation
TABLE_ref_sigs <- cell_type_anno_all %>% 
    left_join(sigs_sym_df) %>% 
    left_join(sigs_ens_df) %>%
    # filter to kept signatures
    filter(Keep) %>% 
    select(-Keep) %>%
    mutate(Cell_type = ifelse(is.na(Cell_type), Cluster, Cell_type)) %>% 
    mutate(Included = ifelse(Cluster %in% sigs_indiscriminant, "N", "Y"))
## Joining with `by = join_by(Cluster)`
## Joining with `by = join_by(Cluster)`
rr_write_tsv(TABLE_ref_sigs,
             glue("{out}/TABLE_reference_signatures.tsv"),
             "Cell type specific reference signatures and associated metadata and dataset info")
## ...writing description of TABLE_reference_signatures.tsv to public/output/06.1/TABLE_reference_signatures.desc

10 Tidy ssGSEA scores

Combine the ssGSEA scores from the two datasets:

bulk_ssgsea_test <- readRDS(glue("{out}/ssgsea_test.Rds"))
bulk_ssgsea_ecnb <- readRDS(glue("{out}/ssgsea_ecnb.Rds"))

scores_all <- full_join(bulk_ssgsea_test, bulk_ssgsea_ecnb, by = "Signature") %>%
    as.data.frame() %>%
    # filter out short signatures
    filter(Signature %in% signatures_keep) %>% 
    # filter out the indiscriminant signature
    filter(!(Signature %in% sigs_indiscriminant))

scores_all$Signature %>% unique %>% length()
## [1] 374
# correct a sample naming issue, where samples with numeric names XX get converted to XX.0
# and save a map between the coerced and cleaned names
sample_naming_map <- data.frame(Sample = colnames(scores_all))
sample_idx <- which(grepl("\\.0$", colnames(scores_all)))
colnames(scores_all)[sample_idx] <- gsub("\\.0$", "", colnames(scores_all)[sample_idx])
sample_naming_map$Sample_clean <- colnames(scores_all)

# get the set of test samples
test_samples <- base::setdiff(colnames(bulk_ssgsea_test), "Signature")
test_samples <- plyr::mapvalues(test_samples, from = sample_naming_map$Sample, to = sample_naming_map$Sample_clean, warn_missing = FALSE)

scores_tidy <- scores_all %>%
    # convert to tidy/long format
    gather("Sample", "Score", 2:ncol(.)) %>%
    # join with sample metadata
    left_join(meta_all, by = c("Sample" = "ID")) %>%
    # join with signature metadata
    left_join(cell_type_anno_all %>% dplyr::rename("Reference_sample" = Sample, "Reference_age" = Age),
              by = c("Signature" = "Cluster")) %>%
    # set ordering
    mutate(Group = factor(Group, levels = names(palette_groups))) %>%
    arrange(Group) %>%
    mutate(Sample = factor(Sample, levels = unique(.$Sample))) %>% 
    mutate(Class = factor(Class, levels = names(palette_class)))

save(scores_all, scores_tidy, test_samples, file = glue("{out}/ssgsea_scores_tidy.Rda"))

# count samples
scores_tidy %>% distinct(Sample, Group) %>% group_by(Group) %>% count()

11 Dimensionality reduction

11.1 Intracranial tumors alone

In-house processed data for PBT alone:

perp <- 10
                             
tsne_plot <- scores_all[, c("Signature", test_samples)] %>%
     tibble::column_to_rownames(var = "Signature") %>% 
     as.matrix() %>% 
     ssgsea_dim_red(method = "tsne", return_df = TRUE, perplexity = perp) %>%
     left_join(meta_all, by = c("Sample" = "ID")) %>%
     plot_dim_red(method = "tsne", palette = palette_groups,
                  point_size = 3, alpha = 0.8, title = glue("Perplexity: {perp}")) + square_theme + large_text
## @ returning data...
tsne_plot

11.2 Intra + extracranial tumors

PBT with Gartlgruber dataset High risk & Stage 4 EC-NB. (EC-NB samples with NA in either Risk or Stage metadata columns are removed.)

Note that this dimensionality reduction tSNE plot does not exactly match the published figure, as the algorithm is stochastic, and a random seed was not set at the time of creating publication figure.

palette_groups_mut <- palette_groups
palette_groups_mut <- c(palette_groups_mut, "EC-NB: MYCN Amp" = "#A575D3")
palette_groups_mut <- c(palette_groups_mut, "EC-NB: FOXR2+, MYCN NonAmp" = "#f382c6")
palette_groups_mut <- c(palette_groups_mut, "EC-NB: FOXR2-, MYCN NonAmp" = "#215ba3")

perplexity <- 30

meta_ecnb <- meta_ecnb %>%
    mutate(Group = case_when(
        MYCN == "Amp" ~ "EC-NB: MYCN Amp",
        MYCN == "NonAmp" & FOXR2_positive == "Y" ~ "EC-NB: FOXR2+, MYCN NonAmp",
        MYCN == "NonAmp" & FOXR2_positive == "N" ~ "EC-NB: FOXR2-, MYCN NonAmp"))

meta_all <- bind_rows(meta_bulk, meta_ecnb)

non_stage4_HR_ecnb <- meta_ecnb %>% 
    filter(Stage != "4" | Risk != "HR" | is.na(Stage) | is.na(Risk)) %>% 
    .$ID

length(non_stage4_HR_ecnb)
## [1] 387
length(filter(meta_ecnb, !(ID %in% non_stage4_HR_ecnb))$ID)
## [1] 192
set.seed(100)
tsne_calc <- scores_all %>%
                select(c("Signature", test_samples, filter(meta_ecnb, !(ID %in% non_stage4_HR_ecnb))$ID)) %>% 
                 tibble::column_to_rownames(var = "Signature") %>% 
                 as.matrix() %>% 
                 ssgsea_dim_red(method = "tsne", return_df = TRUE, perplexity = perplexity)
## @ returning data...
(tsne_plot <- tsne_calc %>% left_join(meta_all, by = c("Sample" = "ID")) %>% 
                            .[which(!(.$Sample %in% non_stage4_HR_ecnb)),] %>% 
                            plot_dim_red(method = "tsne", palette = palette_groups_mut,
                                          point_size = 2, alpha = 0.8, title = glue("High risk & Stage 4 EC-NB only, Perplexity: {perplexity}")) + square_theme + large_text)

12 Absolute scores

Check if distributions are similar:

scores_tidy %>% 
    ggplot(aes(x = Sample, y = Score)) +
    geom_boxplot(aes(fill = Group), outlier.shape = NA) +
    scale_fill_manual(values = palette_groups) +
    rotate_x() +
    ggtitle("Distribution of ssGSEA scores (across signatures) in each sample") +
    theme(axis.text.x = element_blank())

12.1 Top scoring signature per tumor

Pediatric brain tumor cohort:

scores_tidy %>%
    filter(Sample %in% test_samples) %>% 
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(10, Score) %>%
    ggplot(aes(x = Sample, y = Score)) +
    geom_point(aes(colour = Class), size = 2, alpha = 0.8) +
    scale_colour_manual(values = palette_class) +
    facet_grid(~ Group, space = "free", scales = "free_x") +
    rotate_x() +
    ggtitle("Top 10 signatures in each sample, PBTs, coloured by cell class")

Gartlgruber 2021 extracranial neuroblastoma:

scores_tidy %>%
    filter(Sample %in% meta_ecnb$ID) %>% 
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(10, Score) %>%
    ggplot(aes(x = Sample, y = Score)) +
    geom_point(aes(colour = Class), size = 0.5, alpha = 0.8) +
    scale_colour_manual(values = palette_class) +
    facet_grid(~ Group, space = "free", scales = "free_x") +
    rotate_x() +
    theme(axis.text.x = element_blank()) +
    ggtitle("Top 10 signatures in each sample, EC-NB Gartlgruber 2021, coloured by cell class")

Top score in each sample, across pediatric brain tumor groups:

scores_top_pbt <- scores_tidy %>% 
    filter(Sample %in% test_samples) %>%
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(1, Score) %>% 
    arrange(desc(Score)) %>%
    mutate(Sample = factor(Sample, levels = unique(.$Sample)))

scores_top_pbt %>% 
    ggplot(aes(x = Sample, y = Score)) +
    geom_col(aes(fill = Class), width = 0.8) +
    scale_fill_manual(values = palette_class) +
    xlab("Sample") + ylab(NULL) +
    geom_text(aes(label = Signature), angle = 90, nudge_y = 200, hjust = 0) +
    facet_grid(~ Group, scales = "free_x", space = "free") +
    no_legend() +
    theme_min2() +
    theme(axis.text.x = element_blank()) +
    coord_cartesian(ylim = c(22500, 35000))

# Print top signature for each sample
scores_top_pbt %>% group_by(Group) %>%
    mutate(n_group = n()) %>% 
    group_by(Group, Signature) %>% 
    summarize(n = n(), prop = n / n_group) %>%
    distinct() %>% 
    arrange(Group, desc(n))
## `summarise()` has grouped output by 'Group', 'Signature'. You can override
## using the `.groups` argument.

Top matching signature per tumor type:

scores_tidy %>%
    filter(Sample %in% test_samples) %>%
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(1, Score) %>%
    mutate(Group = factor(Group,
                          levels = c("NB-FOXR2",
                                     "ETMR",
                                     "MB-SHH",
                                     "MB-WNT",
                                     "DIPG-H3K27M-FOXR2",
                                     "DIPG-H3K27M",
                                     "HGG-IDH",
                                     "EP-PFA",
                                     "HGG-H3.3G34R/V"))) %>%
    ggplot(aes(x = Group)) +
    geom_bar(aes(fill = Class), size = 2, alpha = 0.8) +
    scale_fill_manual(values = palette_class) +
    theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust=1))

Top signature in NB-FOXR2:

scores_top_nbfoxr2 <- scores_top_pbt %>% filter(Group == "NB-FOXR2")

scores_top_nbfoxr2 %>% 
    ungroup() %>% 
    count(Signature, Class) %>% 
    ggplot(aes(x = reorder(Signature, n),
               y = n,
               fill = Class)) +
    geom_col() +
    scale_fill_manual(values = palette_class) +
    coord_flip() +
    xlab("Signature") +
    ylab("Number of samples") +
    large_text +
    theme(legend.position = "bottom")

Top signature in extracranial neuroblastoma (Gartlgruber cohort):

scores_top_ec <- scores_tidy %>% 
    filter(Sample %in% meta_ecnb$ID) %>%
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(1, Score) %>% 
    arrange(desc(Score)) %>%
    mutate(Sample = factor(Sample, levels = unique(.$Sample)))

scores_top_ec %>% 
    ggplot(aes(x = Sample, y = Score)) +
    geom_col(aes(fill = Class, color = Class)) +
    scale_color_manual(values = palette_class) +
    scale_fill_manual(values = palette_class) +
    xlab("Sample") + ylab(NULL) +
    theme(axis.text.x = element_blank()) +
    coord_cartesian(ylim = c(22500, 30000)) 

# Print percentage of EC-NB samples mapping to each signature
(prop <- scores_top_ec %>%
    group_by(Signature) %>%
    count() %>% 
    mutate(prop = n/sum(.$n)) %>%
    arrange(desc(n)))
# Print total number of EC-NB samples in Gartlgruber cohort
sum(prop$n)
## [1] 579

12.1.1 TABLE: Cell type signature enrichment

Save top scoring signatures for each tumor in a supplementary table:

TABLE_top_scores <- bind_rows(scores_top_pbt, scores_top_ec) %>%
    select(Sample, Group, Signature, Score, Reference_sample,
           Reference_age, Cell_type, Dataset, Species, Class) %>%
    arrange(Group)

dim(TABLE_top_scores)
## [1] 700  10
table(as.character(TABLE_top_scores$Group))
## 
##       DIPG-H3K27M DIPG-H3K27M-FOXR2             EC-NB            EP-PFA 
##                19                 6               579                13 
##              ETMR    HGG-H3.3G34R/V           HGG-IDH            MB-SHH 
##                12                18                10                 8 
##            MB-WNT          NB-FOXR2 
##                10                25
rr_write_tsv(TABLE_top_scores,
             glue("{out}/TABLE_top_scores.tsv"),
             "Top scoring ssGSEA signature for each bulk tumor")
## ...writing description of TABLE_top_scores.tsv to public/output/06.1/TABLE_top_scores.desc

12.2 Top 10 signatures per sample

Define a helper function for per-sample plots of the top 10 signatures:

# Extract top 10 per sample
top_scores <- scores_tidy %>%
    group_by(Sample) %>%
    arrange(desc(Score)) %>%
    top_n(10, Score)

#' @param goi Character, group of interest
plot_scores_per_sample <- function(goi) {
    
    samples <- top_scores %>% filter(Group == goi) %>% pull(Sample) %>% unique()
    
    p1 <- map(samples,
              ~ top_scores %>% 
                  filter(Sample == .x) %>% 
                  arrange(Score) %>%
                  mutate(Signature = factor(Signature, levels = .$Signature)) %>% 
                  ggplot(aes(x = Signature, y = Score)) +
                  geom_col(aes(fill = Class), width = 0.5) +
                  scale_fill_manual(values = palette_class) +
                  theme_min2() +
                  rotate_x() +
                  xlab(NULL) + ylab(NULL) +
                  ggtitle(.x) +
                  no_legend() +
                  coord_flip(ylim = c(20000, 30000)))
    
    plot_grid(plotlist = p1, align = "hv", axis = "rltb", ncol = 4)
    
}

Plot top 10 signatures per NB-FOXR2 sample:

plot_scores_per_sample("NB-FOXR2")

13 Comparative analysis

In this section, we perform a comparative analysis of ssGSEA enrichment scores between tumor types. We'll broadly divide brain tumors into neuronal and glial tumor types, and compare NB-FOXR2 to each of these broad groups. Here, we exclude H3.3G34R/V tumors, since this tumor has an interneuron signal (Chen et al, Cell 2020).

# Label tumors as neuronal or glial, remove G34
scores_tidy_comp <- scores_tidy %>%
    filter(!grepl("EC", Group)) %>% 
    mutate(Group_comp = case_when(
        Group == "NB-FOXR2" ~ "NB-FOXR2",
        grepl("MB|ETMR|HGNET", Group) ~ "Neuronal",
        grepl("HGG|PFA|DIPG", Group) &
            Group != "HGG-H3.3G34R/V" ~ "Glial",
    ))

# sanity check
any(scores_tidy$Signature %in% sigs_indiscriminant)
## [1] FALSE
# compute the mean per signature per group
scores_group_means <- scores_tidy_comp %>% 
    mutate(Score = scales::rescale(Score, to = c(0, 1))) %>% 
    group_by(Signature, Group_comp) %>%
    summarise(group_mean = mean(Score))
## `summarise()` has grouped output by 'Signature'. You can override using the
## `.groups` argument.
# using broom to run pairwise t tests between groups based on the score
# NOTE: see below on differences in defaults between pairwise.t.test and t.test
# https://stackoverflow.com/a/11457871
scores_tidy_ttest <- scores_tidy_comp %>% 
    group_by(Signature) %>% 
    summarise(ttest = list(pairwise.t.test(Score, Group_comp, p.adjust = "none",
                                           paired = FALSE,
                                           pool.sd = FALSE))) %>%
    mutate(ttest = map(ttest, broom::tidy)) %>%
    unnest(cols = c(ttest)) %>% 
    mutate(p_adj = p.adjust(p.value, method = "BH"))

# add the group estimates after pairwise tests for later analysis
scores_tidy_ttest <- scores_tidy_ttest %>% 
    left_join(scores_group_means, by = c("Signature", "group1" = "Group_comp")) %>%
    dplyr::rename(group1_mean = group_mean) %>% 
    left_join(scores_group_means, by = c("Signature", "group2" = "Group_comp")) %>%
    dplyr::rename(group2_mean = group_mean)

write_tsv(scores_tidy_ttest, glue("{out}/ssgsea_scores_tests.tsv"))

Data wrangling to make sure the fold-changes are in the same direction for each comparison:

scores_tidy_ttest_foxr2 <- scores_tidy_ttest %>% 
    filter(group1 == "NB-FOXR2" | group2 == "NB-FOXR2") %>% 
    mutate(log10p = -log10(p_adj),
           log2fc = log2(group1_mean / group2_mean)) %>%
    # put all group1 as NB-FOXR2, and reverse L2FC if needed
    # this is because the pairwise t-tests do not perform two tests for each
    # pair of groups (i.e. only A vs B, not A vs B and B vs A)
    mutate(
        group_tmp = group1,
        # if the Group2 is NB-FOXR2, then the LFC sign needs to be switched
        log2fc = ifelse(group2 == "NB-FOXR2", -log2fc, log2fc),
        # if the Group2 is NB-FOXR2, swap group1 and 2
        group1 = ifelse(group2 == "NB-FOXR2", "NB-FOXR2", group1),
        group2 = ifelse(group2 == "NB-FOXR2", group_tmp, group2)
    ) %>% 
    select(-group_tmp) %>% 
    mutate(signif = ifelse(p_adj < 0.01, TRUE, FALSE)) %>% 
    left_join(cell_type_anno_all, by = c("Signature" = "Cluster"))

Visualize the results as volcano plots:

df1 <- scores_tidy_ttest_foxr2 %>% 
    filter(group2 == "Neuronal")

p1 <- df1 %>%
    mutate(Class = ifelse(log2fc < 0 | p_adj > 0.01, "Other", Class)) %>%
    ggplot(aes(x = log2fc, y = -log10(p_adj))) +
    geom_hline(yintercept = 2, color = "gray50", linetype = 2) +
    geom_vline(xintercept = 0, color = "gray50", linetype = 2) +
    geom_point(aes(color = Class, size = -log10(p_adj)), alpha = 0.8) +
    scale_color_manual(values = palette_class) +
    no_legend() +
    ggtitle("NB-FOXR2 vs neuronal tumors") +
    coord_cartesian(xlim = c(-0.3, 0.3))

df2 <- scores_tidy_ttest_foxr2 %>%
    filter(group2 == "Glial")

p2 <- df2 %>%
    mutate(Class = ifelse(log2fc < 0 | p_adj > 0.01, "Other", Class)) %>%
    ggplot(aes(x = log2fc, y = -log10(p_adj))) +
    geom_hline(yintercept = 2, color = "gray50", linetype = 2) +
    geom_vline(xintercept = 0, color = "gray50", linetype = 2) +
    geom_point(aes(color = Class, size = -log10(p_adj)), alpha = 0.8) +
    scale_color_manual(values = palette_class) +
    no_legend() +
    ggtitle("NB-FOXR2 vs glial tumors") +
    coord_cartesian(xlim = c(-0.3, 0.3))

plot_grid(p1, p2, nrow = 1)

# with signature labels
p1 <- p1 + geom_text_repel(data = df1 %>% filter(-log10(p_adj) > 5 & log2fc > 0),
                    aes(label = Signature, color = Class),
                    size = 4)
p2 <- p2 + geom_text_repel(data = df2 %>% filter(-log10(p_adj) > 10 & log2fc > 0),
                    aes(label = Signature, color = Class),
                    size = 4)

plot_grid(p1, p2, nrow = 1)


14 Reproducibility

This document was last rendered on:

## 2024-11-05 10:36:03

The git repository and last commit:

## Local:    main /project/kleinman/bhavyaa.chandarana/from_hydra/2023-05-NB-FOXR2/public
## Remote:   main @ origin (https://github.com/fungenomics/NB-FOXR2.git)
## Head:     [72d1c5c] 2024-11-04: Add DOI badge

The random seed was set with set.seed(100)

The R session info:

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.1.2 (2021-11-01)
##  os       Rocky Linux 8.10 (Green Obsidian)
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/Toronto
##  date     2024-11-05
##  pandoc   1.19.2.1 @ /cvmfs/soft.computecanada.ca/gentoo/2020/usr/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  ! package        * version date (UTC) lib source
##  P abind            1.4-5   2016-07-21 [?] CRAN (R 4.1.2)
##  P backports        1.4.1   2021-12-13 [?] CRAN (R 4.1.2)
##  P BiocGenerics     0.40.0  2021-10-26 [?] Bioconductor
##  P BiocManager      1.30.15 2021-05-11 [?] CRAN (R 4.1.2)
##  P bit              4.0.4   2020-08-04 [?] CRAN (R 4.1.2)
##  P bit64            4.0.5   2020-08-30 [?] CRAN (R 4.1.2)
##  P broom          * 1.0.1   2022-08-29 [?] CRAN (R 4.1.2)
##  P bslib            0.3.1   2021-10-06 [?] CRAN (R 4.1.2)
##  P cachem           1.0.6   2021-08-19 [?] CRAN (R 4.1.2)
##  P callr            3.7.6   2024-03-25 [?] RSPM
##  P cellranger       1.1.0   2016-07-27 [?] CRAN (R 4.1.2)
##  P circlize         0.4.15  2022-05-10 [?] CRAN (R 4.1.2)
##  P cli              3.6.1   2023-03-23 [?] RSPM (R 4.1.2)
##  P clue             0.3-64  2023-01-31 [?] CRAN (R 4.1.2)
##  P cluster          2.1.2   2021-04-17 [?] CRAN (R 4.1.2)
##  P codetools        0.2-18  2020-11-04 [?] CRAN (R 4.1.2)
##  P colorspace       2.0-2   2021-06-24 [?] CRAN (R 4.1.2)
##  P ComplexHeatmap * 2.10.0  2021-10-26 [?] Bioconductor
##  P conflicted       1.2.0   2023-02-01 [?] CRAN (R 4.1.2)
##  P cowplot        * 1.1.1   2020-12-30 [?] CRAN (R 4.1.2)
##  P crayon           1.4.2   2021-10-29 [?] CRAN (R 4.1.2)
##  P crosstalk        1.2.0   2021-11-04 [?] CRAN (R 4.1.2)
##  P data.table     * 1.14.2  2021-09-27 [?] CRAN (R 4.1.2)
##  P deldir           1.0-6   2021-10-23 [?] CRAN (R 4.1.2)
##  P devtools         2.4.5   2022-10-11 [?] CRAN (R 4.1.2)
##  P digest           0.6.35  2024-03-11 [?] CRAN (R 4.1.2)
##  P doParallel       1.0.16  2020-10-16 [?] CRAN (R 4.1.2)
##  P dplyr          * 1.1.1   2023-03-22 [?] CRAN (R 4.1.2)
##  P DT               0.30    2023-10-05 [?] CRAN (R 4.1.2)
##  P ellipsis         0.3.2   2021-04-29 [?] CRAN (R 4.1.2)
##  P evaluate         0.23    2023-11-01 [?] CRAN (R 4.1.2)
##  P fansi            1.0.2   2022-01-14 [?] CRAN (R 4.1.2)
##  P farver           2.1.0   2021-02-28 [?] CRAN (R 4.1.2)
##  P fastmap          1.1.0   2021-01-25 [?] CRAN (R 4.1.2)
##  P fitdistrplus     1.1-6   2021-09-28 [?] CRAN (R 4.1.2)
##  P foreach          1.5.1   2020-10-15 [?] CRAN (R 4.1.2)
##  P fs               1.5.2   2021-12-08 [?] CRAN (R 4.1.2)
##  P future           1.25.0  2022-04-24 [?] CRAN (R 4.1.2)
##  P future.apply     1.8.1   2021-08-10 [?] CRAN (R 4.1.2)
##  P generics         0.1.3   2022-07-05 [?] CRAN (R 4.1.2)
##  P GetoptLong       1.0.5   2020-12-15 [?] CRAN (R 4.1.2)
##  P ggExtra        * 0.10.0  2022-03-23 [?] CRAN (R 4.1.2)
##  P ggplot2        * 3.4.2   2023-04-03 [?] CRAN (R 4.1.2)
##  P ggrepel        * 0.9.1   2021-01-15 [?] CRAN (R 4.1.2)
##  P ggridges         0.5.3   2021-01-08 [?] CRAN (R 4.1.2)
##  P git2r            0.29.0  2021-11-22 [?] CRAN (R 4.1.2)
##  P GlobalOptions    0.1.2   2020-06-10 [?] CRAN (R 4.1.2)
##  P globals          0.14.0  2020-11-22 [?] CRAN (R 4.1.2)
##  P glue           * 1.6.2   2022-02-24 [?] CRAN (R 4.1.2)
##  P goftest          1.2-3   2021-10-07 [?] CRAN (R 4.1.2)
##  P gridExtra        2.3     2017-09-09 [?] CRAN (R 4.1.2)
##  P gtable           0.3.0   2019-03-25 [?] CRAN (R 4.1.2)
##  P here           * 1.0.1   2020-12-13 [?] CRAN (R 4.1.2)
##  P highr            0.9     2021-04-16 [?] CRAN (R 4.1.2)
##  P hms              1.1.1   2021-09-26 [?] CRAN (R 4.1.2)
##  P htmltools        0.5.2   2021-08-25 [?] CRAN (R 4.1.2)
##  P htmlwidgets      1.5.4   2021-09-08 [?] CRAN (R 4.1.2)
##  P httpuv           1.6.5   2022-01-05 [?] CRAN (R 4.1.2)
##  P httr             1.4.2   2020-07-20 [?] CRAN (R 4.1.2)
##  P ica              1.0-2   2018-05-24 [?] CRAN (R 4.1.2)
##  P igraph           2.0.3   2024-03-13 [?] CRAN (R 4.1.2)
##  P IRanges          2.28.0  2021-10-26 [?] Bioconductor
##  P irlba            2.3.5   2021-12-06 [?] CRAN (R 4.1.2)
##  P iterators        1.0.13  2020-10-15 [?] CRAN (R 4.1.2)
##  P jquerylib        0.1.4   2021-04-26 [?] CRAN (R 4.1.2)
##  P jsonlite         1.8.8   2023-12-04 [?] CRAN (R 4.1.2)
##  P KernSmooth       2.23-20 2021-05-03 [?] CRAN (R 4.1.2)
##  P knitr            1.37    2021-12-16 [?] CRAN (R 4.1.2)
##  P labeling         0.4.2   2020-10-20 [?] CRAN (R 4.1.2)
##  P later            1.3.0   2021-08-18 [?] CRAN (R 4.1.2)
##  P lattice          0.20-45 2021-09-22 [?] CRAN (R 4.1.2)
##  P lazyeval         0.2.2   2019-03-15 [?] CRAN (R 4.1.2)
##  P leiden           0.3.9   2021-07-27 [?] CRAN (R 4.1.2)
##  P lifecycle        1.0.3   2022-10-07 [?] CRAN (R 4.1.2)
##  P listenv          0.8.0   2019-12-05 [?] CRAN (R 4.1.2)
##  P lmtest           0.9-39  2021-11-07 [?] CRAN (R 4.1.2)
##  P magrittr       * 2.0.3   2022-03-30 [?] CRAN (R 4.1.2)
##  P MASS             7.3-54  2021-05-03 [?] CRAN (R 4.1.2)
##  P Matrix           1.3-4   2021-06-01 [?] CRAN (R 4.1.2)
##  P matrixStats      0.61.0  2021-09-17 [?] CRAN (R 4.1.2)
##  P memoise          2.0.1   2021-11-26 [?] CRAN (R 4.1.2)
##  P mgcv             1.8-38  2021-10-06 [?] CRAN (R 4.1.2)
##  P mime             0.12    2021-09-28 [?] CRAN (R 4.1.2)
##  P miniUI           0.1.1.1 2018-05-18 [?] CRAN (R 4.1.2)
##  P munsell          0.5.0   2018-06-12 [?] CRAN (R 4.1.2)
##  P nlme             3.1-153 2021-09-07 [?] CRAN (R 4.1.2)
##  P parallelly       1.30.0  2021-12-17 [?] CRAN (R 4.1.2)
##  P patchwork        1.1.1   2020-12-17 [?] CRAN (R 4.1.2)
##  P pbapply        * 1.5-0   2021-09-16 [?] CRAN (R 4.1.2)
##  P pillar           1.9.0   2023-03-22 [?] RSPM (R 4.1.2)
##  P pkgbuild         1.4.2   2023-06-26 [?] CRAN (R 4.1.2)
##  P pkgconfig        2.0.3   2019-09-22 [?] CRAN (R 4.1.2)
##  P pkgload          1.3.3   2023-09-22 [?] CRAN (R 4.1.2)
##  P plotly           4.10.0  2021-10-09 [?] CRAN (R 4.1.2)
##  P plyr             1.8.6   2020-03-03 [?] CRAN (R 4.1.2)
##  P png              0.1-7   2013-12-03 [?] CRAN (R 4.1.2)
##  P polyclip         1.10-0  2019-03-14 [?] CRAN (R 4.1.2)
##  P prettyunits      1.1.1   2020-01-24 [?] CRAN (R 4.1.2)
##  P processx         3.8.4   2024-03-16 [?] RSPM
##  P profvis          0.3.8   2023-05-02 [?] CRAN (R 4.1.2)
##  P promises         1.2.0.1 2021-02-11 [?] CRAN (R 4.1.2)
##  P ps               1.7.6   2024-01-18 [?] RSPM
##  P purrr          * 1.0.1   2023-01-10 [?] CRAN (R 4.1.2)
##  P R6               2.5.1   2021-08-19 [?] CRAN (R 4.1.2)
##  P RANN             2.6.1   2019-01-08 [?] CRAN (R 4.1.2)
##  P RColorBrewer   * 1.1-2   2014-12-07 [?] CRAN (R 4.1.2)
##  P Rcpp             1.0.8   2022-01-13 [?] CRAN (R 4.1.2)
##  P RcppAnnoy        0.0.19  2021-07-30 [?] CRAN (R 4.1.2)
##  P readr          * 2.1.1   2021-11-30 [?] CRAN (R 4.1.2)
##  P readxl         * 1.3.1   2019-03-13 [?] CRAN (R 4.1.2)
##  P remotes          2.4.2.1 2023-07-18 [?] CRAN (R 4.1.2)
##  P renv             1.0.3   2023-09-19 [?] CRAN (R 4.1.2)
##  P reshape2         1.4.4   2020-04-09 [?] CRAN (R 4.1.2)
##  P reticulate       1.23    2022-01-14 [?] CRAN (R 4.1.2)
##  P rjson            0.2.21  2022-01-09 [?] CRAN (R 4.1.2)
##  P rlang            1.1.3   2024-01-10 [?] CRAN (R 4.1.2)
##  P rmarkdown        2.11    2021-09-14 [?] CRAN (R 4.1.2)
##  P ROCR             1.0-11  2020-05-02 [?] CRAN (R 4.1.2)
##  P rpart            4.1-15  2019-04-12 [?] CRAN (R 4.1.2)
##  P rprojroot        2.0.2   2020-11-15 [?] CRAN (R 4.1.2)
##  P Rtsne          * 0.15    2018-11-10 [?] CRAN (R 4.1.2)
##  P S4Vectors        0.32.4  2022-03-24 [?] Bioconductor
##  P sass             0.4.0   2021-05-12 [?] CRAN (R 4.1.2)
##  P scales           1.2.1   2022-08-20 [?] CRAN (R 4.1.2)
##  P scattermore      0.7     2020-11-24 [?] CRAN (R 4.1.2)
##  P sctransform      0.3.3   2022-01-13 [?] CRAN (R 4.1.2)
##  P sessioninfo      1.2.2   2021-12-06 [?] CRAN (R 4.1.2)
##  P Seurat         * 4.0.0   2021-01-30 [?] CRAN (R 4.1.2)
##  P SeuratObject   * 4.0.4   2021-11-23 [?] CRAN (R 4.1.2)
##  P shape            1.4.6   2021-05-19 [?] CRAN (R 4.1.2)
##  P shiny            1.7.1   2021-10-02 [?] CRAN (R 4.1.2)
##  P spatstat         1.64-1  2020-05-12 [?] CRAN (R 4.1.2)
##  P spatstat.data    2.1-2   2021-12-17 [?] CRAN (R 4.1.2)
##  P spatstat.utils   2.3-0   2021-12-12 [?] CRAN (R 4.1.2)
##  P stringi          1.7.6   2021-11-29 [?] CRAN (R 4.1.2)
##  P stringr        * 1.5.0   2022-12-02 [?] CRAN (R 4.1.2)
##  P survival         3.2-13  2021-08-24 [?] CRAN (R 4.1.2)
##  P tensor           1.5     2012-05-05 [?] CRAN (R 4.1.2)
##  P tibble         * 3.2.1   2023-03-20 [?] RSPM (R 4.1.2)
##  P tidyr          * 1.3.0   2023-01-24 [?] CRAN (R 4.1.2)
##  P tidyselect       1.2.0   2022-10-10 [?] CRAN (R 4.1.2)
##  P tzdb             0.3.0   2022-03-28 [?] CRAN (R 4.1.2)
##  P urlchecker       1.0.1   2021-11-30 [?] CRAN (R 4.1.2)
##  P usethis          2.2.2   2023-07-06 [?] CRAN (R 4.1.2)
##  P utf8             1.2.2   2021-07-24 [?] CRAN (R 4.1.2)
##  P uwot             0.1.11  2021-12-02 [?] CRAN (R 4.1.2)
##  P vctrs            0.6.5   2023-12-01 [?] CRAN (R 4.1.2)
##  P viridis        * 0.5.1   2018-03-29 [?] RSPM (R 4.1.2)
##  P viridisLite    * 0.3.0   2018-02-01 [?] CRAN (R 4.1.2)
##  P vroom            1.5.7   2021-11-30 [?] CRAN (R 4.1.2)
##  P withr            2.5.0   2022-03-03 [?] CRAN (R 4.1.2)
##  P xfun             0.29    2021-12-14 [?] CRAN (R 4.1.2)
##  P xtable           1.8-4   2019-04-21 [?] CRAN (R 4.1.2)
##  P yaml             2.2.1   2020-02-01 [?] CRAN (R 4.1.2)
##  P zoo              1.8-9   2021-03-09 [?] CRAN (R 4.1.2)
## 
##  [1] /project/kleinman/bhavyaa.chandarana/from_hydra/2023-05-NB-FOXR2/public/renv/library/R-4.1/x86_64-pc-linux-gnu
##  [2] /home/kleinman/bhavyaa.chandarana/.cache/R/renv/sandbox/R-4.1/x86_64-pc-linux-gnu/145cef2c
## 
##  P ── Loaded and on-disk path mismatch.
## 
## ──────────────────────────────────────────────────────────────────────────────


A project of the Kleinman Lab at McGill University, using the rr reproducible research template.