1 PREAMBLE

1.1 Purpose of this session

This file describes the different steps to perform the third part of the single cell RNAseq data analysis training course for the EBAII n1 2024, covering these steps :

  • Filtering of bad cells and extreme low-expression features

  • Estimation of the cell cycle phase

  • Identification and removal of cell doublets

2 Start Rstudio

3 Warm-up

  • We set common parameters we will use throughout this session :
## Set your project name
project_name <- "ebaii_sc_teachers"  # Do not copy-paste this ! It's MY project ! Put yours !!

## Seed for the RNG
my_seed <- 1337L

## Minimum cells with counts to keep a feature
min_cells <- 5

4 Prepare the data structure

We will do the same as for former steps, just changing the session name :

4.1 Main directory

## Preparing the path
TD_dir <- paste0("/shared/projects/", project_name, "/SC_TD")

## Creating the root directory
dir.create(path = TD_dir, recursive = TRUE)

## Print the root directory on-screen
print(TD_dir)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD"

4.2 Current session

## Creating the session (Preproc.3) directory
session_dir <- paste0(TD_dir, "/03_Preproc.3")
dir.create(path = session_dir, recursive = TRUE)

## Print the session directory on-screen
print(session_dir)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/03_Preproc.3"

4.3 Input directory

## Creating the INPUT data directory
input_dir <- paste0(session_dir, "/DATA")
dir.create(path = input_dir, recursive = TRUE)

## Print the input directory on-screen
print(input_dir)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/03_Preproc.3/DATA"

4.4 Genelists directory

res_dir <- paste0(TD_dir, "/Resources")
glist_dir <- paste0(res_dir, "/Genelists")

## Print the resources directory on-screen
print(glist_dir)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/Resources/Genelists"

4.5 Output directory

## Creating the OUTPUT data directory
output_dir <- paste0(session_dir, "/RESULTS")
dir.create(path = output_dir, recursive = TRUE)

## Print the output directory on-screen
print(output_dir)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/03_Preproc.3/RESULTS"

5 Reload the Seurat Object

  • We can reload the object we saved at the former step
## The latest Seurat object saved as RDS (name)
sobj_file <- "02_TD3A_S5_Metrics.Bio_31053.4587.RDS"

## The latest Seurat object saved as RDS (full path)
sobj_path <- paste0(TD_dir, 
                    "/02_Preproc.2/RESULTS/",
                    sobj_file)

force <- FALSE  ## To force a re-download of a Zenodo-hosted backup
local <- FALSE  ## To force a loading from a local backup

## In case of error/lost data : force a reload from a Zenodo backup repository
if(force) {
  zen_id <- "14035293"
  zen_backup_file <- paste0("https://zenodo.org/records/",
                            zen_id,
                            "/files/",
                            sobj_file)
  download.file(url = zen_backup_file,
                destfile = sobj_path)
}

## In case of error/lost data : force a reload from a local backup repository
if(local) {
  sobj_path <- paste0(
    "/shared/projects/2422_ebaii_n1/atelier_scrnaseq/TD/BACKUP/RDS/",
    sobj_file)
}

## Load the object
sobj <- readRDS(file = sobj_path)



Now, we finally have all the metrics and bias sources we could use, so we can actually get to the filtering step.

6 Filtering

6.1 Features

  • As already explained, the sparsity of the count matrix is high, very high.

  • To the point that there probably are some features that have so low expression level that they were only counted in very few cells. As such, they have absolutely no chance to characterize any cell type, nor harbor some variation between different cell types.

  • As such, we should discard these features


What are the actual dimensions of our expression data ?

## Dimensions of our Seurat object
dim(sobj)
Show output 
[1] 31053  4587


We want to remove the features that are expressed in less than 5 cells.

First, we quantify, for each feature, the number of cells with 0 count.

## Compute the amount of 0s per feature
nFeat_zero <- sparseMatrixStats::rowCounts(
  x = SeuratObject::GetAssayData(
    object = sobj, 
    assay = "RNA", 
    layer = "counts"), 
  value = 0)

## Inversion : computing the number of cells with at least one count, per feature
nFeat_nonzero <- ncol(sobj) - nFeat_zero

## Identify those with at least 5 cells with expression
nFeat_keep <- nFeat_nonzero >= min_cells

## Quantify the selected ones
table(nFeat_keep)
Show output 
nFeat_keep
FALSE  TRUE 
18545 12508


And now we can discard features expressed in less than 5 cells :

filtc5_pre <- SC.helper::QnD_viz(sobj = sobj, 
                                return_plot = TRUE)
Show plot 


## Restrict the Seurat object
sobj <- sobj[nFeat_keep,]


What are the Seurat object dimensions, now ?

dim(sobj)
Show output 
[1] 12508  4587


We can visualize the cell space after this features filtering

filtc5_post <- SC.helper::QnD_viz(sobj = sobj, 
                                  return_plot = TRUE)
Show plot


Merging plots for ease of visualization :

## Using the patchwork package to merge plots (and ggplot2 to add titles)
patchwork::wrap_plots(
  list(
    filtc5_pre & ggplot2::ggtitle(label = "BEFORE"),
    filtc5_post & ggplot2::ggtitle(label = "AFTER")), 
  nrow = 1)
Show plot


  • Question :

    Do you see any difference when comparing before vs after features filtering ?


    ## . Not much changed, maybe cell groups seem a bit more condensed
##
## . This is expected, as we removed features with almost no expression



6.2 Cells

We are now able to apply all the filtering strategies we established for each QC metric.

  • Identify “good” cells
## Identify the cells to keep
bc_keep <-  sobj$nFeature_RNA_in_range & 
            sobj$nCount_RNA_in_range & 
            sobj$percent_mt_in_range & 
            sobj$percent_rb_in_range & 
            sobj$percent_st_in_range

## Contengency
table(bc_keep)
Show output 
bc_keep
FALSE  TRUE 
  309  4278


  • Visualize metrics filtering effect (as an upset-plot)

    ## Create a list of all cells filtered ou for each criterion
up_list <-list(
  "nFeature" = colnames(sobj)[!sobj$nFeature_RNA_in_range],
  "nCount" = colnames(sobj)[!sobj$nCount_RNA_in_range],
  "%MITO" = colnames(sobj)[!sobj$percent_mt_in_range],
  "%RIBO" = colnames(sobj)[!sobj$percent_rb_in_range],
  "%STRESS" = colnames(sobj)[!sobj$percent_st_in_range])

## Create an upset-plot
UpSetR::upset(data = UpSetR::fromList(up_list), 
              nintersects = NA, 
              sets = rev(names(up_list)),
              keep.order = TRUE,
              order.by = "freq")

    Show plot


  • Apply the filter

    ## Seurat object BEFORE cell filtering
dim(sobj)

    Show output

     
    [1] 12508  4587
    ## Apply cell filtering
sobj <- subset(x = sobj, cells = colnames(sobj)[bc_keep])

## Seurat object AFTER cell filtering
dim(sobj)

    Show output

     
    [1] 12508  4278


We can visualize the cell space since this cell filtering

SC.helper::QnD_viz(sobj = sobj)
Show plot


  • Question

    Do you see any difference when comparing before vs after features filtering ?


    ## . Not much changed as well (we did not discard many cells)
##
## . The biggest cluster structure seems more defined


7 Save the Seurat object

We will save our Seurat object that now contains filtered cells and features :

## Save our Seurat object (rich naming)
out_name <- paste0(
          output_dir, "/", paste(
            c("03", Seurat::Project(sobj), "S5", 
              "Metrics.Filtered", paste(
                dim(sobj), 
                collapse = '.'
              )
            ), collapse = "_"),
            ".RDS")

## Check
print(out_name)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/03_Preproc.3/RESULTS/03_TD3A_S5_Metrics.Filtered_12508.4278.RDS"
## Write on disk
saveRDS(object = sobj, 
        file = out_name)


8 Cell cycle scores

  • We are currently analyzing independent profiles from isolated cells, from a sample dissociation

  • As such, cells were most probably not synchronized, thus the effect of their cell cycle state on genes expression may be strong, to the point that it can bias the data (ie, mask some lower amplitude biological variation).

  • In order to assess (and maybe, remove) this bias, we have to quantify it.

  • We will perform this estimation thanks to heuristics based on knowledge : Seurat includes a method that evaluates the cell cycle phase of cells through scores for the S and G2M phases, each based on phase-specific gene signatures.

  • For this step, we will use additional gene lists from knowledge (cell cycle phase), hosted in a Zenodo respository (Id : 14037355)

8.1 Download gene lists

  • We will directly retrieve data from Zenodo to your input_dir :

    ## Zenodo ID
zen_id <- '14101506'

### Named files (will be used later on !)
cc_file <- "mus_musculus_Seurat_cc.genes_20191031.rds"

## Filename(s) to retrieve
toget_files <- c(cc_file)

## Folder to store retrieved files
local_folder <- glist_dir

## Use local backup ?
backup <- FALSE
if(backup) message("Using local backup !")

## Force download ?
force <- FALSE
if(force) message("Forcing (re)download !")

### Define remote folder
remote_folder <- if (backup) {
  "/shared/projects/2422_ebaii_n1/atelier_scrnaseq/TD/RESOURCES/GENELISTS/"
}  else {
  paste0("https://zenodo.org/records/", zen_id, "/files/")
}

### Reconstruct the input paths
remote_path <- paste0(remote_folder, "/", toget_files)

### Reconstruct the output paths
local_path <- paste0(local_folder, "/", toget_files)

## Retrieve files (if they don't exist), in loop
for (tg in seq_along(toget_files)) {
  ## If the file does not locally exist
  if (!file.exists(local_path[tg]) | force) {
    ## Retrieve data
    if(backup) {
      file.copy(from = remote_path[tg],
                to = local_path[tg])
    } else {
      download.file(url = remote_path[tg], 
                    destfile = local_path[tg])
    }
    ## Check if downloaded files exist locally
    if(file.exists(local_path[tg])) message("\tOK")
  } else message(paste0(toget_files[tg], " already downloaded !"))
}


8.2 Load gene lists

## The cell cycle gene lists file
cc_file <- paste0(glist_dir, 
                  "/", 
                  cc_file)

## Load the cell cycle reference genes lists
cc_genes <- readRDS(file = cc_file)

## Have a look on its content
str(cc_genes)
Show output 
List of 2
 $ s.genes  : chr [1:43] "Mcl5" "Pcna" "Tyms" "Fen1" ...
 $ g2m.genes: chr [1:54] "Hmgb2" "Cdk1" "Nusap1" "Ube2c" ...


  • Question

    As explained, we will use gene lists extracted from community knowledge. 
Our data contain values for genes also, so we will cross them. 
Is there something we should check ?


    ## . We may check if the genes in our gene lists are effectively 
##   present in our Seurat object !
##
## . This is expected, as we removed features with almost no expression


Beyond :

  1. Write a code that performs this check

  2. Add a code to adjust the content of the gene lists accordingly (remove genes from the gene lists that are not present in our dataset)

  3. (NOTE : this is actually not needed as already checked and corrected by the cell cyle estimation method we will use)


## Check if our data genes cover the gene lists
lapply(cc_genes, function(x) x %in% rownames(sobj))

## Remove genes not in sobj
cc_genes <- lapply(cc_genes, function(gl) { gl[gl %in% rownames(sobj)] })

## Check the modification
str(cc_genes)

## Check that all genes are available in sobj
all(unique(unlist(cc_genes)) %in% rownames(sobj))



8.3 Estimation

  • Let’s perform this estimation. But how ?

    ## Reading the function help page
?Seurat::CellCycleScoring


  • We will actually use a helper function to ease up the process :

    ?SC.helper::CC_Seurat




  • Run the cell-cycle estimation

    ## Perform the estimation
## The RNG seed is needed here !
sobj <- SC.helper::CC_Seurat(
  sobj = sobj, 
  assay = "RNA",
  seurat_cc_genes = cc_genes, 
  SmG2M = TRUE, 
  nbin = 20, 
  my_seed = my_seed)


  • Description of the object to see the data added

    SC.helper::SeuratObject_descriptor(
  sobj = sobj,
  describe = "coldata")

    Show output

     
    OBJECT VERSION :    5.0.2 
PROJECT :   [TD3A] 

[BARCODES METADATA]
orig.ident    Freq
-----------  -----
TD3A          4278
NA               0

 nCount_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1000    2054    2446    3737    3213   48866 

 nFeature_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    750    1325    1498    1697    1752    5968 

 log10_nCount_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.000   3.313   3.389   3.458   3.507   4.689 
nCount_RNA_in_range    Freq
--------------------  -----
TRUE                   4278
NA                        0
nFeature_RNA_in_range    Freq
----------------------  -----
TRUE                     4278
NA                          0

 percent_mt 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001801 0.019618 0.024355 0.025143 0.029851 0.050000 

 percent_rb 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02687 0.08067 0.09697 0.11161 0.12280 0.42010 

 percent_st 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01546 0.02991 0.03365 0.03422 0.03773 0.05980 
percent_mt_in_range    Freq
--------------------  -----
TRUE                   4278
NA                        0
percent_rb_in_range    Freq
--------------------  -----
TRUE                   4278
NA                        0
percent_st_in_range    Freq
--------------------  -----
TRUE                   4278
NA                        0

 CC_Seurat_S.Score 
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.2338439 -0.1060737 -0.0584921 -0.0184198 -0.0005651  1.2329054 

 CC_Seurat_G2M.Score 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.22728 -0.11143 -0.07202 -0.03641 -0.02488  1.59231 
CC_Seurat_Phase    Freq
----------------  -----
G1                 2878
G2M                 471
S                   929
NA                    0

 CC_Seurat_SmG2M.Score 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1.56874 -0.04569  0.01532  0.01799  0.08286  0.96435


8.4 Visualization

As usual, we can visualize the results as violins :

Seurat::VlnPlot(object = sobj,
                features = c("CC_Seurat_S.Score",
                             "CC_Seurat_G2M.Score",
                             "CC_Seurat_SmG2M.Score"))
Show plot


But it’s not that easy to interpret…

Let’s plot it in the cell space

8.4.1 “S minus G2M” (SmG2M) score :

## Using the 'features' parameter to plot the SmG2M score (continuous data)
## Here, we will keep the modified Seurat object to speed up further plots
VIZ <- SC.helper::QnD_viz(sobj = sobj, 
                          features = "CC_Seurat_SmG2M.Score", 
                          return_object = TRUE)
Show plot


8.4.2 Estimated cell phase :

## Using the 'group_by' parameter to plot the estimated phases (categorical data)
## Here, we recycle keep the VIZ Seurat object that contains everything to perform the plot without computing it again
SC.helper::QnD_viz(sobj = VIZ, 
                   slot = NULL, 
                   dimred = "umap", 
                   group_by = "CC_Seurat_Phase")
Show plot 

## VIZ is not needed anymore
rm(VIZ)


Questions

Does the structure of the cells in this representation seem to have a link with these cell cycle phases/scores ?


## It's an absolute yes !


Do you think this is the result of an artifact, or
biology-related ?


##
##    ¯\_(ツ)_/¯
##


9 Save the Seurat object

We will save our Seurat object that now contains the cell cycle states/scores :

## Save our Seurat object (rich naming)
out_name <- paste0(
          output_dir, "/", paste(
            c("04", Seurat::Project(sobj), "S5", 
              "CC", paste(
                dim(sobj), 
                collapse = '.'
              )
            ), collapse = "_"),
            ".RDS")

## Check
print(out_name)
[1] "/shared/projects/ebaii_sc_teachers/SC_TD/03_Preproc.3/RESULTS/04_TD3A_S5_CC_12508.4278.RDS"
## Write on disk
saveRDS(object = sobj, 
        file = out_name)



10 Cell doublets

10.1

We will use two different methods to detect and remove cell doublets :

  • scds (in its “hybrid” mode){target=“_blank”} : more efficient at detecting homotypic doublets

  • scDblFinder : more efficient at detecting heterotypic doublets

None of the methods accepts a SeuratObject as input, but a SingleCellExperiment object.

Hopefully :

  • Seurat has a function to perform the conversion
  • The output results can be integrated into our Seurat object with ease

10.2 scds

## Fix seed
set.seed(my_seed)

## Run scds
sobj$doublet_scds.hybrid <- unname(
  scds::cxds_bcds_hybrid(
    Seurat::as.SingleCellExperiment(
      sobj, assay = "RNA"))$hybrid_score > 1)

## Contingencies
table(sobj$doublet_scds.hybrid)
Show output 

FALSE  TRUE 
 4072   206

10.3 scDblFinder

## Fix seed
set.seed(my_seed)

## Run scDblFinder (which needs another object type)
sobj$doublet_scDblFinder <- scDblFinder::scDblFinder(
  sce = Seurat::as.SingleCellExperiment(
    x = sobj, 
    assay = "RNA"),
  returnType = "table")$class == "doublet"

## Contingencies
table(sobj$doublet_scDblFinder)
Show output 

FALSE  TRUE 
 4100   178


10.4 Merge results

We merge results of the two methods

## Logical union of both methods
sobj$doublet_union <- sobj$doublet_scds.hybrid | sobj$doublet_scDblFinder

## Quantify doublets
table(sobj$doublet_union)
Show output 

FALSE  TRUE 
 4035   243


We can assess tool-specific and common doublets

### Singlets by default
sobj$doublet_viz <- "singlet"

### Union
sobj$doublet_viz[sobj$doublet_union] <- "both"

### scds-specific
sobj$doublet_viz[sobj$doublet_scds.hybrid & !sobj$doublet_scDblFinder] <- "scds"

### scDblFinder-specific
sobj$doublet_viz[sobj$doublet_scDblFinder & !sobj$doublet_scds.hybrid] <- "scDblFinder"

## Convert to factor
sobj$doublet_viz <- as.factor(sobj$doublet_viz)

## Contingencies
table(sobj$doublet_viz)
Show output 

       both scDblFinder        scds     singlet 
        141          37          65        4035


Beyond : Create an upset-plot for the doublet status according to the two methods used

## Build a list of cells tagged by one tool, the other, or both
dbl_list <-list(
  "scds" = colnames(sobj)[sobj$doublet_scds.hybrid & !sobj$doublet_scDblFinder],
  "scDblFinder" = colnames(sobj)[!sobj$doublet_scds.hybrid & sobj$doublet_scDblFinder],
  "both" = colnames(sobj)[sobj$doublet_scds.hybrid & sobj$doublet_scDblFinder])

## Draw the upset plot
UpSetR::upset(data = UpSetR::fromList(dbl_list), 
              nintersects = NA, 
              sets = rev(names(dbl_list)),
              keep.order = TRUE,
              order.by = "freq")
Show plot


Now we can remove barcodes identified as cell doublets, and visualize the cell space before and after.

10.5 Doublets filtering

10.5.1 BEFORE

### Doublets viz (before removal)
umap_dbl_unfilt <- SC.helper::QnD_viz(
  sobj = sobj, 
  group_by = "doublet_viz",
  return_plot = TRUE)
Show plot