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

  • Computing and evaluating QC metrics
    • Based on technical characteristics of the count matrix
    • Based on biological characteristics
  • Estimation of the cells cycle phase
  • Filtering of features and cells
  • Identification and removal of cell doublets

1 PREREQUISITES

1.1 The training data set

  • Data from the Paiva et al. publication.
  • The study concerns thymus autonomy :
    • The thymus is an “organ of passage”, critical in its function to the adaptative immune system for the maturation of T cell lymphocytes.
    • This maturation involves two main steps, performed thanks to macrophages :
      • Positive selection : keeping cells that successfully develop react appropriately with MHC immune receptors of the body
      • Negative selection : keeping cells that do not react against natural proteins of the body.
    • Thymus autonomy is a natural mechanism that allows to create T cells in the thymus by differenciation and cell competition, even when normal progenitors from the bone marrow are lacking, in critical conditions.
    • This mechanism is known in its effects, but the cells involved in are not.
    • This study is of importance in the health field, as this mechanism relies on a temporary loss of control of the cell normal functions.
    • The consequence is that if thymus is in autonomy for too long (few weeks), this is a prelude for leukaemia !
  • Organism is : mus musculus
  • Individuals are : mice in development, grafted
  • The design corresponds to two conditions (Test / control)
    • Control : thymus from wild type newborn mouse transplanted into wild type juvenile mouse. In this control case, donor T-cells progenitors (DN3) were replaced by host cells 3 weeks after transplantation.
    • Test : thymus from wild type newborn mouse transplanted in KO Rag-/- type juvenile mouse (the KO partially impairs their ability to produce T-cell progenitors in normal amounts). In this test case, donor T-cells progenitors (DN3) were replaced by host cells 9 weeks after transplantation, showing that the donor DN3 cells outlived their normal lifespan by ~6 weeks.
  • You will mainly work on the KO sample (‘TD3A’)
  • Input data consists in a count matrix, as a gzip’d tabular text file, that contains everything needed to create a basic Seurat object :
    • The expressions counts
    • The feature names (here, gene symbols)
    • The barcode names
  • This matrix has already been filtered for empty droplets.

1.2 Copy data and resources

  • Using the Jupyter cheat sheet, connect to JupyterHub and create a session using the resource requirements for an Rstudio session.
  • In the [Launcher] panel, click on [Terminal] in the [Other] section :


  • You should be in your home directory by default. Check it by typing :
pwd
  • Go to your project DATA directory
cd /shared/projects/<your_project>/TD/DATA
pwd
  • Copy the required input file :
cp -r /shared/projects/2325_ebaii/SingleCell/TD_DATA/DATA_START/TD3A .
  • Get back to your project directory :
cd ..
# or : cd /shared/projects/<your_project>/TD
  • Your directory should look like this (mine is called “golf”)


  • You can check by yourself :
tree -sh $PWD
  • We need to create a new directory in which we will copy the other resources needed for this session. Here we will try a variation : copy without moving !
## Here, we diretly copy the whole RESOURCES directory from the origin (2325_ebaii project) to its destination (your project)
cp -r /shared/projects/2325_ebaii/SingleCell/TD_DATA/RESOURCES /shared/projects/<your_project>/TD/
  • Your directory should look like this


  • You can check by yourself :
tree -sh $PWD
  • You can now create your output directory
mkdir ./RESULTS
  • You can now close the terminal.

2 Load raw data

Start a Rstudio session


Set your working directory in your TD output directory (golf is mine!)

setwd('/shared/projects/golf/TD/RESULTS')

Setting parameters

## Data directory
input_matrix <- '../DATA/TD3A/GSM4861194_gex_2_raw_gene_expression.tsv.gz'

## Sample name
samplename <- 'TD3A'

## Seed for the RNG
my_seed <- 1337

We can load the matrix

## R can read compressed text file without hassle
scmat <- as.matrix(
  read.table(
    file = input_matrix, 
    header = TRUE, 
    sep = '\t'))
## Displaying its size in-memory (this is a basic matrix)
format(object.size(scmat), units = "auto")
[1] "545.7 Mb"

A quick look at its content

## Displaying its structure
str(scmat)
 int [1:31053, 1:4587] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:31053] "Xkr4" "Gm1992" "Gm37381" "Rp1" ...
  ..$ : chr [1:4587] "AAACCTGAGACGCTTT.1" "AAACCTGAGGCATTGG.1" "AAACCTGGTCAACATC.1" "AAACCTGTCGAGGTAG.1" ...


Wait a second …
We have a raw count matrix that’s already been filtered for empty droplets. What may we do with it before creating a Seurat object from it ?

## Run SoupX !
soup_frac <- EBAII.n1.SC.helper::SoupX_auto(scmat_filt = scmat)

Soup-contributing features (Top 20) :


|        |       est| counts|
|:-------|---------:|------:|
|Actb    | 0.0080580| 132061|
|Gm42418 | 0.0069444| 113810|
|Eef1a1  | 0.0055896|  91607|
|Malat1  | 0.0051948|  85137|
|Ppia    | 0.0049767|  81562|
|Tmsb10  | 0.0046515|  76232|
|mt-Co1  | 0.0044872|  73540|
|Ptma    | 0.0044203|  72443|
|mt-Atp8 | 0.0040750|  66784|
|Rplp0   | 0.0038984|  63890|
|Uba52   | 0.0038012|  62297|
|Gm10076 | 0.0037717|  61814|
|Rps24   | 0.0037488|  61438|
|mt-Co2  | 0.0036859|  60407|
|Pfn1    | 0.0035459|  58113|
|Fau     | 0.0034193|  56038|
|Ubb     | 0.0033782|  55364|
|Rps16   | 0.0032381|  53068|
|Rpl13   | 0.0031730|  52002|
|Rps3a1  | 0.0031379|  51426|

cat('Soup fraction : ', soup_frac)
Soup fraction :  0.01

We can create our Seurat object

## Create Seurat object
sobj <- Seurat::CreateSeuratObject(
  counts = scmat, 
  project = samplename, 
  assay = 'RNA')
## Remove the matrix to free some RAM
rm(scmat)
## Displaying sobj size in-memory (matrices are converted to dgCMatrix)
format(object.size(sobj), units = "auto")
[1] "179.4 Mb"

Let’s describe it

## Describe the object
EBAII.n1.SC.helper::seurat4_descriptor(sobj = sobj)
OBJECT VERSION :    5.0.0 
PROJECT :   [TD3A] 

[ASSAYS]
   ASSAY 1 :    [RNA] [ACTIVE] 
      SLOT 1 :  [counts]    Dims:[31053 x 4587]  Range:[0.00-3103.00] 
         Counts :   16388810 
         Sparsity : 94.72482% 
      SLOT 2 :  [data]  Dims:[31053 x 4587]  Range:[0.00-3103.00] 
         Sparsity : 94.72482% 
      SLOT 3 :  [scale.data]    Dims:[0 x 0] 

[DIMREDS]

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

 nCount_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    502    1986    2397    3573    3134   48875 

 nFeature_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     22    1291    1476    1638    1733    5975

3 Count metrics

Here we will focus on the barcode/cell-level count metrics :

  • nCount_RNA : Total number of counts in each cell/barcode
  • nFeature_RNA : Number of features/gens with at least 1 count in each cell/barcode (ie : number of expressed genesin cell)

Both were automatically computed by Seurat when creating the object (but you already observed it, isn’t it ?)

For graphical purpose, we will convert count metrics (nCount_RNA) to log10

## Add l10counts ====
sobj$log10_nCount_RNA <- log10(sobj$nCount_RNA+1)

Question : Why this “+1” ?

## . Because most counts are 0, and log10(0) = -Inf

3.1 Visualization and thresholds

We can visualize these metrics’ distribution as violins …

Seurat::VlnPlot(
  object = sobj, 
  features = c('nFeature_RNA', 
               'nCount_RNA', 
               'log10_nCount_RNA'))

Question : Please describe these distributions (modes, tailing, outliers, …)

… but for my own part, I find them hard to read. I prefer histograms.



nCount_RNA 

hist(x = sobj$nCount_RNA, breaks = 1000)

Still hard to read… let’s focus on the left part

hist(x = sobj$nCount_RNA, 
     breaks = 1000, 
     xlim = c(0,10000))

Question : Which thresholds would you pick to select cells of interest, based on nCount_RNA ?

## . For the left part of the distribution (very low 
##   counts), I would take 1000 as a lower bound
## . For the right part ... this is highly debatable 
##   due to the high spread of this right tailing. 
##   Removing it infers the risk of discarding somme 
##   cell subtypes with different characteristics than 
##   the mode of the distribution. Keeping them may 
##   just pollute the dataset... I'm feeling 
##   optimistic, I would take all values (so the upper 
##   threshold would be the distribution max), but will 
##   keep in mind that I may have to adapt my anaylsis 
##   to this variation in global expression levels.

Let’s set this nCount_RNA selection range and display it

## Setting nCount conservation range
ncount_range <- c(1000, max(sobj$nCount_RNA))
## Replot histogram with cutoffs
hist(x = sobj$nCount_RNA, 
     breaks = 1000, 
     xlim = c(0,10000))
abline(v = ncount_range, lty = 3, col = "red")

nFeature_RNA 

hist(x = sobj$nFeature_RNA, breaks = 1000)

Question : Which thresholds would you pick to select cells of interest, based on nFeature_RNA ?

## . For the left part of the distribution (very low 
##   counts), I would take 750 as a lower bound.
## . For the right part ... for the same reason, 
##   I would take all values (so the distribution 
##   max)

Let’s set this nFeature_RNA selection range and display it

## Setting nFeature conservation range
nfeature_range <- c(750, max(sobj$nFeature_RNA))
## Replot histogram with cutoffs
hist(x = sobj$nFeature_RNA, breaks = 1000)
abline(v = nfeature_range, lty = 3, col = "red")



However,
we won’t apply these filtering criteria right now, because we have other metrics to consider.

4 Technical and biological metrics

Here we will build 3 new metrics that correspond to the fraction of expression (in counts) of three categories of gene signatures among the global expression. These are

  • Mitochondrial genes : A good proxy to display cells in bad condition (high level of expression for those genes indicates that the cell outer membrane lysis went to far and also began for internal organites).
  • Genes coding for ribosomal proteins (and not ribosomal genes as often used as a bad shortcut) : High expression of these witnesses of higher metabolic rates for concerned cells, thus it may be useful to regroup cells of a same cell type but with different metabolic states (linked to cell division, by example). However, different cell types may naturally have default different metabolic rates by default …
  • Mechanical stress genes : In droplet-based protocols, during the separation, cells flow through a microfluidic circuit, during which they pass through hard bottlenecks compared to their diameter. For very fragile cells, this may have an impact in their gene expression. Actually, fragile cells are most often completely lost/destroyed in the process, so the interest of this metric is often very limited … but it may be in some cases !

4.1 Compute metrics

We will need a function that eases the computational step for you :

## Reading the function help page
?EBAII.n1.SC.helper::count_rate_metric



%MITO 

## Load the genelist
mito_symbols <- readRDS(file = '../RESOURCES/GENELISTS/mus_musculus_mito_symbols_20191015.rds')
## Compute the metric
sobj$percent_mt <- EBAII.n1.SC.helper::count_rate_metric(
  sobj = sobj,
  features = mito_symbols,
  assay = 'RNA')
## Summarize the computed values
summary(object = sobj$percent_mt)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.01985 0.02489 0.02997 0.03125 0.97250
%RIBO 

## Load the genelist
ribo_symbols <- readRDS(file = '../RESOURCES/GENELISTS/mus_musculus_cribo_symbols_20191015.rds')
## Compute the metric
sobj$percent_rb <- EBAII.n1.SC.helper::count_rate_metric(
  sobj = sobj,
  features = ribo_symbols,
  assay = 'RNA')
## Summarize the computed values
summary(object = sobj$percent_rb)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.07884 0.09577 0.10880 0.12167 0.42010
%STRESS 

## Load the genelist
stress_symbols <- readRDS(file = '../RESOURCES/GENELISTS/mus_musculus_stress_symbols_20200224.rds')
## Compute the metric
sobj$percent_st <- EBAII.n1.SC.helper::count_rate_metric(
  sobj = sobj,
  features = stress_symbols,
  assay = 'RNA')
## Summarize the computed values
summary(object = sobj$percent_st)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.02946 0.03338 0.03368 0.03749 0.10198

4.2 Visualization and thresholds

We can then visualize these metrics as violins (as well as older ones)

Seurat::VlnPlot(object = sobj, features = c(
  'nFeature_RNA', 'nCount_RNA', 
  'log10_nCount_RNA', 'percent_mt', 
  'percent_rb', 'percent_st'), 
  ncol = 3)

Did you notice I preferred histograms ? :)



%MITO 

hist(x = sobj$percent_mt, breaks = 1000)

Hard to read, let’s slash it to focus on lower values

## Restrict to first 10%
hist(x = sobj$percent_mt, 
     breaks = 1000, 
     xlim = c(0,.1))

Let’s observe the %mito distribution in the data space

## Performing the "quick viz" and saving 
## it (to speed up next plots)
VIZ <- EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  features = 'percent_mt',
  return_object = TRUE)

Let’s set this %mito selection range to [0, 5%] and display it

## Setting %mito conservation range
pc_mito_range <- c(0, .05)
# pc_mito_range <- c(0, .04)
## Replot histogram with cutoffs
hist(x = sobj$percent_mt, 
     breaks = 1000, 
     xlim = c(0,.1))
abline(v = pc_mito_range, 
       lty = 3, 
       col = "red")

%RIBO 

hist(x = sobj$percent_rb, breaks = 1000)

Let’s observe the %ribo distribution in the data space

## Performing the "quick viz" faster with 
## the "VIZ" Seurat object
EBAII.n1.SC.helper::QnD_viz(
  sobj = VIZ, 
  slot = NULL, dimred = 'umap',
  features = 'percent_rb')

We don’t have any clue to consider thise %ribo variations as biological or technical artefacts… So we will keep all values : [0, 1]

## Setting %ribo conservation range
pc_ribo_range <- c(0, 1)
%STRESS 

hist(x = sobj$percent_st, breaks = 1000)

Let’s observe the %stress distribution in the data space

## Performing the "quick viz" faster with 
## the "VIZ" Seurat object
EBAII.n1.SC.helper::QnD_viz(
  sobj = VIZ, 
  slot = NULL, dimred = 'umap',
  features = 'percent_st')

Let’s set this %stress selection range to [0, 6%] and display it

## Setting %stress conservation range
pc_stress_range <- c(0, .06)
## Replot histogram with cutoffs
hist(x = sobj$percent_st, 
     breaks = 1000, 
     xlim = c(0,.1))
abline(v = pc_stress_range, lty = 3, 
       col = "red")

We can delete the VIZ object

rm(VIZ)



But we’re not done yet !

Single cell data from droplet-based technologies can be biased in so many ways…

Evaluating their quality, increasing the knowledge of biases affecting the values is mandatory to perform an analysis of quality !

5 Cell cycle scores

As we are analysis cells independently, they may each be in a different phase of their cell cycle. Interestingly, the effect of this cell cycle step on the cell genes expression is strong, and can bias the data. 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.

Some of these require the evaluation of genes which expression is expected to be null : that’s the reason why we did not filter out for unexpressed genes till now.

5.1 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 :

?EBAII.n1.SC.helper::CC_Seurat

Run it !

## Load the cell cycle reference genes lists
seu.cc <- readRDS(file = '../RESOURCES/GENELISTS/mus_musculus_Seurat_cc.genes_20191031.rds')
## Perform the estimation
sobj <- EBAII.n1.SC.helper::CC_Seurat(
  sobj = sobj, assay = 'RNA', 
  seurat_cc_genes = seu.cc, SmG2M = TRUE, 
  nbin = 20, my_seed = my_seed)

Description of the object to see the data added

EBAII.n1.SC.helper::seurat4_descriptor(
  sobj = sobj,
  describe = 'coldata'
)
OBJECT VERSION :    5.0.0 
PROJECT :   [TD3A] 

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

 nCount_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    502    1986    2397    3573    3134   48875 

 nFeature_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     22    1291    1476    1638    1733    5975 

 log10_nCount_RNA 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.702   3.298   3.380   3.429   3.496   4.689 

 percent_mt 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.01985 0.02489 0.02997 0.03125 0.97250 

 percent_rb 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.07884 0.09577 0.10880 0.12167 0.42010 

 percent_st 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.02946 0.03338 0.03368 0.03749 0.10198 

 CC_Seurat_S.Score 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-2.6505 -0.3444 -0.2557 -0.2278 -0.1815  4.3548 

 CC_Seurat_G2M.Score 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-2.04357 -0.21453 -0.15637 -0.02719 -0.10224 12.76613 
CC_Seurat_Phase    Freq
----------------  -----
G1                 4259
G2M                 209
S                   119
NA                    0

 CC_Seurat_SmG2M.Score 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-13.3524  -0.1735  -0.1036  -0.2006  -0.0439   3.2976

5.2 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'))

But it’s not that easy to interpret… Let’s plot it in the cell space

  • For the S minus G2M 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 <- EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  features = 'CC_Seurat_SmG2M.Score', 
  return_object = TRUE)

  • For the 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
EBAII.n1.SC.helper::QnD_viz(sobj = VIZ, slot = NULL, dimred = 'umap', 
                            group_by = 'CC_Seurat_Phase')