Mouse Muscle Stem Cell Project

Part 5a: run differential expression analysis on stem 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(EnhancedVolcano)
  library(BiocParallel)
  library(zinbwave)
  theme_set(theme_bw())
})

Define file paths

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

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

  umi_counts_max_1 <-
    gene_counts_by_cluster[c(2, 4, 5), -1] %>%
    summarize_all(list(~ max(.))) %>%
    mutate(dummy = "dummy") %>%
    pivot_longer(-dummy, names_to = "gene_name", values_to = "max_yng") %>%
    select(-dummy)

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



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

  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_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) {
  EnhancedVolcano(lfc,
    lab = lfc %>% pull(gene_id),
    x = paste0("lfc_", var_name),
    xlim = c(-4, 4),
    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)
    )
  )
}

Load data

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

sce10x_stem <- sce10x[, colData(sce10x)$celltype == "stem"]
rm(sce10x)

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

        
         yng1 yng2 yng3 aged1 aged2 aged3 aged4
  stem_1  397  661  541   370   404   185   435
  stem_2  549 1030  379   151   278   141   304
  stem_3  253  222  129     0     0     0     0

Discard low counts genes

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

keep
FALSE  TRUE 
30184  3290 

sce10x_stem <- sce10x_stem[keep, ]

batch_genes_remove <- c("Gm12918", "Rps23-ps1", "Diaph3", "Top2a", "Ecrg4", "Rpl9-ps6")

map(batch_genes_remove, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

sce10x_stem <- sce10x_stem[!(rowData(sce10x_stem)$external_gene_name %in% batch_genes_remove), ]

Discard low cluster mean expression genes

stats <- get_counts_stats(sce10x_stem)
stats

max_thresh <- .3
stats <-
  stats %>%
  mutate(
    highexpr_yng = max_yng > max_thresh,
    highexpr_stem12 = max_stem12 > max_thresh,
    genes_stem3 = highexpr_yng & !highexpr_stem12,
    genes_stem12 = !highexpr_yng & highexpr_stem12,
    genes_neither = !highexpr_yng & !highexpr_stem12,
    genes_both = highexpr_yng & highexpr_stem12
  )

stats %>%
  count(highexpr_yng, highexpr_stem12)

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

genes_to_plot <-
  stats %>%
  filter(genes_stem3) %>%
  pull(gene_name) %>%
  sort()
genes_to_plot

 [1] "Adam12"  "Atp2a3"  "Bcl2l11" "Ccnb2"   "Cd24a"   "Cdh20"   "Cdkn1c"  "Cenpa"   "Cenpp"   "Ckap2"   "Cks2"    "Ctnna2"  "Dlk1"    "Dpp6"   
[15] "Efna5"   "Ehd4"    "Epb41"   "Fgd4"    "Gm28653" "Gm49397" "H1f10"   "H1f2"    "Hmox1"   "Igfbp3"  "Il1rap"  "Jpx"     "Khdrbs3" "L3mbtl4"
[29] "Lmnb1"   "Lmnb2"   "Lockd"   "Lrrn1"   "Mcm3"    "Mcm5"    "Mcm6"    "Mcm7"    "Mitf"    "Mthfd1l" "Myog"    "Nes"     "Nr2f2"   "Oca2"   
[43] "Ophn1"   "Pak3"    "Pclaf"   "Pola1"   "Prune2"  "Rrm2"    "Sdk1"    "Slc12a2" "Smc2"    "Stmn1"   "Whrn"    "Ypel2"  

map(genes_to_plot[1:5], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

genes_to_plot <-
  stats %>%
  filter(genes_stem12) %>%
  pull(gene_name) %>%
  sort()
genes_to_plot

 [1] "0610043K17Rik"         "AL954352.1"            "Atxn1"                 "Cadps2"                "Ccl11"                 "Ccn2"                 
 [7] "Cdh4"                  "Col8a1"                "Cped1"                 "Cyp26a1"               "Dab2"                  "ENSMUSG00000031842.14"
[13] "Eya2"                  "Frmpd4"                "Gm10288"               "Gm26834"               "Gm29216"               "Gm6225"               
[19] "Gm7536"                "Grid2"                 "Iigp1"                 "Macrod2"               "Matn4"                 "Mbd1"                 
[25] "Mkx"                   "Myo1d"                 "Nap1l5"                "Nr4a2"                 "P2ry14"                "Pak1"                 
[31] "Parp14"                "Pid1"                  "Rnf180"                "Rnf213"                "Samhd1"                "Sgip1"                
[37] "Slpi"                  "Spock3"                "Thsd7a"                "Tshr"                 

map(genes_to_plot[1:5], plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

genes_keep <-
  stats %>%
  filter(!genes_neither) %>%
  pull(gene_name)
length(genes_keep)

[1] 3196

sce10x_stem <- sce10x_stem[genes_keep, ]

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

Create Design Matrix

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

design_stem[1:3, ]

                    stem_1_aged stem_1_yng stem_2_aged stem_2_yng stem_3_yng yng2 yng3 aged2 aged3 aged4
CTACAGAGTGGCCCAT_d1           1          0           0          0          0    0    0     0     0     0
CTAACCCCACACCTGG_d1           1          0           0          0          0    0    0     0     0     0
AGGAGGTCACAGGATG_d1           1          0           0          0          0    0    0     0     0     0

Compute Observational Weights

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

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

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

   user  system elapsed 
  203.0     2.7   105.7 

After right optimization = -71574036.0998483
After orthogonalization = -71574036.0998483

   user  system elapsed 
  122.9     1.5    76.8 

After left optimization = -64813825.7211657
After orthogonalization = -64813825.7211657
Iteration 2
penalized log-likelihood = -64813825.7211657
After dispersion optimization = -64813825.7667293

   user  system elapsed 
  277.5     3.2   141.0 

After right optimization = -64813722.3999213
After orthogonalization = -64813722.3999213

   user  system elapsed 
   27.5     1.3    15.2 

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

   user  system elapsed 
   1412      62     669 

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

convert to SCE to DESeqDataSet object

dds_stem <-
  convertTo(sce10x_stem, type = c("DESeq2"))

converting counts to integer mode

design(dds_stem) <- design_stem
assay(dds_stem, "weights") <- assay(sce10x_stem, "weights")

dds_stem <- estimateSizeFactors(dds_stem, type = "poscounts")
dds_stem <-
  DESeq(dds_stem,
    test = "LRT",
    useT = TRUE,
    reduced = design_stem[, 1:5],
    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_stem)

resultsNames(dds_stem)

 [1] "stem_1_aged" "stem_1_yng"  "stem_2_aged" "stem_2_yng"  "stem_3_yng"  "yng2"        "yng3"        "aged2"       "aged3"       "aged4"      

Plot batch effects

batch_dt <-
  rowData(dds_stem) %>%
  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
  )

batch_genes <- c(
  "Cyp26a1", "Cenpa", "Dgkg", "Sorbs2",
  "Gm28653", "Pclaf", "Myog", "Prune2"
)
map(batch_genes, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

Test contrasts

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


aged_yng_contrast_vec <- list(
  stem_aging = c1,
  stem2_1 = c2,
  stem1_aging = c3,
  stem2_aging = c4,
  stem2_1a = c5,
  stem2_1y = c6
)

stem3_contrast_vec <- list(
  stem3_12y = c7,
  stem3_1y = c8,
  stem3_2y = c9
)

genes_stem3_contrasts <-
  stats %>%
  filter(genes_stem3 | genes_both) %>%
  pull(gene_name)

genes_aged_yng_contrasts <-
  stats %>%
  filter(genes_stem12 | genes_both) %>%
  pull(gene_name)

lfc_stem3 <-
  imap_dfc(stem3_contrast_vec,
    get_lfc,
    dds = dds_stem[genes_stem3_contrasts, ]
  ) %>%
  bind_cols(
    gene_id = rownames(dds_stem[genes_stem3_contrasts, ]),
    .
  )

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_aged_yng <-
  imap_dfc(aged_yng_contrast_vec,
    get_lfc,
    dds = dds_stem[genes_aged_yng_contrasts, ]
  ) %>%
  bind_cols(
    gene_id = rownames(dds_stem[genes_aged_yng_contrasts, ]),
    .
  )

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 <- full_join(lfc_aged_yng, lfc_stem3, by = "gene_id")

lfc <-
  left_join(lfc,
    rowData(dds_stem) %>%
      as_tibble(rownames = "gene_id"),
    by = "gene_id"
  )

lfc

Make volcano plots

lfc_thresh <- 1

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls

svalue_thresh <- 10e-8
volcano_plots <-
  map(names(c(aged_yng_contrast_vec, stem3_contrast_vec)),
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = lfc
  )

One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

volcano_plots

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

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

Saving 7 x 7 in image

ggsave(file.path(figures_dir, "volcano_stem_cluster2vs1.pdf"), volcano_plots[[2]])
ggsave(file.path(figures_dir, "volcano_stem_cluster3vs1_2.pdf"), volcano_plots[[7]])

Make scatterplots

aging_high_de <- c("Dgkg", "Sorbs2")

map(aging_high_de, plot_expression)

[[1]]

[[2]]

p1 <-
  ggplot(
    lfc %>% filter(!(gene_id %in% aging_high_de)),
    aes(
      x = lfc_stem_aging,
      y = lfc_stem2_1,
      color = lfc_stem3_12y
    )
  ) +
  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() +
  xlim(c(-6, 3))
p1

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

Saving 7 x 7 in image

p2 <-
  ggplot(
    lfc,
    aes(
      x = lfc_stem2_1,
      y = lfc_stem3_12y,
      color = lfc_stem_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"
  ) +
  ylim(c(-2.2, 1.8)) +
  theme_classic()
p2

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

Saving 7 x 7 in image

p3 <-
  ggplot(
    lfc %>% filter(!(gene_id %in% aging_high_de)),
    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()
p3

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

Saving 7 x 7 in image

Save data

write_tsv(
  lfc,
  file.path(
    data_dir,
    "preprocessed",
    "stem_scrnaseq_moderated_lfc.txt"
  )
)

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

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

saveRDS(
  sce10x_stem,
  file.path(
    data_dir,
    "preprocessed",
    "sce10x_stem_filtered_final.rds"
  )
)

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             BiocParallel_1.22.0         EnhancedVolcano_1.6.0       ggrepel_0.8.2               DESeq2_1.28.1              
 [6] forcats_0.5.0               stringr_1.4.0               dplyr_1.0.2                 purrr_0.3.4                 readr_1.3.1                
[11] tidyr_1.1.2                 tibble_3.0.3                tidyverse_1.3.0             scran_1.16.0                scater_1.16.2              
[16] ggplot2_3.3.2               SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1         matrixStats_0.56.0         
[21] Biobase_2.48.0              GenomicRanges_1.40.0        GenomeInfoDb_1.24.2         IRanges_2.22.2              S4Vectors_0.26.1           
[26] BiocGenerics_0.34.0        

loaded via a namespace (and not attached):
 [1] ggbeeswarm_0.6.0          colorspace_1.4-1          ellipsis_0.3.1            rprojroot_1.3-2           XVector_0.28.0           
 [6] BiocNeighbors_1.6.0       fs_1.5.0                  rstudioapi_0.11           farver_2.0.3              bit64_4.0.5              
[11] AnnotationDbi_1.50.3      fansi_0.4.1               lubridate_1.7.9           xml2_1.3.2                splines_4.0.2            
[16] geneplotter_1.66.0        knitr_1.29                jsonlite_1.7.1            packrat_0.5.0             broom_0.7.0              
[21] annotate_1.66.0           ashr_2.2-47               dbplyr_1.4.4              compiler_4.0.2            httr_1.4.2               
[26] dqrng_0.2.1               backports_1.1.9           assertthat_0.2.1          Matrix_1.2-18             limma_3.44.3             
[31] cli_2.0.2                 BiocSingular_1.4.0        tools_4.0.2               rsvd_1.0.3                igraph_1.2.5             
[36] gtable_0.3.0              glue_1.4.2                GenomeInfoDbData_1.2.3    Rcpp_1.0.5                softImpute_1.4           
[41] cellranger_1.1.0          vctrs_0.3.4               DelayedMatrixStats_1.10.1 xfun_0.17                 rvest_0.3.6              
[46] lifecycle_0.2.0           irlba_2.3.3               statmod_1.4.34            XML_3.99-0.5              edgeR_3.30.3             
[51] zlibbioc_1.34.0           scales_1.1.1              hms_0.5.3                 RColorBrewer_1.1-2        memoise_1.1.0            
[56] gridExtra_2.3             SQUAREM_2020.4            stringi_1.5.3             RSQLite_2.2.0             genefilter_1.70.0        
[61] truncnorm_1.0-8           rlang_0.4.7               pkgconfig_2.0.3           bitops_1.0-6              invgamma_1.1             
[66] lattice_0.20-41           labeling_0.3              cowplot_1.1.0             bit_4.0.4                 tidyselect_1.1.0         
[71] magrittr_1.5              R6_2.4.1                  generics_0.0.2            DBI_1.1.0                 pillar_1.4.6             
[76] haven_2.3.1               withr_2.2.0               mixsqp_0.3-43             survival_3.1-12           RCurl_1.98-1.2           
[81] modelr_0.1.8              crayon_1.3.4              viridis_0.5.1             locfit_1.5-9.4            grid_4.0.2               
[86] readxl_1.3.1              blob_1.2.1                reprex_0.3.0              digest_0.6.25             xtable_1.8-4             
[91] munsell_0.5.0             beeswarm_0.2.3            viridisLite_0.3.0         vipor_0.4.5