Cell quality control

This is the analysis of two scRNAseq dataset containing WT or Gmnc KO P0 cortex + hippocampus. 4 WT and 4 KO were used. Cells were prepared by Vicente Elorriaga & Frédéric Causeret
Libraries were generated by the Imagine scRNAseq platform
Sequencing was achieved at the genomics platform of Imagine
Reads were aligned on the mm10 genome

Load libraries

library(Seurat)
library(cowplot)
library(dplyr)
library(ggplot2)
library(ggExtra)
library(ggrepel)
library(reticulate)
library(Matrix)
library(viridis)
library(RColorBrewer)
library(MetBrewer)
library(wesanderson)
library(R.utils)

# Set ggplot theme as classic
theme_set(theme_classic())

Load the dataset and calculate QC metrics

Initialize a Seurat object from the raw filtered gene/bc matrix

# Load the WT filtered_gene_bc_matrix outputed by Cell Ranger v7.1 (including intronic reads!)
Countdata <- Read10X(data.dir = "./outs/P0_Gmnc_WT/filtered_feature_bc_matrix/")

# Initialize the Seurat object
Gmnc.WT <- CreateSeuratObject(counts = Countdata,
                               min.cells = 3,
                               min.features = 800,
                               project = "Gmnc.WT")
Gmnc.WT$Barcodes <- colnames(x=Gmnc.WT)

# Load the KO filtered_gene_bc_matrix outputed by Cell Ranger v7.1 (including intronic reads!)
Countdata <- Read10X(data.dir = "./outs/P0_Gmnc_KO/filtered_feature_bc_matrix/")

# Initialize the Seurat object
Gmnc.KO <- CreateSeuratObject(counts = Countdata,
                               min.cells = 3,
                               min.features = 800,
                               project = "Gmnc.KO")
Gmnc.KO$Barcodes <- colnames(x=Gmnc.KO)

dim(Gmnc.WT)

## [1] 20768 11363

dim(Gmnc.KO)

## [1] 21383 25346

rm("Countdata")

Calculate percentage of mitochondrial and ribosomal counts

# Percent of mitochondrial counts
Gmnc.WT[["percent.mt"]] <- PercentageFeatureSet(Gmnc.WT, pattern = "^mt-")
Gmnc.KO[["percent.mt"]] <- PercentageFeatureSet(Gmnc.KO, pattern = "^mt-")
# Percent of ribosomal counts
Gmnc.WT[["percent.rb"]] <- PercentageFeatureSet(Gmnc.WT, pattern = "(^Rpl|^Rps|^Mrp)")
Gmnc.KO[["percent.rb"]] <- PercentageFeatureSet(Gmnc.KO, pattern = "(^Rpl|^Rps|^Mrp)")

Cell Quality according to total number and % mitochondrial reads

# Filter WT cells with more than 10% mito reads and/or less than 3000 reads
VlnPlot(Gmnc.WT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), ncol = 4, pt.size = 0.2)

Gmnc.WT@meta.data$Cell.quality <- ifelse(Gmnc.WT@meta.data$percent.mt>10 & Gmnc.WT@meta.data$nCount_RNA<3000,
                                         "Low Quality",
                                         ifelse(Gmnc.WT@meta.data$percent.mt>10,
                                                "High.mt",
                                                ifelse(Gmnc.WT@meta.data$nCount_RNA<3000,
                                                       "Low.counts", "High.quality")))
table(Gmnc.WT$Cell.quality)

## 
##      High.mt High.quality  Low Quality   Low.counts 
##          272         7646          779         2666

# Filter KO cells with more than 5% mito reads and/or less than 3000 reads

VlnPlot(Gmnc.KO, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), ncol = 4, pt.size = 0.2)

Gmnc.KO@meta.data$Cell.quality <- ifelse(Gmnc.KO@meta.data$percent.mt>10 & Gmnc.KO@meta.data$nCount_RNA<3000,
                                         "Low Quality",
                                         ifelse(Gmnc.KO@meta.data$percent.mt>10,
                                                "High.mt",
                                                ifelse(Gmnc.KO@meta.data$nCount_RNA<3000,
                                                       "Low.counts", "High.quality")))
table(Gmnc.KO$Cell.quality)

## 
##      High.mt High.quality  Low Quality   Low.counts 
##            4        19632            8         5702

Plot WT QC metrics

# Relation between nCount_RNA and nFeatures_RNA detected with cell quality parameter
p1 <- ggplot(Gmnc.WT@meta.data, aes(x=nCount_RNA, y=nFeature_RNA)) + geom_point(aes(color=Cell.quality), size=0.1) 
p1 <- ggMarginal(p1, type = "histogram", fill="lightgrey")

p2 <- ggplot(Gmnc.WT@meta.data, aes(x=log10(nCount_RNA), y=log10(nFeature_RNA))) + geom_point(aes(color=Cell.quality), size=0.1) + geom_smooth(method="lm")
p2 <- ggMarginal(p2, type = "histogram", fill="lightgrey")

# Relation between nFeatures_RNA and the mitochondrial RNA percentage detected with cell quality parameter
p3 <- ggplot(Gmnc.WT@meta.data, aes(x=nFeature_RNA, y=percent.mt, color=Cell.quality)) + geom_point(size=0.1)
p3 <- ggMarginal(p3, type = "histogram", fill="lightgrey", bins=100) 
    
plot_grid(plotlist = list(p1,p2,p3), ncol=3, align='h', rel_widths = c(1, 1, 1))

Plot KO QC metrics

# Relation between nCount_RNA and nFeatures_RNA detected with cell quality parameter
p1 <- ggplot(Gmnc.KO@meta.data, aes(x=nCount_RNA, y=nFeature_RNA)) + geom_point(aes(color=Cell.quality), size=0.1) 
p1 <- ggMarginal(p1, type = "histogram", fill="lightgrey")

p2 <- ggplot(Gmnc.KO@meta.data, aes(x=log10(nCount_RNA), y=log10(nFeature_RNA))) + geom_point(aes(color=Cell.quality), size=0.1) + geom_smooth(method="lm")
p2 <- ggMarginal(p2, type = "histogram", fill="lightgrey")

# Relation between nFeatures_RNA and the mitochondrial RNA percentage detected with cell quality parameter
p3 <- ggplot(Gmnc.KO@meta.data, aes(x=nFeature_RNA, y=percent.mt, color=Cell.quality)) + geom_point(size=0.1)
p3 <- ggMarginal(p3, type = "histogram", fill="lightgrey", bins=100) 
    
plot_grid(plotlist = list(p1,p2,p3), ncol=3, align='h', rel_widths = c(1, 1, 1))

Cell filtering

# Subset the data
Gmnc.WT <- subset(x = Gmnc.WT, subset =  Cell.quality == "High.quality")
dim(Gmnc.WT)

## [1] 20768  7646

Gmnc.KO <- subset(x = Gmnc.KO, subset =  Cell.quality == "High.quality")
dim(Gmnc.KO)

## [1] 21383 19632

Normalization

# Subset the data
Gmnc.WT <- NormalizeData(Gmnc.WT, normalization.method = "LogNormalize", scale.factor = 10000)
Gmnc.KO <- NormalizeData(Gmnc.KO, normalization.method = "LogNormalize", scale.factor = 10000)

Cell Cycle Scoring

# Assign cell-cycle scores
s.genes <- c("Mcm5", "Pcna", "Tym5", "Fen1", "Mcm2", "Mcm4", "Rrm1", "Ung", "Gins2", "Mcm6", "Cdca7", "Dtl", "Prim1", "Uhrf1", "Mlf1ip", "Hells", "Rfc2", "Rap2", "Nasp", "Rad51ap1", "Gmnn", "Wdr76", "Slbp", "Ccne2", "Ubr7", "Pold3", "Msh2", "Atad2", "Rad51", "Rrm2", "Cdc45", "Cdc6", "Exo1", "Tipin", "Dscc1", "Blm", " Casp8ap2", "Usp1", "Clspn", "Pola1", "Chaf1b", "Brip1", "E2f8")

g2m.genes <- c("Hmgb2", "Ddk1","Nusap1", "Ube2c", "Birc5", "Tpx2", "Top2a", "Ndc80", "Cks2", "Nuf2", "Cks1b", "Mki67", "Tmpo", " Cenpk", "Tacc3", "Fam64a", "Smc4", "Ccnb2", "Ckap2l", "Ckap2", "Aurkb", "Bub1", "Kif11", "Anp32e", "Tubb4b", "Gtse1", "kif20b", "Hjurp", "Cdca3", "Hn1", "Cdc20", "Ttk", "Cdc25c", "kif2c", "Rangap1", "Ncapd2", "Dlgap5", "Cdca2", "Cdca8", "Ect2", "Kif23", "Hmmr", "Aurka", "Psrc1", "Anln", "Lbr", "Ckap5", "Cenpe", "Ctcf", "Nek2", "G2e3", "Gas2l3", "Cbx5", "Cenpa")

Gmnc.WT <- CellCycleScoring(Gmnc.WT,
                             s.features = s.genes,
                             g2m.features = g2m.genes,
                             set.ident = T)

Gmnc.KO <- CellCycleScoring(Gmnc.KO,
                             s.features = s.genes,
                             g2m.features = g2m.genes,
                             set.ident = T)
table(Gmnc.WT$Phase)

## 
##   G1  G2M    S 
## 5430 1254  962

table(Gmnc.KO$Phase)

## 
##    G1   G2M     S 
## 13512  3382  2738

rm(s.genes, g2m.genes)

Gmnc.WT$CC.Difference <- Gmnc.WT$S.Score - Gmnc.WT$G2M.Score
Gmnc.KO$CC.Difference <- Gmnc.KO$S.Score - Gmnc.KO$G2M.Score

Merge datasets and perform dimensionality reduction

Gmnc.combined <- merge(Gmnc.WT, y = Gmnc.KO, add.cell.ids = c("WT", "KO"), project = "P0_Gmnc")
rm(Gmnc.KO, Gmnc.WT)

Gmnc.combined <- FindVariableFeatures(Gmnc.combined, selection.method = "vst", nfeatures = 2000)
Gmnc.combined <- ScaleData(Gmnc.combined)
Gmnc.combined <- RunPCA(Gmnc.combined, features = VariableFeatures(object = Gmnc.combined))
Gmnc.combined <- RunUMAP(Gmnc.combined, dims = 1:20)
Gmnc.combined <- JoinLayers(Gmnc.combined)

Plot UMAP projection

DimPlot(Gmnc.combined, reduction = "umap", label = F, label.size = 2, pt.size = 0.1, group.by = "orig.ident", cols = met.brewer("Egypt", 4)) + NoAxes() +ggtitle(paste0(table(Gmnc.combined$orig.ident)[2], " WT cells + ", table(Gmnc.combined$orig.ident)[1], " KO cells"))

Cluster the cells

Gmnc.combined <- FindNeighbors(Gmnc.combined, dims = 1:20)
Gmnc.combined <- FindClusters(Gmnc.combined, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 27278
## Number of edges: 915606
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9268
## Number of communities: 23
## Elapsed time: 5 seconds

DimPlot(Gmnc.combined, reduction = "umap", label = T, label.size = 4, pt.size = 0.1, group.by = "seurat_clusters",  cols = met.brewer("Klimt", 23)) + NoAxes()

FeaturePlot(Gmnc.combined, features=c("Trp73", "Tbr1", "Gad2", "Pdgfra", "Vim", "Foxc1", "C1qa"),ncol=3, reduction = "umap", order = T) & scale_color_gradientn(colors=c("grey90", brewer.pal(9,"YlGnBu"))) & NoLegend() & NoAxes()

# Save Seurat object

saveRDS(Gmnc.combined, "Gmnc.combined.RDS")

Session Info

#date
format(Sys.time(), "%d %B, %Y, %H,%M")

## [1] "08 January, 2024, 22,08"

#Packages used
sessionInfo()

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.1.1/lib/libopenblasp-r0.3.18.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] R.utils_2.11.0     R.oo_1.24.0        R.methodsS3_1.8.1  wesanderson_0.3.6 
##  [5] MetBrewer_0.2.0    RColorBrewer_1.1-3 viridis_0.6.2      viridisLite_0.4.1 
##  [9] Matrix_1.6-4       reticulate_1.24    ggrepel_0.9.1      ggExtra_0.9       
## [13] ggplot2_3.3.6      dplyr_1.0.10       cowplot_1.1.1      Seurat_5.0.1      
## [17] SeuratObject_5.0.1 sp_2.1-2          
## 
## loaded via a namespace (and not attached):
##   [1] spam_2.7-0             plyr_1.8.7             igraph_1.3.1          
##   [4] lazyeval_0.2.2         splines_4.1.1          RcppHNSW_0.3.0        
##   [7] listenv_0.8.0          scattermore_1.2        digest_0.6.30         
##  [10] htmltools_0.5.2        fansi_1.0.3            magrittr_2.0.3        
##  [13] tensor_1.5             cluster_2.1.2          ROCR_1.0-11           
##  [16] globals_0.14.0         matrixStats_0.62.0     spatstat.sparse_3.0-3 
##  [19] colorspace_2.0-3       xfun_0.34              crayon_1.5.2          
##  [22] jsonlite_1.8.2         progressr_0.10.0       spatstat.data_3.0-3   
##  [25] survival_3.2-13        zoo_1.8-10             glue_1.6.2            
##  [28] polyclip_1.10-0        gtable_0.3.1           leiden_0.3.10         
##  [31] future.apply_1.9.0     abind_1.4-5            scales_1.2.1          
##  [34] DBI_1.1.2              spatstat.random_3.2-2  miniUI_0.1.1.1        
##  [37] Rcpp_1.0.9             xtable_1.8-4           dotCall64_1.0-1       
##  [40] htmlwidgets_1.5.4      httr_1.4.4             ellipsis_0.3.2        
##  [43] ica_1.0-2              pkgconfig_2.0.3        farver_2.1.1          
##  [46] sass_0.4.1             uwot_0.1.11            deldir_1.0-6          
##  [49] utf8_1.2.2             tidyselect_1.2.0       labeling_0.4.2        
##  [52] rlang_1.0.6            reshape2_1.4.4         later_1.3.0           
##  [55] munsell_0.5.0          tools_4.1.1            cli_3.4.1             
##  [58] generics_0.1.3         ggridges_0.5.3         evaluate_0.17         
##  [61] stringr_1.4.1          fastmap_1.1.0          yaml_2.3.6            
##  [64] goftest_1.2-3          knitr_1.40             fitdistrplus_1.1-8    
##  [67] purrr_0.3.5            RANN_2.6.1             pbapply_1.5-0         
##  [70] future_1.25.0          nlme_3.1-153           mime_0.12             
##  [73] ggrastr_1.0.1          compiler_4.1.1         rstudioapi_0.13       
##  [76] beeswarm_0.4.0         plotly_4.10.0          png_0.1-7             
##  [79] spatstat.utils_3.0-4   tibble_3.1.8           bslib_0.3.1           
##  [82] stringi_1.7.8          highr_0.9              RSpectra_0.16-1       
##  [85] lattice_0.20-45        vctrs_0.4.2            pillar_1.8.1          
##  [88] lifecycle_1.0.3        spatstat.geom_3.2-7    lmtest_0.9-40         
##  [91] jquerylib_0.1.4        RcppAnnoy_0.0.19       data.table_1.14.2     
##  [94] irlba_2.3.5.1          httpuv_1.6.5           patchwork_1.1.1       
##  [97] R6_2.5.1               promises_1.2.0.1       KernSmooth_2.23-20    
## [100] gridExtra_2.3          vipor_0.4.5            parallelly_1.31.1     
## [103] codetools_0.2-18       fastDummies_1.7.3      MASS_7.3-54           
## [106] assertthat_0.2.1       withr_2.5.0            sctransform_0.4.1     
## [109] mgcv_1.8-38            parallel_4.1.1         grid_4.1.1            
## [112] tidyr_1.2.1            rmarkdown_2.11         Rtsne_0.16            
## [115] spatstat.explore_3.2-5 shiny_1.7.1            ggbeeswarm_0.6.0

---

1. IPNP & Imagine Institute, Paris, France, frederic.causeret@inserm.fr↩︎