Mouse Muscle Stem Cell Project

Part 5c: run differential expression analysis on macrophage 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(zinbwave)
  theme_set(theme_bw())
})

Define file paths

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

get_counts_stats <- function(sce10x_macrophage, cluster_to_remove) {
  gene_counts_by_cluster <-
    t(assay(sce10x_macrophage)) %>%
    as.matrix() %>%
    as_tibble() %>%
    group_by(cluster = colData(sce10x_macrophage)$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_macrophage,
      exprs_values = c("counts"),
      flatten = TRUE
    ) %>%
    as_tibble() %>%
    bind_cols(., rowData(sce10x_macrophage) %>%
      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_macrophage,
    features = gene,
    x = "sample",
    exprs_values = "logcounts",
    colour_by = "condition",
    point_size = 1
  ) +
    facet_wrap(~ colData(sce10x_macrophage)$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_macrophage <- sce10x[, colData(sce10x)$celltype == "macrophage"]
rm(sce10x)

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

            
             yng1 yng2 yng3 aged1 aged2 aged3 aged4
  macrophage  282  274  398    58   245   149   189

Discard low counts genes

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

keep
FALSE  TRUE 
30428  3046 

sce10x_macrophage <- sce10x_macrophage[keep, ]

batch_genes_remove <-
  c(
    "Ccr7", "Cd209c", "Ifit1", "Slc40a1",
    "Fgf13", "Igfbp5", "Gpnmb", "Upp1", "Irgm2",
    "S100a8", "Chil3", "Saa3", "S100a9", "Apol7c", "F5", "Upp1"
  )

map(batch_genes_remove, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

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

Discard low cluster mean expression genes

stats <- get_counts_stats(sce10x_macrophage)
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_macrophage) +
  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] 2977

sce10x_macrophage <- sce10x_macrophage[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_macrophage <-
  model.matrix(~ -1 + condition + sample,
    data = colData(sce10x_macrophage)
  )[, -c(5)]
colnames(design_macrophage) <- str_replace(colnames(design_macrophage), "condition|sample", "")

design_macrophage[1:3, ]

                    yng aged yng2 yng3 aged2 aged3 aged4
ATGGTTGGTTGTAAAG_d1   0    1    0    0     0     0     0
GATGAGGGTTCAGTAC_d1   0    1    0    0     0     0     0
GCAACCGTCACCTTGC_d1   0    1    0    0     0     0     0

Compute Observational Weights

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

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

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

   user  system elapsed 
   45.7     2.2    24.6 

After right optimization = -23318260.1799767
After orthogonalization = -23318260.1799767

   user  system elapsed 
   30.8     1.3    17.4 

After left optimization = -20653340.6726067
After orthogonalization = -20653340.6726067
Iteration 2
penalized log-likelihood = -20653340.6726067
After dispersion optimization = -20653340.6757374

   user  system elapsed 
   32.3     1.1    33.8 

After right optimization = -20653321.6706135
After orthogonalization = -20653321.6706135

   user  system elapsed 
   6.78    0.74    3.76 

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

   user  system elapsed 
    333      22     158 

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

convert to SCE to DESeqDataSet object

dds_macrophage <-
  convertTo(sce10x_macrophage, type = c("DESeq2"))

converting counts to integer mode

design(dds_macrophage) <- design_macrophage
assay(dds_macrophage, "weights") <- assay(sce10x_macrophage, "weights")

dds_macrophage <- estimateSizeFactors(dds_macrophage, type = "poscounts")
dds_macrophage <-
  DESeq(dds_macrophage,
    test = "LRT",
    useT = TRUE,
    reduced = design_macrophage[, 1:2],
    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
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates, fitting model and testing: 3 workers

plotDispEsts(dds_macrophage)

resultsNames(dds_macrophage)

[1] "yng"   "aged"  "yng2"  "yng3"  "aged2" "aged3" "aged4"

Plot batch effects

batch_dt <-
  rowData(dds_macrophage) %>%
  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, 0, 0, 0, 0, 0)

macrophage_contrast_vec <- list(
  macrophage_aging = c1
)

lfc_macrophage <-
  imap_dfc(macrophage_contrast_vec,
    get_lfc,
    dds = dds_macrophage
  ) %>%
  bind_cols(
    gene_id = rownames(dds_macrophage),
    .
  )

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_macrophage <-
  left_join(lfc_macrophage,
    rowData(dds_macrophage) %>%
      as_tibble(rownames = "gene_id"),
    by = "gene_id"
  )

lfc_macrophage

Save data

write_tsv(
  lfc_macrophage,
  file.path(
    data_dir,
    "preprocessed",
    "macrophage_scrnaseq_moderated_lfc.txt"
  )
)

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

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

saveRDS(
  sce10x_macrophage,
  file.path(
    data_dir,
    "preprocessed",
    "sce10x_macrophage_filtered_final.rds"
  )
)

Make volcano plots

lfc_thresh <- 1
svalue_thresh <- 10e-8
volcano_plots <-
  map(names(macrophage_contrast_vec),
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = lfc_macrophage,
    suffix = "macrophage_labeled"
  )

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

volcano_plots

[[1]]

lfc_thresh <- 1
svalue_thresh <- 10e-8
volcano_plots_unlabeled <-
  map(names(macrophage_contrast_vec),
    plot_volcano,
    lfc_thresh = lfc_thresh,
    svalue_thresh = svalue_thresh,
    lfc = lfc_macrophage,
    suffix = "macrophage_unlabeled",
    label = c("")
  )

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

volcano_plots_unlabeled

[[1]]

genes_hits_aging_up <-
  lfc_macrophage %>%
  filter(lfc_macrophage_aging > 1) %>%
  arrange(-lfc_macrophage_aging) %>%
  pull(gene_id)
genes_hits_aging_up

 [1] "mt-Nd3"                "Zfp36"                 "Ifi27l2a"              "Rnase4"                "Folr2"                 "Apoe"                 
 [7] "F13a1"                 "mt-Co2"                "Cst3"                  "Cd83"                  "ENSMUSG00000034708.11" "Rps23"                

map(genes_hits_aging_up, plot_expression)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

genes_hits_aging_down <-
  lfc_macrophage %>%
  filter(lfc_macrophage_aging < -1) %>%
  arrange(lfc_macrophage_aging) %>%
  pull(gene_id)
genes_hits_aging_down

 [1] "Plcb1"                          "Lyz1"                           "Gsr"                            "Tmsb10"                        
 [5] "Vcan"                           "Msrb1"                          "Thbs1"                          "Prdx5"                         
 [9] "Smpdl3a"                        "Grk5"                           "Slc41a2"                        "Arhgap26_ENSMUSG00000036452.18"
[13] "Adgre4"                         "Rbm3"                           "Ern1"                           "Samhd1"                        
[17] "Cytip"                          "Emilin2"                        "Atp11b"                         "Sipa1l1"                       
[21] "Cers6"                          "Pbx1"                           "Abtb2"                          "Gda"                           
[25] "Jarid2"                         "Nrg1"                           "Runx2"                         

map(genes_hits_aging_down, 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]]

[[25]]

[[26]]

[[27]]

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               ashr_2.2-47                
 [6] DESeq2_1.28.1               forcats_0.5.0               stringr_1.4.0               dplyr_1.0.2                 purrr_0.3.4                
[11] readr_1.3.1                 tidyr_1.1.2                 tibble_3.0.3                tidyverse_1.3.0             scran_1.16.0               
[16] scater_1.16.2               ggplot2_3.3.2               SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1        
[21] matrixStats_0.56.0          Biobase_2.48.0              GenomicRanges_1.40.0        GenomeInfoDb_1.24.2         IRanges_2.22.2             
[26] S4Vectors_0.26.1            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] base64enc_0.1-3           BiocNeighbors_1.6.0       fs_1.5.0                  rstudioapi_0.11           farver_2.0.3             
 [11] bit64_4.0.5               AnnotationDbi_1.50.3      fansi_0.4.1               lubridate_1.7.9           xml2_1.3.2               
 [16] splines_4.0.2             R.methodsS3_1.8.1         geneplotter_1.66.0        knitr_1.29                jsonlite_1.7.1           
 [21] packrat_0.5.0             broom_0.7.0               annotate_1.66.0           dbplyr_1.4.4              R.oo_1.24.0              
 [26] compiler_4.0.2            httr_1.4.2                dqrng_0.2.1               backports_1.1.9           assertthat_0.2.1         
 [31] Matrix_1.2-18             limma_3.44.3              cli_2.0.2                 BiocSingular_1.4.0        htmltools_0.5.0          
 [36] tools_4.0.2               rsvd_1.0.3                igraph_1.2.5              gtable_0.3.0              glue_1.4.2               
 [41] GenomeInfoDbData_1.2.3    Rcpp_1.0.5                softImpute_1.4            cellranger_1.1.0          styler_1.3.2             
 [46] vctrs_0.3.4               DelayedMatrixStats_1.10.1 xfun_0.17                 rvest_0.3.6               lifecycle_0.2.0          
 [51] irlba_2.3.3               statmod_1.4.34            XML_3.99-0.5              edgeR_3.30.3              zlibbioc_1.34.0          
 [56] scales_1.1.1              hms_0.5.3                 RColorBrewer_1.1-2        rematch2_2.1.2            yaml_2.2.1               
 [61] memoise_1.1.0             gridExtra_2.3             SQUAREM_2020.4            stringi_1.5.3             RSQLite_2.2.0            
 [66] genefilter_1.70.0         truncnorm_1.0-8           rlang_0.4.7               pkgconfig_2.0.3           bitops_1.0-6             
 [71] invgamma_1.1              evaluate_0.14             lattice_0.20-41           labeling_0.3              cowplot_1.1.0            
 [76] bit_4.0.4                 tidyselect_1.1.0          magrittr_1.5              R6_2.4.1                  generics_0.0.2           
 [81] DBI_1.1.0                 pillar_1.4.6              haven_2.3.1               withr_2.2.0               mixsqp_0.3-43            
 [86] survival_3.1-12           RCurl_1.98-1.2            modelr_0.1.8              crayon_1.3.4              rmarkdown_2.3            
 [91] viridis_0.5.1             locfit_1.5-9.4            grid_4.0.2                readxl_1.3.1              blob_1.2.1               
 [96] reprex_0.3.0              digest_0.6.25             xtable_1.8-4              R.cache_0.14.0            R.utils_2.10.1           
[101] munsell_0.5.0             beeswarm_0.2.3            viridisLite_0.3.0         vipor_0.4.5