TL;DR

You may click on the following links to fast forward to the codes directly producing the Figures:

Back to index page.

About this Vignette

One of the key questions asked by the reviewers is whether the transduction level (as measured by mScarlet expression mScarlet) is associated with activation of microglia, i.e. relationship between mScarlet level and signature gene. To address this point, we use the Smart-seq2 dataset to make violin plots of signature genes across groups, with overlaying jitter points colored by mScarlet levels.

Initialization

suppressPackageStartupMessages({
  library(Seurat)
  library(SeuratDisk)
  library(tidyverse)
  library(ggplot2)
  library(ggpubr)
  library(patchwork)
  library(RColorBrewer)
})
querydata <- LoadH5Seurat('../data/AAV_MG_scRNAseq/query_clean.h5Seurat', verbose = F)
query_transferred <- LoadH5Seurat('../data/AAV_MG_scRNAseq/query_transferred.h5Seurat', verbose = F)
queryaav_sub <- subset(querydata, subset = group == "AAV")
query_transduced <- subset(queryaav_sub, subset = transduction == "Transduced")

Making Plots

Extended Data Figure 8a

features_hom = c("P2ry12","P2ry13","Cx3cr1","Csf1r","Lrp1","Ptpre","Lipa","Ets1")
features_hom_lab = c("Sall1", "Sirpa")

vplots <- list()
for(feature in features_hom){
  vplots[[feature]] <- FetchData(queryaav_sub, vars = c("treatment", "transduction", feature, "mScarlet_norm")) %>%
    ggplot(aes_string(x = "treatment", y = feature)) + 
    geom_jitter(aes_string(color = "mScarlet_norm"), size = 1.5, alpha = 1, width = 0.20) + 
    geom_violin(aes_string(fill = "treatment"), alpha = 0.3, width = 0.6) +
    scale_color_distiller('normalized\nmScarlet\nexpression', palette = "Reds", direction = 1) +
    scale_fill_manual(values = brewer.pal(n = 12, name = "Paired")[c(11,7)]) +
    guides(fill = "none") +
    theme_pubr() +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_blank(),
      axis.title.y = element_text(size = 15, family = "Arial"),
      axis.text.y = element_text(size = 12, family = "Arial"),
      legend.position = "none",
      legend.title = element_text(size = 10, family = "Arial"),
      legend.text = element_text(size = 9, family = "Arial")
    )
}

for(feature in features_hom_lab){
  vplots[[feature]] <- FetchData(queryaav_sub, vars = c("treatment", "transduction", feature, "mScarlet_norm")) %>%
    ggplot(aes_string(x = "treatment", y = feature)) + 
    geom_jitter(aes_string(color = "mScarlet_norm"), size = 1.5, alpha = 1, width = 0.20) + 
    geom_violin(aes_string(fill = "treatment"), alpha = 0.3, width = 0.6) +
    scale_color_distiller('normalized\nmScarlet\nexpression', palette = "Reds", direction = 1) +
    scale_fill_manual(values = brewer.pal(n = 12, name = "Paired")[c(11,7)]) +
    guides(fill = "none") +
    theme_pubr() +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_text(size = 12, family = "Arial"),
      axis.title.y = element_text(size = 15, family = "Arial"),
      axis.text.y = element_text(size = 12, family = "Arial"),
      legend.position = "none",
      legend.title = element_text(size = 10, family = "Arial"),
      legend.text = element_text(size = 9, family = "Arial"),
    )
}
ggarrange(vplots[["P2ry12"]], vplots[["P2ry13"]], vplots[["Cx3cr1"]], vplots[["Csf1r"]], vplots[["Sall1"]], vplots[["Sirpa"]], 
          ncol=2, nrow=3, common.legend = TRUE, legend="right")

Extended Data Figure 8b

Alternatively, we could calculate the Spearmann correlation between the expression levels of signature genes and mScarlet transgene. We did not find a correlative relationship between the signature genes and transgene, suggesting that AAV transduction does not induce microglia activation. These are not required by the reviewers but I think they do better at illustrating the points.

scplot <- list()
for(feature in c(features_hom, features_hom_lab)){
  scplot[[feature]] <- FetchData(query_transduced, vars = c(feature, "mScarlet_norm")) %>%
    ggscatter(
      x = "mScarlet_norm",
      y = feature,
      xlab = "Normalized mScarlet expression",
      add = c("reg.line"),
      add.params = list(color = "blue", linetype = "dashed"),
      size = .8
    ) +
    stat_cor(
      label.x.npc = 0.70,
      label.y.npc = 0.8,
      label.sep = "\n",
      color = "blue",
      cor.coef.name = "rho",
      method = "spearman",
      size = 2
    ) +
    theme(
      axis.title.x = element_text(size = 9),
      axis.text.x = element_text(size = 12),
      axis.title.y = element_text(size = 15),
      axis.text.y = element_text(size = 12),
    )
}
(scplot[["P2ry12"]] | scplot[["P2ry13"]]) / (scplot[["Cx3cr1"]] | scplot[["Csf1r"]]) / (scplot[["Sall1"]] | scplot[["Sirpa"]])

Extended Data Figure 8c

features_rea = c("Stat3","Gbp7","Ptgds","Irf7","Irf9","Saa3","Tlr2","Mt1","Mt2","Cd74","Lcn2","Tnf","Il1b","Cst7","Cxcl10",
                 "Ifitm3","Spp1")
features_rea_lab = c("Ccl5","Ccl2")

vplots <- list()
for(feature in features_rea){
  vplots[[feature]] <- FetchData(queryaav_sub, vars = c("treatment", "transduction", feature, "mScarlet_norm")) %>%
    ggplot(aes_string(x = "treatment", y = feature)) + 
    geom_jitter(aes_string(color = "mScarlet_norm"), size = 1.5, alpha = 1, width = 0.20) + 
    geom_violin(aes_string(fill = "treatment"), alpha = 0.3, width = 0.6) +
    scale_color_distiller('normalized\nmScarlet\nexpression', palette = "Reds", direction = 1) +
    scale_fill_manual(values = brewer.pal(n = 12, name = "Paired")[c(11,7)]) +
    guides(fill = "none") +
    theme_pubr() +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_blank(),
      axis.title.y = element_text(size = 15, family = "Arial"),
      axis.text.y = element_text(size = 12, family = "Arial"),
      legend.position = "none",
      legend.title = element_text(size = 10, family = "Arial"),
      legend.text = element_text(size = 9, family = "Arial")
    )
}

for(feature in features_rea_lab){
  vplots[[feature]] <- FetchData(queryaav_sub, vars = c("treatment", "transduction", feature, "mScarlet_norm")) %>%
    ggplot(aes_string(x = "treatment", y = feature)) + 
    geom_jitter(aes_string(color = "mScarlet_norm"), size = 1.5, alpha = 1, width = 0.20) + 
    geom_violin(aes_string(fill = "treatment"), alpha = 0.3, width = 0.6) +
    scale_color_distiller('normalized\nmScarlet\nexpression', palette = "Reds", direction = 1) +
    scale_fill_manual(values = brewer.pal(n = 12, name = "Paired")[c(11,7)]) +
    guides(fill = "none") +
    theme_pubr() +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_text(size = 12, family = "Arial"),
      axis.title.y = element_text(size = 15, family = "Arial"),
      axis.text.y = element_text(size = 12, family = "Arial"),
      legend.position = "none",
      legend.title = element_text(size = 10, family = "Arial"),
      legend.text = element_text(size = 9, family = "Arial"),
    )
}
ggarrange(vplots[["Tnf"]], vplots[["Il1b"]], vplots[["Cd74"]], vplots[["Cst7"]], vplots[["Spp1"]], vplots[["Irf9"]], vplots[["Ccl2"]], vplots[["Ccl5"]], 
          ncol=2, nrow=4, common.legend = TRUE, legend="right")

Extended Data Figure 8d

scplot <- list()
for(feature in c(features_rea, features_rea_lab)){
  scplot[[feature]] <- FetchData(query_transduced, vars = c(feature, "mScarlet_norm")) %>%
    ggscatter(
      x = "mScarlet_norm",
      y = feature,
      xlab = "Normalized mScarlet expression",
      add = c("reg.line"),
      add.params = list(color = "blue", linetype = "dashed"),
      size = .8
    ) +
    stat_cor(
      label.x.npc = 0.70,
      label.y.npc = 0.8,
      label.sep = "\n",
      color = "blue",
      cor.coef.name = "rho",
      method = "spearman",
      size = 2
    ) +
    theme(
      axis.title.x = element_text(size = 9),
      axis.text.x = element_text(size = 12),
      axis.title.y = element_text(size = 15),
      axis.text.y = element_text(size = 12),
    )
}
(scplot[["Tnf"]] | scplot[["Il1b"]]) / (scplot[["Cd74"]] | scplot[["Cst7"]]) / (scplot[["Spp1"]] | scplot[["Irf9"]]) / (scplot[["Ccl2"]] | scplot[["Ccl5"]])

Session info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=zh_CN.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=zh_CN.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=zh_CN.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3    patchwork_1.1.1       ggpubr_0.4.0         
##  [4] forcats_0.5.1         stringr_1.4.0         dplyr_1.0.9          
##  [7] purrr_0.3.4           readr_2.1.2           tidyr_1.2.0          
## [10] tibble_3.1.7          ggplot2_3.3.6         tidyverse_1.3.1      
## [13] SeuratDisk_0.0.0.9020 sp_1.5-0              SeuratObject_4.1.0   
## [16] Seurat_4.1.1         
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.0          backports_1.4.1       plyr_1.8.7           
##   [4] igraph_1.3.2          lazyeval_0.2.2        splines_4.0.5        
##   [7] listenv_0.8.0         scattermore_0.8       digest_0.6.29        
##  [10] htmltools_0.5.2       fansi_1.0.3           magrittr_2.0.3       
##  [13] tensor_1.5            cluster_2.1.3         ROCR_1.0-11          
##  [16] tzdb_0.3.0            globals_0.15.0        modelr_0.1.8         
##  [19] matrixStats_0.62.0    spatstat.sparse_2.1-1 colorspace_2.0-3     
##  [22] rvest_1.0.2           ggrepel_0.9.1         haven_2.5.0          
##  [25] xfun_0.31             crayon_1.5.1          jsonlite_1.8.0       
##  [28] progressr_0.10.1      spatstat.data_2.2-0   survival_3.3-1       
##  [31] zoo_1.8-10            glue_1.6.2            polyclip_1.10-0      
##  [34] gtable_0.3.0          leiden_0.4.2          car_3.1-0            
##  [37] future.apply_1.9.0    abind_1.4-5           scales_1.2.0         
##  [40] DBI_1.1.2             rstatix_0.7.0         spatstat.random_2.2-0
##  [43] miniUI_0.1.1.1        Rcpp_1.0.8.3          viridisLite_0.4.0    
##  [46] xtable_1.8-4          reticulate_1.25       spatstat.core_2.4-4  
##  [49] bit_4.0.4             htmlwidgets_1.5.4     httr_1.4.3           
##  [52] ellipsis_0.3.2        ica_1.0-2             farver_2.1.0         
##  [55] pkgconfig_2.0.3       sass_0.4.1            uwot_0.1.11          
##  [58] dbplyr_2.2.0          deldir_1.0-6          utf8_1.2.2           
##  [61] labeling_0.4.2        tidyselect_1.1.2      rlang_1.0.2          
##  [64] reshape2_1.4.4        later_1.3.0           munsell_0.5.0        
##  [67] cellranger_1.1.0      tools_4.0.5           cli_3.3.0            
##  [70] generics_0.1.2        broom_0.8.0           ggridges_0.5.3       
##  [73] evaluate_0.15         fastmap_1.1.0         yaml_2.3.5           
##  [76] goftest_1.2-3         fs_1.5.2              knitr_1.39           
##  [79] bit64_4.0.5           fitdistrplus_1.1-8    RANN_2.6.1           
##  [82] pbapply_1.5-0         future_1.26.1         nlme_3.1-158         
##  [85] mime_0.12             xml2_1.3.3            hdf5r_1.3.5          
##  [88] compiler_4.0.5        rstudioapi_0.13       plotly_4.10.0        
##  [91] png_0.1-7             ggsignif_0.6.3        spatstat.utils_2.3-1 
##  [94] reprex_2.0.1          bslib_0.3.1           stringi_1.7.6        
##  [97] highr_0.9             rgeos_0.5-9           lattice_0.20-45      
## [100] Matrix_1.4-1          vctrs_0.4.1           pillar_1.7.0         
## [103] lifecycle_1.0.1       spatstat.geom_2.4-0   lmtest_0.9-40        
## [106] jquerylib_0.1.4       RcppAnnoy_0.0.19      data.table_1.14.2    
## [109] cowplot_1.1.1         irlba_2.3.5           httpuv_1.6.5         
## [112] R6_2.5.1              promises_1.2.0.1      KernSmooth_2.23-20   
## [115] gridExtra_2.3         parallelly_1.32.0     codetools_0.2-18     
## [118] MASS_7.3-57           assertthat_0.2.1      withr_2.5.0          
## [121] sctransform_0.3.3     mgcv_1.8-40           parallel_4.0.5       
## [124] hms_1.1.1             grid_4.0.5            rpart_4.1-15         
## [127] rmarkdown_2.14        carData_3.0-5         Rtsne_0.16           
## [130] shiny_1.7.1           lubridate_1.8.0