4.1 Cell type density along normal interneuron differentiation trajectory

First we load a table of pseudotime values per cell:

# Clean up the cell names, which got mangled a bit due to additional prefixes
# when the Seurat objects were merged
tidy_cellnames <- function(pseudotime) {
    
    pseudotime %>%
        mutate(cell = case_when(
            grepl("e12|e15", cell) ~ paste0("__", cell),
            grepl("p0", cell) ~ paste0("__", cell),
            TRUE ~ cell
        )) %>%
        mutate(pseudotime = Pseudotime)
    
}

inhibitory_pseudotime <- read_tsv(here("reference_datasets/2019_Jessa/Jessa2019_inhibitory_neurons_pseudotime.txt")) %>%
    tidy_cellnames()

## Parsed with column specification:
## cols(
##   cell = col_character(),
##   Pseudotime = col_double()
## )

Next load the cell type identities for each cell:

inhibitory_pseudotime$pseudotime %<>% scales::rescale(to = c(0, 100))

# Load the metadata for the joint forebrain and set the identities
forebrain_meta <- read_tsv(here("reference_datasets/2019_Jessa/Jessa2019_joint_forebrain_metadata.tsv"))

## Parsed with column specification:
## cols(
##   Cell = col_character(),
##   nGene = col_double(),
##   nUMI = col_double(),
##   orig.ident = col_character(),
##   percent_mito = col_double(),
##   res.0.8 = col_double(),
##   ID_20190715_with_blacklist_and_refined = col_character(),
##   ID_20190715_joint_clustering = col_character(),
##   ID_20190730_with_blacklist_and_refined = col_character(),
##   Cluster = col_double()
## )

inhibitory_pseudotime <- left_join(inhibitory_pseudotime, forebrain_meta, by = c("cell" = "Cell"))

Finally, generate a density plot of each cell type along pseudotime, to represent the relative proportions of each cell type as differentiation proceeds:

# Grab cells in the clusters included in this Monocle analysis,
# and plot their density in pseuodtime
inhibitory_pseudotime %>%
    filter(grepl("VRGC|INIP|MGINH|CINH", ID_20190730_with_blacklist_and_refined)) %>%
    separate(ID_20190730_with_blacklist_and_refined,
             sep = "_",
             into = c("Sample", "Cell_type")) %>%
    rr_ggplot(aes(x = Pseudotime), plot_num = 1) +
    geom_density(aes(fill = Cell_type), alpha = 0.8, size = 0.4) +
    scale_fill_manual(values = c(
        "VRGC"   = "#ffb219",
        "INIP-P" = "#8798C8",
        "MGINH"  = "#4a1cc9",
        "CINHN"  = "#135ca0"
    )) +
    theme_min()

_{[figure/source data @ G34-gliomas/singlecell_normal/figures/01//interneuron_differentiation_density_plot…]}

4.2 Plot genes of interest along pseudotime

load(here("reference_datasets/2019_Jessa/Jessa2019_inhibitory_neurons.monocle.Rda"))

pal_pseudotime <- colorRampPalette(c("navy", "white", "red"))(n = 200)

One way to evaluate whether the genes which drive enrichment of fetal neuronal signatures vs. depletion of pediatric/adult signatures in G34 gliomas by GSEA are related to maturity is by checking their expression during normal differentiation. Get the genes from the GSEA analysis:

signatures <- c("MGE newborn neurons 1", "Parvalbumin interneurons")

fgsea_df_idh_filt <- read_tsv(here("bulk_transcriptome_epigenome/output/02/fgsea_results_padj<0.01_IDH_top_le.tsv"))

## Parsed with column specification:
## cols(
##   Comparison = col_character(),
##   Signature = col_character(),
##   Sample = col_character(),
##   Age = col_character(),
##   Cell_type = col_character(),
##   pval = col_double(),
##   padj = col_double(),
##   ES = col_double(),
##   NES = col_double(),
##   nMoreExtreme = col_double(),
##   size = col_double(),
##   leadingEdge = col_character(),
##   Dataset = col_character(),
##   leadingEdge_top = col_character()
## )

leading_edge <- fgsea_df_idh_filt %>% 
    filter(Cell_type %in% signatures) %>% 
    select(Cell_type, leadingEdge_top) %>% 
    tibble::deframe() %>% 
    map(~ .x %>% stringr::str_split(", ") %>% unlist() %>% hg2mm() %>% filter_NA())

leading_edge

## $`MGE newborn neurons 1`
##  [1] "Arx"           "Dlx6"          "Smarcd1"       "Plxna2"       
##  [5] "Dcx"           "St8sia5"       "E030030I06Rik" "Dlx5"         
##  [9] "Fam65b"        "Dlx2"          "Tmem2"         "Gad2"         
## [13] "Dlx1"          "Pdgfd"         "Metap1d"       "Arl4d"        
## [17] "Sox1"          "Ccser1"       
## 
## $`Parvalbumin interneurons`
##  [1] "Epha6"    "Fgf12"    "Nxph1"    "Gria4"    "Sparcl1"  "Cntnap5c"
##  [7] "Cntnap5b" "Cntnap5a" "Adcy8"    "Kcnab1"   "Frmd5"    "Plcxd3"  
## [13] "Caln1"    "Kcnd2"    "Dner"     "Vwc2"     "Fam19a2"  "Shisa9"  
## [19] "Ptchd4"

However, we only want to plot genes with reasonable mean exp in the clusters of this lineage, so that the smoothing makes sense:

all_genes <- read_tsv(here("reference_datasets/2019_Jessa/Jessa2019_forebrain_gene_names.tsv")) %>%
    pull(genes)

## Parsed with column specification:
## cols(
##   genes = col_character()
## )

filter_by_mean_expr <- function(genes) {
    
    genes <- genes[genes %in% all_genes] %>% unique()
    
    # Check expression in all the clusters which comprise the cortical
    # inhibitory neuron lineages
    feather::read_feather(here("reference_datasets/2019_Jessa/Jessa2019_mean_expression.feather"),
                          c("Cluster", genes)) %>%
        gather(Gene, Mean_expression, 2:ncol(.)) %>%
        filter(grepl("^F", Cluster) & grepl("e12|e15|p0", Cluster)) %>%
        filter(grepl("INIP|INH|VRGC", Cluster)) %>%
        filter(Mean_expression > 0.1) %>%
        pull(Gene) %>%
        unique() %>%
        sort()
    
}

(leading_edge_filt <- map(leading_edge, ~ filter_by_mean_expr(.x)))

## $`MGE newborn neurons 1`
##  [1] "Arl4d"   "Arx"     "Dcx"     "Dlx1"    "Dlx2"    "Dlx5"    "Dlx6"   
##  [8] "Fam65b"  "Gad2"    "Metap1d" "Plxna2"  "Smarcd1" "Sox1"    "Tmem2"  
## 
## $`Parvalbumin interneurons`
## [1] "Caln1"   "Dner"    "Fam19a2" "Frmd5"   "Gria4"   "Kcnd2"   "Nxph1"  
## [8] "Plcxd3"  "Sparcl1"

names(leading_edge_filt) <- gsub(" ", "_", names(leading_edge_filt))

Extract genes to examine in pseudotime from the transcriptome/epigenome analyses – we’ll gate genes based on high zK27Ac and low zK27me3, or vice versa:

# Read in output of ChIP-seq analysis
multiomic_genes <- read_tsv(here("bulk_transcriptome_epigenome/output/03/DGE_and_histone_marks.tsv"))

## Parsed with column specification:
## cols(
##   ENSID = col_character(),
##   symbol = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double(),
##   G34R.zK27ac = col_double(),
##   G34R.zK27m3 = col_double(),
##   G34RV.K27ac.median = col_double(),
##   G34RV.K27me3.median = col_double(),
##   nonG34RV.K27ac.median = col_double(),
##   nonG34RV.K27me3.median = col_double()
## )

multiomic_genes_top <- list("UP" = multiomic_genes %>%
                             filter(G34R.zK27ac > 0.6) %>%
                             pull(symbol) %>% hg2mm() %>% filter_NA(),
                         "DOWN" = multiomic_genes %>%
                             filter(G34R.zK27m3 > 0.6) %>%
                             pull(symbol) %>% hg2mm() %>% filter_NA())

(multiomic_genes_top_filt <- map(multiomic_genes_top, ~ filter_by_mean_expr(.x)))

## $UP
##  [1] "Cdca7l"  "Chic2"   "Dlx1"    "Dlx2"    "Dlx6"    "Efnb1"   "Elovl2" 
##  [8] "Ephb2"   "Foxg1"   "Foxm1"   "Gad2"    "Gsx2"    "Ildr2"   "Insm1"  
## [15] "Jag1"    "Npy"     "Nr2e1"   "Otx1"    "Pax6"    "Plk2"    "Prdm16" 
## [22] "Sel1l3"  "Skida1"  "Slc10a4" "Sox1"    "Sp8"     "Tacc3"   "Tead2"  
## 
## $DOWN
##  [1] "Arhgdig" "Cadps"   "Chd5"    "Dlgap3"  "Igsf21"  "Jph3"    "Jph4"   
##  [8] "Npas1"   "Ntrk2"   "Nxph1"   "Sertad4" "Sez6l2"  "Shank1"  "St8sia1"

Intersect RNA-seq leading genes and multi-omic genes:

pseudotime_union <- list(
    # Newborn interneuron & upregulated genes
    "UP"   = unique(c(leading_edge_filt[[1]], multiomic_genes_top_filt[[1]])),
                         
    # Postnatal interneuron & downregulated genes
    "DOWN" = unique(c(leading_edge_filt[[2]], multiomic_genes_top_filt[[2]])))

imap(pseudotime_union, ~ ordered_pseudotime_heatmap(monocle[.x,],
                                                    show_rownames = TRUE,
                                                    order_by_pseudotime = TRUE,
                                                    hmcols = pal_pseudotime,
                                                    prefix = glue("{figout}/pseudotime_heatmap{.y}")))

## $UP
## 
## $DOWN

knitr::include_graphics(glue("{figout}/pseudotime_heatmapUP.png"))

knitr::include_graphics(glue("{figout}/pseudotime_heatmapDOWN.png"))

_{[figure/source data @ G34-gliomas/singlecell_normal/figures/01//pseudotime_heatmap…]}