Mouse Muscle Stem Cell Project

Part 6: compare with bulk rna-seq and run enrichment analysis

Génome Québec Innovation Centre, McGill University, Montreal, Canada

Set filepaths and parameters

set.seed(42)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
options(
  readr.show_progress = FALSE,
  digits = 2
)

Load packages

suppressPackageStartupMessages({
  library(scater)
  library(scran)
  library(tidyverse)
  library(pheatmap)
  library(hypeR)
  library(EnhancedVolcano)
  theme_set(theme_bw())
})

Define file paths

data_dir <- "./data"
figures_dir <- file.path("./figures")

plot_expression <- function(gene) {
  plotExpression(sce10x_stem,
    features = gene,
    x = "sample",
    exprs_values = "logcounts",
    colour_by = "condition",
    point_size = 1
  ) +
    facet_wrap(~ colData(sce10x_stem)$clusters, ncol = 1)
}

plot_volcano <- function(var_name, lfc_thresh, svalue_thresh, lfc, label = NULL) {
  EnhancedVolcano(lfc,
    lab = lfc %>% pull(gene_id),
    x = paste0("sc_", var_name),
    selectLab = label,
    y = paste0("svalue_sc_", var_name),
    xlab = bquote(~ italic(Moderated) ~ Log[2] ~ FC),
    ylab = bquote(~ -Log[10] ~ italic(svalue)),
    col = c("grey30", "forestgreen", "red2", "royalblue"),
    pCutoff = svalue_thresh,
    FCcutoff = lfc_thresh,
    legendLabels = c(
      "NS",
      expression(Log[2] ~ FC),
      "s-value",
      expression(s - value ~ and ~ log[2] ~ FC)
    )
  )
}

sc_lfc <-
  read_tsv(file.path(
    data_dir,
    "preprocessed",
    "stem_scrnaseq_moderated_lfc.txt"
  ))

Parsed with column specification:
cols(
  .default = col_double(),
  gene_id = col_character(),
  ensembl_gene_id_version = col_character(),
  external_gene_name = col_character(),
  gene_biotype = col_character(),
  description = col_character(),
  nonzero_stem = col_character(),
  keep_stem = col_logical(),
  top_hvg_stem = col_logical(),
  nonzero_fap = col_character(),
  keep_fap = col_logical(),
  top_hvg_fap = col_logical(),
  nonzero_macrophage = col_character(),
  keep_macrophage = col_logical(),
  top_hvg_macrophage = col_logical(),
  allZero = col_logical(),
  dispOutlier = col_logical(),
  fullBetaConv = col_logical(),
  reducedBetaConv = col_logical()
)
See spec(...) for full column specifications.
124 parsing failures.
 row col expected actual                                                  file
3023 chr a double      X './data/preprocessed/stem_scrnaseq_moderated_lfc.txt'
3024 chr a double      X './data/preprocessed/stem_scrnaseq_moderated_lfc.txt'
3025 chr a double      X './data/preprocessed/stem_scrnaseq_moderated_lfc.txt'
3026 chr a double      X './data/preprocessed/stem_scrnaseq_moderated_lfc.txt'
3027 chr a double      X './data/preprocessed/stem_scrnaseq_moderated_lfc.txt'
.... ... ........ ...... .....................................................
See problems(...) for more details.

sce10x_stem <-
  readRDS(
    file.path(
      data_dir,
      "preprocessed",
      "sce10x_stem_filtered_final.rds"
    )
  )

Combine with bulk rnaseq results

bulk_lfc <-
  read_tsv(file.path(
    data_dir,
    "preprocessed",
    "leg_trunk_moderated_lfc_1_svalue_05_marked.txt")) 

Parsed with column specification:
cols(
  .default = col_double(),
  gene_name = col_character(),
  gene_id = col_character(),
  EAN_leg = col_character(),
  keep_niche_leg = col_logical(),
  keep_age_leg = col_logical(),
  keep_engraf_leg = col_logical(),
  EAN_trunk = col_character(),
  keep_niche_trunk = col_logical(),
  keep_age_trunk = col_logical(),
  keep_engraf_trunk = col_logical(),
  seqnames = col_character(),
  gene_biotype = col_character(),
  description = col_character()
)
See spec(...) for full column specifications.

gene_metadata_filepath <-
  "~/Documents/genome_metadata/data/mm10_ens97_gene_meta.txt"
gene_metadata <- read_csv(gene_metadata_filepath)

Parsed with column specification:
cols(
  gene_name = col_character(),
  gene_len = col_double(),
  perc_gene_gc = col_double(),
  gene_id = col_character(),
  seqnames = col_character(),
  gene_biotype = col_character(),
  gene_id_version = col_character(),
  description = col_character()
)

bulk_lfc <-
  left_join(bulk_lfc, 
            gene_metadata %>% 
              select(gene_id, gene_id_version), 
            by="gene_id")  %>%
  rename(ensembl_gene_id_version="gene_id_version")

#write_tsv(bulk_lfc, file.path(
#    data_dir,
 #   "preprocessed",
 #   "leg_trunk_moderated_lfc_1_svalue_05_marked.txt"))

lfc <-
  left_join(sc_lfc, bulk_lfc %>%
    drop_na(EAN_leg)%>%
    select(c(ensembl_gene_id_version, ends_with("_leg"))),
  by = "ensembl_gene_id_version"
  )

moderated_lfc_leg_muscle <-
  lfc[, c(
    "gene_id", "baseMean",
    "lfc_stem_aging", "lfc_stem2_1", "lfc_stem3_12y",
    "lfc_engraf_leg", "lfc_age_leg", "lfc_niche_leg",
    "svalue_stem_aging", "svalue_stem2_1", "svalue_stem3_12y",
    "svalue_engraf_leg", "svalue_age_leg", "svalue_niche_leg",
    "description", "gene_biotype", "chr", "ensembl_gene_id_version"
  )] %>%
  set_names(c(
    "gene_id", "sc_mean",
    "sc_aging", "sc_clust2vs1", "sc_clust3vs12y",
    "bulk_engraf", "bulk_aging", "bulk_niche",
    "svalue_sc_aging", "svalue_sc_clust2vs1", "svalue_sc_clust3vs12y",
    "svalue_bulk_engraf", "svalue_bulk_aging", "svalue_bulk_niche",
    "description", "gene_biotype", "chr", "ensembl_gene_id_version"
  )) %>%
  arrange(-sc_mean)
moderated_lfc_leg_muscle

write_tsv(
  moderated_lfc_leg_muscle,
  file.path(
    data_dir,
    "preprocessed",
    "moderated_lfc_leg_muscle.txt"
  )
)

cor(moderated_lfc_leg_muscle %>%
  select(c("sc_aging", "bulk_aging")) %>%
  drop_na(),
method = "kendall"
)

           sc_aging bulk_aging
sc_aging       1.00       0.27
bulk_aging     0.27       1.00

p1 <-
  ggplot(
    moderated_lfc_leg_muscle,
    aes(
      x = sc_aging,
      y = bulk_aging,
      color = bulk_niche
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.9
  ) +
  geom_smooth() +
  coord_fixed(ratio = 1 / 2) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  scale_color_gradient2(
    midpoint = 2,
    low = "orange",
    mid = "black",
    high = "blue"
  ) +
  theme_classic() +
  xlim(-6, 3)

p1

ggsave(file.path(figures_dir, "scatterplot_bulk_vs_singlecell_aging_effect_labeled.pdf"), p1)

Saving 7 x 7 in image

p2 <-
  ggplot(
    moderated_lfc_leg_muscle,
    aes(
      x = sc_aging,
      y = sc_clust2vs1,
      color = sc_clust3vs12y
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  scale_color_gradient2(
    midpoint = 0,
    low = "yellow",
    mid = "black",
    high = "blue"
  ) +
  theme_classic()
p2

ggsave(file.path(figures_dir, "scatterplot_stem_aging_vs_cluster_effect_labeled.pdf"), p2)

Saving 7 x 7 in image

p2 <-
  ggplot(
    moderated_lfc_leg_muscle,
    aes(
      x = sc_aging,
      y = sc_clust2vs1,
      color = sc_clust3vs12y
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  scale_color_gradient2(
    midpoint = 0,
    low = "yellow",
    mid = "black",
    high = "blue"
  ) +
  theme_classic() +
  xlim(c(-6, 3))
p2

ggsave(file.path(figures_dir, "scatterplot_stem_aging_vs_cluster_effect_unlabeled.pdf"), p2)

Saving 7 x 7 in image

p3 <-
  ggplot(
    moderated_lfc_leg_muscle,
    aes(
      x = sc_clust2vs1,
      y = sc_clust3vs12y,
      color = sc_aging
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  scale_color_gradient2(
    midpoint = 0,
    limits = c(-3, 2),
    low = "yellow",
    mid = "black",
    high = "blue"
  ) +
  theme_classic() +
  ylim(c(-2.2, 1.8))
p3

ggsave(file.path(figures_dir, "scatterplot_stem3_vs_stem21_cluster_effect_labeled.pdf"), p3)

Saving 7 x 7 in image

p3 <-
  ggplot(
    moderated_lfc_leg_muscle,
    aes(
      x = sc_clust2vs1,
      y = sc_clust3vs12y,
      color = sc_aging
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  scale_color_gradient2(
    midpoint = 0,
    limits = c(-3, 2),
    low = "yellow",
    mid = "black",
    high = "blue"
  ) +
  theme_classic() +
  ylim(c(-2.2, 1.8))
p3

ggsave(file.path(figures_dir, "scatterplot_stem3_vs_stem21_effect_unlabeled.pdf"), p3)

Saving 7 x 7 in image

p4 <-
  ggplot(
    lfc,
    aes(
      x = lfc_stem1_aging,
      y = lfc_stem2_aging
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 4
  ) +
  theme_classic()
p4

ggsave(file.path(figures_dir, "scatterplot_stem_1vs2_aging_effect_labeled.pdf"), p4)

Saving 7 x 7 in image

p4 <-
  ggplot(
    lfc,
    aes(
      x = lfc_stem1_aging,
      y = lfc_stem2_aging
    )
  ) +
  geom_point(
    size = .8,
    alpha = 0.8
  ) +
  theme_classic() +
  xlim(c(-6, 3)) +
  ylim(c(-6, 3))
p4

ggsave(file.path(figures_dir, "scatterplot_stem_1vs2_aging_effect_unlabeled.pdf"), p4)

Saving 7 x 7 in image

Make volcano plots

lfc_thresh <- 1
svalue_thresh <- 10e-8
volcano_plots <-
  map(c("aging", "clust2vs1", "clust3vs12y"),
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = moderated_lfc_leg_muscle
  )
volcano_plots

[[1]]

[[2]]

[[3]]

ggsave(file.path(figures_dir, "volcano_stem_aging_effect_labeled.pdf"), volcano_plots[[1]])

Saving 7 x 7 in image

ggsave(file.path(figures_dir, "volcano_stem_cluster2vs1_labeled.pdf"), volcano_plots[[2]])
ggsave(file.path(figures_dir, "volcano_stem_cluster3vs1_2_labeled.pdf"), volcano_plots[[3]])

volcano_plots_unlabeled <-
  map(c("aging", "clust2vs1", "clust3vs12y"),
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = moderated_lfc_leg_muscle,
    label = c("")
  )
volcano_plots_unlabeled

[[1]]

[[2]]

[[3]]

ggsave(file.path(figures_dir, "volcano_stem_aging_effect_unlabeled.pdf"), volcano_plots_unlabeled[[1]])

Saving 7 x 7 in image

ggsave(file.path(figures_dir, "volcano_stem_cluster2vs1_unlabeled.pdf"), volcano_plots_unlabeled[[2]])
ggsave(file.path(figures_dir, "volcano_stem_cluster3vs1_2_unlabeled.pdf"), volcano_plots_unlabeled[[3]])

Histograms

ggplot(moderated_lfc_leg_muscle, aes(
  x = sc_aging
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

ggplot(moderated_lfc_leg_muscle, aes(
  x = sc_clust2vs1
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

ggplot(moderated_lfc_leg_muscle, aes(
  x = sc_clust3vs12y
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

ggplot(moderated_lfc_leg_muscle, aes(
  x = bulk_engraf
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

ggplot(moderated_lfc_leg_muscle, aes(
  x = bulk_aging
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

ggplot(moderated_lfc_leg_muscle, aes(
  x = bulk_niche
)) +
  geom_histogram(bins = 200) +
  geom_vline(xintercept = c(-1, 1))

Correlation analysis

cor_pairs <-
  correlatePairs(sce10x_stem,
    block = factor(colData(sce10x_stem)$sample)
  ) %>%
  as_tibble()
cor_pairs

Clusters 2 vs 1

genes_hits_clust2vs1_up <-
  moderated_lfc_leg_muscle %>%
  filter(sc_clust2vs1 > .75 & svalue_sc_clust2vs1 < 10^-8) %>%
  arrange(-sc_clust2vs1) %>%
  pull(gene_id)
genes_hits_clust2vs1_up

[1] "ENSMUSG00000031842.14" "Fos"                   "Gm26802"               "Jun"                   "Egr1"                  "ENSMUSG00000052837.6" 
[7] "Ier2"                  "Socs3"                 "Ier5l"                

map(genes_hits_clust2vs1_up, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

genes_hits_clust2vs1_down <-
  moderated_lfc_leg_muscle %>%
  filter(sc_clust2vs1 < -0.75 & svalue_sc_clust2vs1 < 10^-8) %>%
  arrange(sc_clust2vs1) %>%
  pull(gene_id)
genes_hits_clust2vs1_down

 [1] "Col8a1"                      "Iigp1"                       "Pakap_ENSMUSG00000038729.24" "Nfkb1"                      
 [5] "Mt1"                         "Dnah7a"                      "Mpp7"                        "9030624G23Rik"              
 [9] "Arid5b"                      "Samd4"                       "Asb5"                        "Mt2"                        
[13] "Emp1"                        "Gm48228"                     "Pmepa1"                      "Sdc4"                       
[17] "Runx1"                       "Gal"                        

map(genes_hits_clust2vs1_down, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

genes_hits_clust2vs1 <- c(genes_hits_clust2vs1_down, genes_hits_clust2vs1_up)

cor_pairs_clust2vs1 <-
  cor_pairs %>%
  filter((gene1 %in% genes_hits_clust2vs1 & gene2 %in% genes_hits_clust2vs1) & abs(rho) >= .3) %>%
  arrange(gene1, rho)
cor_pairs_clust2vs1

cor_pairs_clust2vs1_genes <- sort(union(cor_pairs_clust2vs1$gene1, cor_pairs_clust2vs1$gene2))
cor_pairs_clust2vs1_genes

 [1] "9030624G23Rik"               "Arid5b"                      "Egr1"                        "ENSMUSG00000052837.6"       
 [5] "Fos"                         "Gm48228"                     "Ier2"                        "Jun"                        
 [9] "Mt1"                         "Mt2"                         "Nfkb1"                       "Pakap_ENSMUSG00000038729.24"
[13] "Pmepa1"                      "Runx1"                       "Sdc4"                        "Socs3"                      

setdiff(genes_hits_clust2vs1, cor_pairs_clust2vs1_genes)

 [1] "Col8a1"                "Iigp1"                 "Dnah7a"                "Mpp7"                  "Samd4"                 "Asb5"                 
 [7] "Emp1"                  "Gal"                   "ENSMUSG00000031842.14" "Gm26802"               "Ier5l"                

plotExpression(sce10x_stem,
  features = setdiff(cor_pairs_clust2vs1_genes, "Fos"),
  x = c("Fos"),
  exprs_values = "logcounts",
  colour_by = "clusters",
  point_size = .7,
  point_alpha = 0.7,
  ncol = 4
) + coord_fixed()

plotExpression(sce10x_stem,
  features = setdiff(genes_hits_clust2vs1, cor_pairs_clust2vs1_genes),
  x = c("Fos"),
  exprs_values = "logcounts",
  colour_by = "clusters",
  point_size = .7,
  point_alpha = 0.7,
  ncol = 4
) + coord_fixed()

Clusters 3 vs 1 and 2

genes_hits_clust3vs12y_up <-
  moderated_lfc_leg_muscle %>%
  filter(sc_clust3vs12y > 1 & svalue_sc_clust3vs12y < 10^-8) %>%
  arrange(-sc_clust3vs12y) %>%
  pull(gene_id)
genes_hits_clust3vs12y_up

 [1] "Myog"    "Pclaf"   "Cdkn1c"  "Ccnb2"   "Cenpa"   "Stmn1"   "Acta2"   "Lgals1"  "Hmgn2"   "H2az1"   "Ctnna2"  "Cfap77"  "Ankrd10" "Lsm6"   

map(genes_hits_clust3vs12y_up, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

genes_hits_clust3vs12y_down <-
  moderated_lfc_leg_muscle %>%
  filter(sc_clust3vs12y < -1 & svalue_sc_clust3vs12y < 10^-8) %>%
  arrange(sc_clust3vs12y) %>%
  pull(gene_id)
genes_hits_clust3vs12y_down

 [1] "Hs6st3"   "Igfbp5"   "Gpx3"     "Nppc"     "Ltbp1"    "Apoe"     "Chodl"    "Crlf1"    "Mt1"      "Pdzd2"    "Sdc4"     "Kcnma1"   "Tnxb"    
[14] "Meg3"     "Rora"     "Igfbp4"   "Serping1" "Tmtc2"   

map(genes_hits_clust3vs12y_down, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

genes_hits_clust3vs12y <- c(genes_hits_clust3vs12y_up, genes_hits_clust3vs12y_down)

cor_pairs_clust3vs12y <-
  cor_pairs %>%
  filter((gene1 %in% genes_hits_clust3vs12y &
    gene2 %in% genes_hits_clust3vs12y) &
    abs(rho) >= .3) %>%
  arrange(gene1, rho)
cor_pairs_clust3vs12y

cor_pairs_clust3vs12y_genes <-
  sort(union(
    cor_pairs_clust3vs12y$gene1,
    cor_pairs_clust3vs12y$gene2
  ))
cor_pairs_clust3vs12y_genes

[1] "Chodl" "Crlf1" "Gpx3"  "Mt1"   "Sdc4" 

setdiff(genes_hits_clust3vs12y, cor_pairs_clust3vs12y_genes)

 [1] "Myog"     "Pclaf"    "Cdkn1c"   "Ccnb2"    "Cenpa"    "Stmn1"    "Acta2"    "Lgals1"   "Hmgn2"    "H2az1"    "Ctnna2"   "Cfap77"   "Ankrd10" 
[14] "Lsm6"     "Hs6st3"   "Igfbp5"   "Nppc"     "Ltbp1"    "Apoe"     "Pdzd2"    "Kcnma1"   "Tnxb"     "Meg3"     "Rora"     "Igfbp4"   "Serping1"
[27] "Tmtc2"   

plotExpression(sce10x_stem,
  features = setdiff(cor_pairs_clust3vs12y_genes, "Gpx3"),
  x = c("Gpx3"),
  exprs_values = "logcounts",
  colour_by = "clusters",
  point_size = .7,
  point_alpha = 0.8,
  ncol = 4
) + coord_fixed()

plotExpression(sce10x_stem,
  features = setdiff(
    genes_hits_clust3vs12y,
    cor_pairs_clust3vs12y_genes
  ),
  x = c("Gpx3"),
  exprs_values = "logcounts",
  colour_by = "clusters",
  point_size = .7,
  point_alpha = 0.8,
  ncol = 4
) + coord_fixed()

Aged vs young

genes_hits_aging_up <-
  moderated_lfc_leg_muscle %>%
  filter((sc_aging > 1 & svalue_sc_aging < 10^-8) |
    (sc_aging > 0.5 & svalue_sc_aging < 10^-8 & bulk_aging > 2)) %>%
  arrange(-sc_aging) %>%
  pull(gene_id)
genes_hits_aging_up

 [1] "Gm7536"        "Ccl11"         "Frmpd4"        "Grid2"         "Csmd1"         "P2ry14"        "Mt1"           "Mbd1"          "Glis3"        
[10] "Timp3"         "Eya2"          "Aff1"          "9530026P05Rik" "Sugct"         "Spock3"        "Ntn4"          "Myo1d"         "9030624G23Rik"
[19] "Anxa1"         "Thsd7a"        "Smim3"         "Dcn"           "Ell2"          "Cfh"          

map(genes_hits_aging_up, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

genes_hits_aging_down <-
  moderated_lfc_leg_muscle %>%
  filter(sc_aging < -1 & svalue_sc_aging < 10^-8 |
    (sc_aging < -.5 & svalue_sc_aging < 10^-8 & bulk_aging < -2)) %>%
  arrange(sc_aging) %>%
  pull(gene_id)
genes_hits_aging_down

  [1] "Dgkg"                  "Sorbs2"                "Sparc"                 "Ccnd1"                 "Col3a1"                "6030407O03Rik"        
  [7] "Itm2a"                 "Kif21a"                "mt-Co2"                "Col4a1"                "mt-Nd3"                "Lgals1"               
 [13] "Gnas"                  "Rps18"                 "mt-Atp6"               "Meg3"                  "mt-Co3"                "Rps14"                
 [19] "Tmsb10"                "Rps23"                 "Apoe"                  "Hsp90ab1"              "Dag1"                  "Rpl13a"               
 [25] "Rpl5"                  "Rps6"                  "Gas1"                  "Gm28438"               "Serpinh1"              "Rps9"                 
 [31] "Gsn"                   "Rpl10a"                "Rpl7a-ps12"            "Uba52"                 "Pdlim4"                "Marcks"               
 [37] "Rpl4"                  "Rps25"                 "Ppp1r14b"              "Tmsb4x"                "Rps17"                 "Tagln"                
 [43] "Airn"                  "Crip2"                 "Rpl10"                 "Rack1"                 "Rpl18a"                "Cfap77"               
 [49] "Calr"                  "Atp5b"                 "Crip1"                 "Rpl34"                 "Msc"                   "Rps11"                
 [55] "Rpl27"                 "Rps15"                 "Rpl17"                 "Rpl26"                 "mt-Nd1"                "Rplp2"                
 [61] "Rplp0"                 "Mest"                  "Rpl41"                 "H3f3a"                 "Zbtb20"                "Rpl29"                
 [67] "Vim"                   "Galnt2l"               "Cox8a"                 "Cox4i1"                "Rpl11"                 "mt-Nd4"               
 [73] "mt-Co1"                "Igfbp5"                "Ubc"                   "Anxa6"                 "Rpl12"                 "Rpl37a"               
 [79] "S100a10"               "Hspa8"                 "Rps27a"                "Rpl10-ps3"             "Actb"                  "Ppia"                 
 [85] "Eef1b2"                "Rps29"                 "Cd63"                  "Col4a2"                "Gas5"                  "Uqcrh"                
 [91] "Cnn3"                  "Tpt1"                  "mt-Nd2"                "Rps27"                 "Anxa2"                 "Atp5l"                
 [97] "Ptma"                  "Kalrn"                 "Ptms"                  "Rpl23a"                "Oaz1"                  "Cp"                   
[103] "Gapdh"                 "Rpl24"                 "Mbnl2"                 "Igf1"                  "Gm15500"               "Nfib"                 
[109] "Rpl7a"                 "Rpl32"                 "Serf2"                 "Zeb2"                  "Rps3"                  "Rpl31"                
[115] "Rps7"                  "Rhoj"                  "Col6a1"                "Rps13"                 "Filip1l"               "Rpl35"                
[121] "Gxylt2"                "Myl6"                  "Rpl15"                 "Rps12"                 "Zfp36l1"               "Hmgb1"                
[127] "Rpl14"                 "Eef1g"                 "mt-Cytb"               "Cd81"                  "Eef2"                  "Atp5e"                
[133] "Cox6c"                 "Eef1a1"                "Rpl28"                 "Rpl36"                 "Crlf1"                 "mt-Rnr2"              
[139] "Anxa5"                 "Itm2b"                 "Rpl27a"                "Sh3glb1"               "Rps15a"                "Cox7c"                
[145] "Hspa5"                 "Kxd1"                  "Chodl"                 "Rpl3"                  "Grb10"                 "Cd82"                 
[151] "Cavin2"                "Nfia"                  "Peg3"                  "Lama2"                 "Rpl9"                  "ENSMUSG00000039607.16"
[157] "Rpl35a"                "Gm11808"               "Chchd2"                "Rpl19"                 "Naca"                  "Gm28661"              
[163] "Capns1"                "Pfn1"                  "Rplp1"                 "Atp1a2"                "Rps3a1"                "Cav1"                 
[169] "Rpl22"                 "Hsp90b1"               "Egr1"                  "Mmp2"                  "Gm10275"               "Lamc1"                
[175] "Nfix"                  "Nfic"                  "Hnrnpa0"              

map(genes_hits_aging_down[170:177], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

genes_hits_aging <- c(genes_hits_aging_down, genes_hits_aging_up)

cor_pairs_aging <-
  cor_pairs %>%
  filter((gene1 %in% genes_hits_aging & gene2 %in% genes_hits_aging) &
    abs(rho) >= .3) %>%
  arrange(gene1, rho)
cor_pairs_aging

cor_pairs_aging_genes <- sort(union(cor_pairs_aging$gene1, cor_pairs_aging$gene2))
cor_pairs_aging_genes

 [1] "Eef1a1"     "Eef1b2"     "Eef2"       "Gm11808"    "Gm28438"    "Kalrn"      "Kxd1"       "mt-Atp6"    "mt-Co1"     "mt-Co2"     "mt-Co3"    
[12] "mt-Cytb"    "mt-Nd1"     "mt-Nd2"     "mt-Nd3"     "mt-Nd4"     "Naca"       "Rack1"      "Rpl10"      "Rpl10-ps3"  "Rpl10a"     "Rpl11"     
[23] "Rpl12"      "Rpl13a"     "Rpl14"      "Rpl15"      "Rpl17"      "Rpl18a"     "Rpl19"      "Rpl22"      "Rpl24"      "Rpl26"      "Rpl27a"    
[34] "Rpl28"      "Rpl29"      "Rpl3"       "Rpl32"      "Rpl34"      "Rpl35"      "Rpl35a"     "Rpl36"      "Rpl37a"     "Rpl4"       "Rpl41"     
[45] "Rpl5"       "Rpl7a-ps12" "Rpl9"       "Rplp0"      "Rplp1"      "Rplp2"      "Rps11"      "Rps12"      "Rps13"      "Rps14"      "Rps15a"    
[56] "Rps18"      "Rps23"      "Rps27a"     "Rps29"      "Rps3"       "Rps3a1"     "Rps6"       "Rps7"       "Rps9"       "Tpt1"       "Uba52"     

setdiff(genes_hits_aging, cor_pairs_aging_genes)

  [1] "Dgkg"                  "Sorbs2"                "Sparc"                 "Ccnd1"                 "Col3a1"                "6030407O03Rik"        
  [7] "Itm2a"                 "Kif21a"                "Col4a1"                "Lgals1"                "Gnas"                  "Meg3"                 
 [13] "Tmsb10"                "Apoe"                  "Hsp90ab1"              "Dag1"                  "Gas1"                  "Serpinh1"             
 [19] "Gsn"                   "Pdlim4"                "Marcks"                "Rps25"                 "Ppp1r14b"              "Tmsb4x"               
 [25] "Rps17"                 "Tagln"                 "Airn"                  "Crip2"                 "Cfap77"                "Calr"                 
 [31] "Atp5b"                 "Crip1"                 "Msc"                   "Rpl27"                 "Rps15"                 "Mest"                 
 [37] "H3f3a"                 "Zbtb20"                "Vim"                   "Galnt2l"               "Cox8a"                 "Cox4i1"               
 [43] "Igfbp5"                "Ubc"                   "Anxa6"                 "S100a10"               "Hspa8"                 "Actb"                 
 [49] "Ppia"                  "Cd63"                  "Col4a2"                "Gas5"                  "Uqcrh"                 "Cnn3"                 
 [55] "Rps27"                 "Anxa2"                 "Atp5l"                 "Ptma"                  "Ptms"                  "Rpl23a"               
 [61] "Oaz1"                  "Cp"                    "Gapdh"                 "Mbnl2"                 "Igf1"                  "Gm15500"              
 [67] "Nfib"                  "Rpl7a"                 "Serf2"                 "Zeb2"                  "Rpl31"                 "Rhoj"                 
 [73] "Col6a1"                "Filip1l"               "Gxylt2"                "Myl6"                  "Zfp36l1"               "Hmgb1"                
 [79] "Eef1g"                 "Cd81"                  "Atp5e"                 "Cox6c"                 "Crlf1"                 "mt-Rnr2"              
 [85] "Anxa5"                 "Itm2b"                 "Sh3glb1"               "Cox7c"                 "Hspa5"                 "Chodl"                
 [91] "Grb10"                 "Cd82"                  "Cavin2"                "Nfia"                  "Peg3"                  "Lama2"                
 [97] "ENSMUSG00000039607.16" "Chchd2"                "Gm28661"               "Capns1"                "Pfn1"                  "Atp1a2"               
[103] "Cav1"                  "Hsp90b1"               "Egr1"                  "Mmp2"                  "Gm10275"               "Lamc1"                
[109] "Nfix"                  "Nfic"                  "Hnrnpa0"               "Gm7536"                "Ccl11"                 "Frmpd4"               
[115] "Grid2"                 "Csmd1"                 "P2ry14"                "Mt1"                   "Mbd1"                  "Glis3"                
[121] "Timp3"                 "Eya2"                  "Aff1"                  "9530026P05Rik"         "Sugct"                 "Spock3"               
[127] "Ntn4"                  "Myo1d"                 "9030624G23Rik"         "Anxa1"                 "Thsd7a"                "Smim3"                
[133] "Dcn"                   "Ell2"                  "Cfh"                  

plotExpression(sce10x_stem,
  features = setdiff(cor_pairs_aging_genes, "Sparc"),
  x = c("Sparc"),
  exprs_values = "logcounts",
  colour_by = "condition",
  point_size = .7,
  point_alpha = 0.7,
  ncol = 4
) + coord_fixed()

Enrichment analysis

msigdb_info()

|------------------------------------------------------------------|
| Via: R package msigdbr                                    v7.1.1 |
|------------------------------------------------------------------|
| Available Species                                                |
|------------------------------------------------------------------|
| Homo sapiens                                                     |
| Mus musculus                                                     |
| Drosophila melanogaster                                          |
| Gallus gallus                                                    |
| Saccharomyces cerevisiae                                         |
| Bos taurus                                                       |
| Caenorhabditis elegans                                           |
| Canis lupus familiaris                                           |
| Danio rerio                                                      |
| Rattus norvegicus                                                |
| Sus scrofa                                                       |
|------------------------------------------------------------------|
| Available Genesets                                               |
|------------------------------------------------------------------|
| Category Subcategory | Description                               |
|------------------------------------------------------------------|
| C1                   | Positional                                |
| C2 CGP               | Chemical and Genetic Perturbations        |
| C2 CP                | Canonical Pathways                        |
| C2 CP:BIOCARTA       | Canonical BIOCARTA                        |
| C2 CP:KEGG           | Canonical KEGG                            |
| C2 CP:PID            | Canonical PID                             |
| C2 CP:REACTOME       | Canonical REACTOME                        |
| C3 MIR               | Motif miRNA Targets                       |
| C3 TFT               | Motif Transcription Factor Targets        |
| C4 CGN               | Cancer Gene Neighborhoods                 |
| C4 CM                | Cancer Modules                            |
| C5 BP                | GO Biological Process                     |
| C5 CC                | GO Cellular Component                     |
| C5 MF                | GO Molecular Function                     |
| C6                   | Oncogenic Signatures                      |
| C7                   | Immunologic Signatures                    |
| H                    | Hallmark                                  |
|------------------------------------------------------------------|
| Source: http://software.broadinstitute.org/gsea/msigdb           |
|------------------------------------------------------------------|

HALLMARK <- msigdb_gsets(species = "Mus musculus", "H")

weighted_signatures <- vector("list", length = 6)
weighted_signatures[[1]] <- genes_hits_aging_down
weighted_signatures[[2]] <- genes_hits_aging_up
weighted_signatures[[3]] <- genes_hits_clust2vs1_down
weighted_signatures[[4]] <- genes_hits_clust2vs1_up
weighted_signatures[[5]] <- genes_hits_clust3vs12y_down
weighted_signatures[[6]] <- genes_hits_clust3vs12y_up
names(weighted_signatures) <- c("aging_down", "aging_up ", "clust2vs1_down", "clust2vs1_up", "clust3vs12y_down", "clust3vs12y_up")

hyp_obj <- hypeR(weighted_signatures, HALLMARK, test = "hypergeometric", fdr = .9, plotting = TRUE)

aging_down
aging_up 
clust2vs1_down
clust2vs1_up
clust3vs12y_down
clust3vs12y_up

hyp_obj$data

$aging_down
(hyp) 

  data: 

  plots: 30 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet
$`aging_up `
(hyp) 

  data: 

  plots: 14 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet
$clust2vs1_down
(hyp) 

  data: 

  plots: 15 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet
$clust2vs1_up
(hyp) 

  data: 

  plots: 10 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet
$clust3vs12y_down
(hyp) 

  data: 

  plots: 17 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet
$clust3vs12y_up
(hyp) 

  data: 

  plots: 13 Figures

  args: signature
        genesets
        test
        background
        power
        absolute
        pval
        fdr
        plotting
        quiet

hyp_dots(hyp_obj, val = "fdr")

$aging_down

$`aging_up `

$clust2vs1_down

$clust2vs1_up

$clust3vs12y_down

$clust3vs12y_up

hyp_to_excel(hyp_obj,
  file_path = file.path(data_dir, "preprocessed", "hallmark_enrichment_stem.xlsx")
)

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] EnhancedVolcano_1.6.0       ggrepel_0.8.2               hypeR_1.4.0                 pheatmap_1.0.12             forcats_0.5.0              
 [6] stringr_1.4.0               dplyr_1.0.2                 purrr_0.3.4                 readr_1.3.1                 tidyr_1.1.2                
[11] tibble_3.0.3                tidyverse_1.3.0             scran_1.16.0                scater_1.16.2               ggplot2_3.3.2              
[16] SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1         matrixStats_0.56.0          Biobase_2.48.0             
[21] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2         IRanges_2.22.2              S4Vectors_0.26.1            BiocGenerics_0.34.0        

loaded via a namespace (and not attached):
  [1] readxl_1.3.1              backports_1.1.9           igraph_1.2.5              splines_4.0.2             BiocParallel_1.22.0      
  [6] inline_0.3.16             digest_0.6.25             htmltools_0.5.0           viridis_0.5.1             fansi_0.4.1              
 [11] magrittr_1.5              openxlsx_4.1.5            limma_3.44.3              modelr_0.1.8              RcppParallel_5.0.2       
 [16] R.utils_2.10.1            prettyunits_1.1.1         colorspace_1.4-1          blob_1.2.1                rvest_0.3.6              
 [21] haven_2.3.1               xfun_0.17                 callr_3.4.4               crayon_1.3.4              RCurl_1.98-1.2           
 [26] jsonlite_1.7.1            glue_1.4.2                kableExtra_1.2.1          polyclip_1.10-0           gtable_0.3.0             
 [31] zlibbioc_1.34.0           XVector_0.28.0            webshot_0.5.2             V8_3.2.0                  R.cache_0.14.0           
 [36] pkgbuild_1.1.0            BiocSingular_1.4.0        rstan_2.21.2              scales_1.1.1              msigdbr_7.1.1            
 [41] DBI_1.1.0                 edgeR_3.30.3              Rcpp_1.0.5                viridisLite_0.3.0         dqrng_0.2.1              
 [46] rsvd_1.0.3                StanHeaders_2.21.0-6      htmlwidgets_1.5.1         httr_1.4.2                RColorBrewer_1.1-2       
 [51] ellipsis_0.3.1            loo_2.3.1                 pkgconfig_2.0.3           R.methodsS3_1.8.1         farver_2.0.3             
 [56] dbplyr_1.4.4              locfit_1.5-9.4            labeling_0.3              tidyselect_1.1.0          rlang_0.4.7              
 [61] munsell_0.5.0             cellranger_1.1.0          tools_4.0.2               visNetwork_2.0.9          cli_2.0.2                
 [66] generics_0.0.2            broom_0.7.0               evaluate_0.14             yaml_2.2.1                processx_3.4.4           
 [71] knitr_1.29                fs_1.5.0                  zip_2.1.1                 packrat_0.5.0             nlme_3.1-149             
 [76] reactable_0.2.1           R.oo_1.24.0               xml2_1.3.2                compiler_4.0.2            rstudioapi_0.11          
 [81] beeswarm_0.2.3            curl_4.3                  reprex_0.3.0              statmod_1.4.34            tweenr_1.0.1             
 [86] stringi_1.5.3             ps_1.3.4                  lattice_0.20-41           Matrix_1.2-18             styler_1.3.2             
 [91] vctrs_0.3.4               pillar_1.4.6              lifecycle_0.2.0           BiocNeighbors_1.6.0       cowplot_1.1.0            
 [96] bitops_1.0-6              irlba_2.3.3               R6_2.4.1                  gridExtra_2.3             vipor_0.4.5              
[101] codetools_0.2-16          MASS_7.3-52               assertthat_0.2.1          rprojroot_1.3-2           withr_2.2.0              
[106] GenomeInfoDbData_1.2.3    mgcv_1.8-33               hms_0.5.3                 grid_4.0.2                rmarkdown_2.3            
[111] DelayedMatrixStats_1.10.1 ggforce_0.3.2             lubridate_1.7.9           base64enc_0.1-3           ggbeeswarm_0.6.0