Mouse Muscle Stem Cell Project

Part 5b: run differential expression analysis on fap cells

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

Prepare analysis workflow

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(DESeq2)
  library(ashr)
  library(EnhancedVolcano)
  library(BiocParallel)
  library(hypeR)
  library(zinbwave)
  theme_set(theme_bw())
})

Define file paths

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

get_counts_stats <- function(sce10x_fap, cluster_to_remove) {
  gene_counts_by_cluster <-
    t(assay(sce10x_fap)) %>%
    as.matrix() %>%
    as_tibble() %>%
    group_by(cluster = colData(sce10x_fap)$clusters_condition) %>%
    summarise_all(list(~ mean(.)))

  umi_counts_max<-
    gene_counts_by_cluster[, -1] %>%
    summarize_all(list(~ max(.))) %>%
    mutate(dummy = "dummy") %>%
    pivot_longer(-dummy, names_to = "gene_name", values_to = "max_expr") %>%
    select(-dummy)

  stats <-
    perFeatureQCMetrics(sce10x_fap,
      exprs_values = c("counts"),
      flatten = TRUE
    ) %>%
    as_tibble() %>%
    bind_cols(., rowData(sce10x_fap) %>%
      as_tibble()) %>%
    bind_cols(., umi_counts_max)

  return(stats)
}

get_lfc <- function(contrast, name, dds) {
  lfcShrink(dds,
    contrast = contrast,
    type = "ashr",
    svalue = T
  ) %>%
    as_tibble() %>%
    dplyr::select(-baseMean) %>%
    mutate(svalue = abs(svalue)) %>%
    set_names(paste(c("lfc", "lfc_se", "svalue"), name, sep = "_"))
}

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

plot_volcano <- function(var_name, lfc_thresh, svalue_thresh, lfc, suffix, label=NULL) {
 p <- EnhancedVolcano(lfc,
    lab =lfc %>% pull(gene_id),
    x = paste0("lfc_", var_name),
    selectLab = label,
    y = paste0("svalue_", 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)
    )
  )
  
  ggsave(file.path(figures_dir, paste0("volcano_sc", var_name,"_effect_", suffix, ".pdf")), p)
  return(p)
}

Load data

sce10x <-
  readRDS(file.path(
    data_dir,
    "preprocessed",
    "sce10x_filtered_final.rds"
  ))

sce10x_fap <- sce10x[, colData(sce10x)$celltype == "fap"]
rm(sce10x)

table(colData(sce10x_fap)$clusters, colData(sce10x_fap)$sample)

       
        yng1 yng2 yng3 aged1 aged2 aged3 aged4
  fap_1  295  284  389  1630  1585  1134   833
  fap_2  203  106  248   390   424   434   194

Discard low counts genes

n_exprs_genes <-
  nexprs(sce10x_fap,
    detection_limit = 5,
    byrow = TRUE
  )
keep <- n_exprs_genes >= 10
table(keep)

keep
FALSE  TRUE 
29144  4330 

sce10x_fap <- sce10x_fap[keep, ]

batch_genes_remove <- c("Chad", "Cxcl13", "Spp1", "Cilp2", "Kcnq5", "Tnc", "Krtdap")

map(batch_genes_remove, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

sce10x_fap <- sce10x_fap[!(rownames(sce10x_fap) %in% batch_genes_remove), ]

Discard low cluster mean expression genes

stats <- get_counts_stats(sce10x_fap)
stats

max_thresh <- .3
stats %>%
  count(max_expr > max_thresh)

ggplot(stats) +
  geom_histogram(aes(max_expr), bins = 100) +
  scale_x_log10() +
  facet_wrap(~top_hvg_fap) +
  geom_vline(xintercept = max_thresh)

genes_to_plot <-
  stats %>%
  filter(max_expr <= max_thresh) %>%
  arrange(-max_expr)%>%
  pull(gene_name) 

map(genes_to_plot[1:5], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

genes_keep <-
  stats %>%
  filter(max_expr > max_thresh) %>%
  pull(gene_name)
length(genes_keep)

[1] 4133

sce10x_fap <- sce10x_fap[genes_keep, ]

ggplot(
  stats,
  aes(
    mean,
    detected
  )
) +
  scale_x_log10() +
  geom_point(size = 0.3, aes(color = max_expr <1)) +
  geom_text(aes(
    label = gene_name
  ),
  check_overlap = TRUE, nudge_y = -0.1, size = 2.5
  )

Create Design Matrix

design_fap <-
  model.matrix(~ -1 + clusters_condition + sample,
    data = colData(sce10x_fap)
  )[, -c(7)]
colnames(design_fap) <- str_replace(colnames(design_fap), "clusters_condition|sample", "")
# colnames(design_fap) <- str_replace(colnames(design_fap), ":condition", "")

design_fap[1:3, ]

                    fap_1_aged fap_1_yng fap_2_aged fap_2_yng yng2 yng3 aged2 aged3 aged4
GACTATGTCCGGCTTT_d1          1         0          0         0    0    0     0     0     0
GATTGGTAGGGAGTGG_d1          0         0          1         0    0    0     0     0     0
AAGACAATCTTCCAGC_d1          0         0          1         0    0    0     0     0     0

Compute Observational Weights

assay(sce10x_fap, "counts") <- round(assay(sce10x_fap, "counts"))

system.time({
  zinb <-
    zinbFit(sce10x_fap,
      K = 0,
      X = design_fap,
      verbose = TRUE,
      BPPARAM = MulticoreParam(3),
      epsilon = 1e12
    )
})

Create model:
ok
Initialize parameters:
ok
Optimize parameters:
Iteration 1
penalized log-likelihood = -59958111.6221251
After dispersion optimization = -107505848.923134

   user  system elapsed 
    244       3     125 

After right optimization = -107495622.381593
After orthogonalization = -107495622.381593

   user  system elapsed 
  208.3     3.5   114.5 

After left optimization = -99050659.2138885
After orthogonalization = -99050659.2138885
Iteration 2
penalized log-likelihood = -99050659.2138885
After dispersion optimization = -99050659.1494286

   user  system elapsed 
  229.7     2.5   121.4 

After right optimization = -99050483.8513361
After orthogonalization = -99050483.8513361

   user  system elapsed 
   44.0     1.8    23.4 

After left optimization = -99050483.4421824
After orthogonalization = -99050483.4421824
Iteration 3
penalized log-likelihood = -99050483.4421824
ok

   user  system elapsed 
   1821     105     948 

weights <- computeObservationalWeights(zinb, as.matrix(assay(sce10x_fap)))
dimnames(weights) <- dimnames(sce10x_fap)
assay(sce10x_fap, "weights") <- weights

convert to SCE to DESeqDataSet object

dds_fap <-
  convertTo(sce10x_fap, type = c("DESeq2"))

converting counts to integer mode

design(dds_fap) <- design_fap
assay(dds_fap, "weights") <- assay(sce10x_fap, "weights")

dds_fap <- estimateSizeFactors(dds_fap, type = "poscounts")
dds_fap <-
  DESeq(dds_fap,
    test = "LRT",
    useT = TRUE,
    reduced = design_fap[, 1:4],
    minmu = 1e-6,
    parallel = TRUE,
    BPPARAM = MulticoreParam(3),
    minRep = Inf
  )

using supplied model matrix
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates: 3 workers
mean-dispersion relationship
final dispersion estimates, fitting model and testing: 3 workers

plotDispEsts(dds_fap)

resultsNames(dds_fap)

[1] "fap_1_aged" "fap_1_yng"  "fap_2_aged" "fap_2_yng"  "yng2"       "yng3"       "aged2"      "aged3"      "aged4"     

Plot batch effects

batch_dt <-
  rowData(dds_fap) %>%
  as_tibble(rownames = "gene_id") %>%
  select("gene_id", "yng2":"aged4")

ggplot(
  batch_dt,
  aes(
    x = yng2,
    y = yng3
  )
) +
  geom_point(
    size = .2,
    alpha = 0.3
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  )

ggplot(
  batch_dt,
  aes(
    x = aged2,
    y = aged3
  )
) +
  geom_point(
    size = .2,
    alpha = 0.3
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  )

ggplot(
  batch_dt,
  aes(
    x = aged2,
    y = aged4
  )
) +
  geom_point(
    size = .2,
    alpha = 0.3
  ) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  )

Test contrasts

c1 <- c(1, -1, 1, -1, 0, 0, 0, 0, 0) / 2
c2 <- c(-1, -1, 1, 1, 0, 0, 0, 0, 0) / 2
c3 <- c(1, -1, 0, 0, 0, 0, 0, 0, 0)
c4 <- c(0, 0, 1, -1, 0, 0, 0, 0, 0)
c5 <- c(-1, 0, 1, 0, 0, 0, 0, 0, 0)
c6 <- c(0, -1, 0, 1, 0, 0, 0, 0, 0)


fap_contrast_vec <- list(
  fap_aging = c1,
  fap2_1 = c2,
  fap1_aging = c3,
  fap2_aging = c4,
  fap2_1a = c5,
  fap2_1y = c6
)

lfc_fap <-
  imap_dfc(fap_contrast_vec,
    get_lfc,
    dds = dds_fap
  ) %>%
  bind_cols(
    gene_id = rownames(dds_fap),
    .
  )

using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041

lfc_fap <-
  left_join(lfc_fap,
    rowData(dds_fap) %>%
      as_tibble(rownames = "gene_id"),
    by = "gene_id"
  )

lfc_fap

Save data

write_tsv(
  lfc_fap,
  file.path(
    data_dir,
    "preprocessed",
    "fap_scrnaseq_moderated_lfc.txt"
  )
)

rowData(sce10x_fap) <-
  left_join(rowData(sce10x_fap) %>%
    as_tibble(rownames = "gene_id") %>%
    select(gene_id), lfc_fap, by = "gene_id") %>%
  DataFrame(.)

rownames(rowData(sce10x_fap)) <- rowData(sce10x_fap)$gene_id

saveRDS(
  sce10x_fap,
  file.path(
    data_dir,
    "preprocessed",
    "sce10x_fap_filtered_final.rds"
  )
)

Make volcano plots

lfc_thresh <- 1
svalue_thresh <- 10e-8
volcano_plots <-
  map(names(fap_contrast_vec)[1:2],
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = lfc_fap,
  suffix="fap_labeled"
  )

Saving 7 x 7 in image

volcano_plots

[[1]]

[[2]]

lfc_thresh <- 1
svalue_thresh <- 10e-8
volcano_plots_unlabeled <-
  map(names(fap_contrast_vec)[1:2],
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = lfc_fap,
  suffix="fap_unlabeled",
  label=c("")
  )

Saving 7 x 7 in image

volcano_plots_unlabeled

[[1]]

[[2]]

genes_hits_aging_up <-
  lfc_fap %>%
  filter(lfc_fap_aging > 1) %>%
  arrange(-lfc_fap_aging) %>%
  pull(gene_id)
genes_hits_aging_up

 [1] "Apod"          "Csmd1"         "Fetub"         "Trf"           "Sbno2"         "Apoe"          "9030624G23Rik" "Cebpb"         "Runx1"        
[10] "C4b"           "Arih1"         "H2-Q4"         "Dbp"           "Cstf3"         "Rtn1"          "Gm13307"       "Galnt15"       "Nr1d2"        
[19] "Prg4"          "0610043K17Rik" "Gpx3"          "Spats2"        "Col12a1"       "Dram1"         "Dpyd"          "Marchf3"       "Sugct"        
[28] "Irak3"         "Il1r1"         "Gm20400"       "Cxcl1"         "Gm49767"       "Nop58"         "Aff1"          "Ell2"          "Abca6"        

map(genes_hits_aging_up[1:10], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

genes_hits_aging_down <-
  lfc_fap %>%
  filter(lfc_fap_aging < -1) %>%
  arrange(lfc_fap_aging) %>%
  pull(gene_id)
genes_hits_aging_down

  [1] "Col3a1"                "Col1a1"                "Cpz"                   "Sparc"                 "Col1a2"                "H19"                  
  [7] "C430049B03Rik"         "Nrk"                   "Itm2a"                 "Col6a1"                "Col6a2"                "Col5a1"               
 [13] "Col6a3"                "Col5a2"                "Unc13c"                "Adamts16"              "Col5a3"                "Col4a1"               
 [19] "Adamts17"              "Serpinh1"              "Sfrp2"                 "Mfap5"                 "mt-Nd3"                "Plac9a"               
 [25] "Tppp3"                 "Meg3"                  "Pcolce"                "mt-Co2"                "C1qtnf3"               "S100a10"              
 [31] "Lpl"                   "Nrep"                  "Hspg2"                 "Serpinf1"              "Nid1"                  "mt-Atp6"              
 [37] "Fbn1"                  "C1qtnf6"               "Marcks"                "Gm49450"               "Fstl1"                 "Dpt"                  
 [43] "Cd248"                 "B830012L14Rik"         "Ppic"                  "Sbsn"                  "Fam155a"               "Col15a1"              
 [49] "Col4a2"                "mt-Co3"                "Gm28438"               "S100a11"               "Maged2"                "Stmn4"                
 [55] "Adamts19"              "Fn1"                   "Mmp2"                  "Gm37899"               "Tmsb10"                "Adamtsl2"             
 [61] "Adamts2"               "Tmsb4x"                "Stmn1"                 "Calr"                  "Thy1"                  "mt-Rnr2"              
 [67] "Clec3b"                "Gm9780"                "Serf2"                 "Sphkap"                "Lamc1"                 "Ret"                  
 [73] "Ace"                   "Nupr1"                 "Ptn"                   "Ccdc80"                "Anxa2"                 "Adam12"               
 [79] "Cd34"                  "Angptl1"               "Gap43"                 "mt-Rnr1"               "Igfbp6"                "Lgals1"               
 [85] "Cd81"                  "Bcat1"                 "Ednra"                 "Gas1"                  "Tuba1a"                "Cavin3"               
 [91] "Islr"                  "Crip1"                 "Mfap4"                 "Htra1"                 "Vim"                   "Rcn3"                 
 [97] "Myl6"                  "Hhip"                  "Diaph3"                "Oaf"                   "Bmp1"                  "Ifi27l2a"             
[103] "Ppib"                  "Spon2"                 "Esrrg"                 "Gnas"                  "Rps18"                 "Ucp2"                 
[109] "Nme2"                  "H3f3a"                 "Pmp22"                 "Itih5"                 "Loxl2"                 "Ifi27"                
[115] "Robo2"                 "ENSMUSG00000002900.16" "Ppp2r2b"               "Timp2"                 "Rps14"                 "Rian"                 
[121] "Rps6"                  "mt-Nd4"                "Anxa5"                 "Fam167a"               "Rps23"                 "Ptms"                 
[127] "Rpl13a"                "Rps17"                 "Gm45213"               "Tmem119"               "Hsp90b1"               "Nrxn3"                
[133] "Lamb2"                 "Adarb2"                "Rpl41"                 "Calu"                  "ENSMUSG00000001175.15" "Tmed3"                
[139] "Cdkn1c"                "Creb3l1"              

map(genes_hits_aging_down[1:10], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

genes_hits_fap2_up <-
  lfc_fap %>%
  filter(lfc_fap2_1 > 1) %>%
  arrange(-lfc_fap2_1) %>%
  pull(gene_id)
genes_hits_fap2_up

 [1] "Sema3e"                "Adamts16"              "Sema3c"                "Dact2"                 "Pi16"                  "Cmah"                 
 [7] "Car8"                  "Sbsn"                  "Fbn1"                  "Cd55"                  "Ugdh"                  "Opcml"                
[13] "Limch1"                "Anxa3"                 "Sv2c"                  "Gpr1"                  "Fam167a"               "Pdgfc"                
[19] "Ano3"                  "Gap43"                 "Efemp1"                "Rorb"                  "Efhd1"                 "H19"                  
[25] "Stmn4"                 "Fn1"                   "Duox1"                 "Creb5"                 "Efna5"                 "Gfpt2"                
[31] "Krt80"                 "Uap1"                  "Robo1"                 "Cd248"                 "Adarb2"                "Adgrd1"               
[37] "Pcolce2"               "Heg1"                  "Tmem100"               "Dpp4"                  "Ildr2"                 "Igfbp5"               
[43] "Ackr3"                 "Grm7"                  "Flnb"                  "Axl"                   "Mfap5"                 "Ntrk3"                
[49] "ENSMUSG00000071984.10" "Adamts5"               "Pcsk6"                 "Emilin2"               "Ptprj"                 "Ly6c1"                
[55] "Pla1a"                 "Lurap1l"               "Itgb7"                 "Procr"                 "Kank1"                 "Ppp1r14b"             
[61] "Sdk1"                  "Pde8a"                 "Smpd3"                 "Ppp2r2b"               "Il18"                  "Aldh1a3"              

map(genes_hits_fap2_up[1:10], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

genes_hits_fap1_up <-
  lfc_fap %>%
  filter(lfc_fap2_1 < -1) %>%
  arrange(lfc_fap2_1) %>%
  pull(gene_id)
genes_hits_fap1_up

 [1] "Mgp"                   "Cxcl14"                "Dlk1"                  "Kcnb2"                 "Apod"                  "Fmod"                 
 [7] "Ednra"                 "Smoc2"                 "Cst3"                  "Clu"                   "Mfap4"                 "Col15a1"              
[13] "Fetub"                 "Thbs4"                 "Nrk"                   "Lum"                   "Necab1"                "Fbln7"                
[19] "9530026P05Rik"         "Frmpd4"                "Sparcl1"               "Adamtsl2"              "Ptn"                   "Gm4804"               
[25] "Sphkap"                "Igfbp7"                "Cilp"                  "Mylk"                  "Cfh"                   "Rasgrp2"              
[31] "Comp"                  "Crlf1"                 "Trf"                   "Gdf10"                 "Cdh11"                 "Apoe"                 
[37] "Bgn"                   "Piezo2"                "Hsd11b1"               "G0s2"                  "Meox1"                 "Cxcl9"                
[43] "Steap4"                "Lsamp"                 "Lpl"                   "Srpx"                  "Pde4d"                 "Csmd1"                
[49] "Cpe"                   "Alpl"                  "Mdk"                   "Col4a1"                "Cpxm2"                 "Olfml3"               
[55] "Adamts17"              "Inpp4b"                "Sorl1"                 "Col4a2"                "Angptl7"               "Kcnma1"               
[61] "Pdgfrl"                "Spry1"                 "Actn1"                 "Gas1"                  "Lepr"                  "Nrxn2"                
[67] "Dgkb"                  "Crispld1"              "Nup210l"               "Lrrtm3"                "Adamts9"               "Tmem176b"             
[73] "Ccl2"                  "Cp"                    "Gfra2"                 "Fmo2"                  "Gm29216"               "ENSMUSG00000031842.14"
[79] "Lama2"                 "Cygb"                  "Slc1a3"                "Adam12"               

map(genes_hits_fap1_up[1:10], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

p1 <- ggplot(lfc_fap,
  aes(
    x = lfc_fap_aging,
    y = lfc_fap2_1
  )
) +
  geom_point(
    size = 0.5,
    alpha=0.5
  )  +
  geom_smooth(alpha=0.2, colour="grey")+
  coord_fixed(ratio = 1) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  theme_classic()
p1

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

Saving 12 x 7.41 in image

p2 <-ggplot(
  lfc_fap,
  aes(
    x = lfc_fap_aging,
    y = lfc_fap2_1
  )
) +
  geom_point(
    size = 0.5,
    alpha=0.5
  )  +
  geom_smooth(alpha=0.2, colour="grey")+
  coord_fixed(ratio = 1) +
  theme_classic()
p2

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

Saving 12 x 7.41 in image

ggplot(
  lfc_fap,
  aes(
    x = lfc_fap1_aging,
    y = lfc_fap2_aging
  )
) +
  geom_point(
    size = 0.5,
    alpha=0.5
  )  +
  geom_smooth(alpha=0.2, colour="grey")+
  coord_fixed(ratio = 1) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  theme_classic()

ggplot(
  lfc_fap,
  aes(
    x = lfc_fap2_1a,
    y = lfc_fap2_1y
  )
) +
  geom_point(
    size = 0.5,
    alpha=0.5
  )  +
  geom_smooth(alpha=0.2, colour="grey")+
  coord_fixed(ratio = 1) +
  geom_text(aes(
    label = gene_id
  ),
  check_overlap = TRUE, nudge_y = -0.15, size = 3
  ) +
  theme_classic()

map(c("Nrk", "Clu", "Sbsn"), plot_expression)

[[1]]

[[2]]

[[3]]

Enrichment analysis

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

weighted_signatures <- vector("list", length = 4)
weighted_signatures[[1]] <- genes_hits_aging_down
weighted_signatures[[2]] <- genes_hits_aging_up
weighted_signatures[[3]] <- genes_hits_fap2_up
weighted_signatures[[4]] <- genes_hits_fap1_up
names(weighted_signatures) <- c("aging_down", "aging_up ", "fap2_up", "fap1_up")

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

aging_down
aging_up 
fap2_up
fap1_up

hyp_obj$data

$aging_down
(hyp) 

  data: 

  plots: 26 Figures

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

  data: 

  plots: 23 Figures

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

  data: 

  plots: 28 Figures

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

  data: 

  plots: 18 Figures

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

hyp_dots(hyp_obj, val = "fdr")

$aging_down

$`aging_up `

$fap2_up

$fap1_up

hyp_to_excel(hyp_obj, 
             file_path=file.path(data_dir, "preprocessed", "hallmark_enrichment_fap.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] zinbwave_1.10.0             hypeR_1.4.0                 BiocParallel_1.22.0         EnhancedVolcano_1.6.0       ggrepel_0.8.2              
 [6] ashr_2.2-47                 DESeq2_1.28.1               forcats_0.5.0               stringr_1.4.0               dplyr_1.0.2                
[11] purrr_0.3.4                 readr_1.3.1                 tidyr_1.1.2                 tibble_3.0.3                tidyverse_1.3.0            
[16] scran_1.16.0                scater_1.16.2               ggplot2_3.3.2               SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2
[21] DelayedArray_0.14.1         matrixStats_0.56.0          Biobase_2.48.0              GenomicRanges_1.40.0        GenomeInfoDb_1.24.2        
[26] 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             digest_0.6.25            
  [6] invgamma_1.1              htmltools_0.5.0           viridis_0.5.1             SQUAREM_2020.4            fansi_0.4.1              
 [11] magrittr_1.5              memoise_1.1.0             openxlsx_4.1.5            limma_3.44.3              annotate_1.66.0          
 [16] modelr_0.1.8              colorspace_1.4-1          blob_1.2.1                rvest_0.3.6               haven_2.3.1              
 [21] xfun_0.17                 crayon_1.3.4              RCurl_1.98-1.2            jsonlite_1.7.1            genefilter_1.70.0        
 [26] survival_3.1-12           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             BiocSingular_1.4.0        scales_1.1.1             
 [36] msigdbr_7.1.1             DBI_1.1.0                 edgeR_3.30.3              Rcpp_1.0.5                viridisLite_0.3.0        
 [41] xtable_1.8-4              dqrng_0.2.1               bit_4.0.4                 rsvd_1.0.3                truncnorm_1.0-8          
 [46] htmlwidgets_1.5.1         httr_1.4.2                RColorBrewer_1.1-2        ellipsis_0.3.1            pkgconfig_2.0.3          
 [51] XML_3.99-0.5              farver_2.0.3              dbplyr_1.4.4              locfit_1.5-9.4            labeling_0.3             
 [56] tidyselect_1.1.0          rlang_0.4.7               softImpute_1.4            AnnotationDbi_1.50.3      munsell_0.5.0            
 [61] cellranger_1.1.0          tools_4.0.2               visNetwork_2.0.9          cli_2.0.2                 generics_0.0.2           
 [66] RSQLite_2.2.0             broom_0.7.0               evaluate_0.14             yaml_2.2.1                knitr_1.29               
 [71] bit64_4.0.5               fs_1.5.0                  zip_2.1.1                 packrat_0.5.0             nlme_3.1-149             
 [76] reactable_0.2.1           xml2_1.3.2                compiler_4.0.2            rstudioapi_0.11           beeswarm_0.2.3           
 [81] reprex_0.3.0              statmod_1.4.34            tweenr_1.0.1              geneplotter_1.66.0        stringi_1.5.3            
 [86] lattice_0.20-41           Matrix_1.2-18             vctrs_0.3.4               pillar_1.4.6              lifecycle_0.2.0          
 [91] BiocNeighbors_1.6.0       cowplot_1.1.0             bitops_1.0-6              irlba_2.3.3               R6_2.4.1                 
 [96] gridExtra_2.3             vipor_0.4.5               MASS_7.3-52               assertthat_0.2.1          rprojroot_1.3-2          
[101] withr_2.2.0               GenomeInfoDbData_1.2.3    mgcv_1.8-33               hms_0.5.3                 grid_4.0.2               
[106] rmarkdown_2.3             DelayedMatrixStats_1.10.1 mixsqp_0.3-43             ggforce_0.3.2             lubridate_1.7.9          
[111] base64enc_0.1-3           ggbeeswarm_0.6.0