EBAII n1 2023 : SINGLE CELL ANALYSIS TRAINING
PREPROCESSING (II)
Technical and biological metrics, cell cyle , doublets and filtering

bastien.job@gustaveroussy.fr

Show TOC

1
2
3

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

## VIZ is not needed anymore
rm(VIZ)

Question1 : 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 !

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

Answer

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 sparcity 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 a 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 can discard these features

What are the actual dimensions of our count matrix ?

dim(sobj)

[1] 31053  4587

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

## Computing the 0s per feature
nFeat_zero <- sparseMatrixStats::rowCounts(
  x = Seurat::GetAssayData(
    object = sobj, 
    assay = 'RNA', 
    slot = 'counts'), 
  value = 0)
nFeat_nonzero <- ncol(sobj) - nFeat_zero
## Summary
summary(nFeat_nonzero)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0     242     166    4575 

## Identifying those with at least than 5 cells with expression
nFeat_keep <- nFeat_nonzero >= 5

## Restrict the Seurat object
sobj <- sobj[nFeat_keep,]
## Quantify the selected ones
table(nFeat_keep)

nFeat_keep
FALSE  TRUE 
18545 12508 

What are the Seurat object dimensions, now ?

dim(sobj)

[1] 12508  4587

We can visualize the cell space after this features filtering

EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  group_by = 'CC_Seurat_Phase')

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.

As a reminder, here are the metrics thresholds we set

## Feature_RNA
print(nfeature_range)

[1]  750 5975

## nCount_RNA
print(ncount_range)

[1]  1000 48875

## %mito
print(pc_mito_range)

[1] 0.00 0.05

## %ribo
print(pc_ribo_range)

[1] 0 1

## %stress
print(pc_stress_range)

[1] 0.00 0.06

Selecting “good” cells

## Applying QC metrics thresholds
bc_keep <- sobj$nFeature_RNA > nfeature_range[1] & 
  sobj$nFeature_RNA < nfeature_range[2] & 
  sobj$nCount_RNA > ncount_range[1] & 
  sobj$nCount_RNA < ncount_range[2] & 
  sobj$percent_mt > pc_mito_range[1] & 
  sobj$percent_mt < pc_mito_range[2] & 
  sobj$percent_mt > pc_ribo_range[1] & 
  sobj$percent_mt < pc_ribo_range[2] & 
  sobj$percent_st > pc_stress_range[1] & 
  sobj$percent_st < pc_stress_range[2]
table(bc_keep)

bc_keep
FALSE  TRUE 
  311  4276 

Applying the filter

sobj <- sobj[, bc_keep]
dim(sobj)

[1] 12508  4276

We can visualize the cell space since this cells filtering

EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  group_by = 'CC_Seurat_Phase')

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 Cell doublets

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

• scds (in its “hybrid” mode)
• scDblFinder

7.1 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)

7.2 scDblFinder

## Fix seed
set.seed(my_seed)
## Run scDblFinder
sobj$doublet_scDblFinder <- unname(
  Seurat::as.Seurat(
    scDblFinder::scDblFinder(
      Seurat::as.SingleCellExperiment(
        x = sobj, 
        assay = 'RNA')))$scDblFinder.class == 'doublet')

7.3 Merge results

We merge results to identify tool-specific and common doublets

sobj$doublet_union <- sobj$doublet_scds.hybrid | sobj$doublet_scDblFinder
sobj$doublet_viz <- 'singlet'
sobj$doublet_viz[sobj$doublet_union] <- 'both'
sobj$doublet_viz[sobj$doublet_scds.hybrid & !sobj$doublet_scDblFinder] <- 'scds'
sobj$doublet_viz[sobj$doublet_scDblFinder & !sobj$doublet_scds.hybrid] <- 'scDblFinder'
sobj$doublet_viz <- as.factor(sobj$doublet_viz)
table(sobj$doublet_viz)

       both scDblFinder        scds     singlet 
        142          28          68        4038 

7.4 Removal and visualization

Now we can remove and visualize the cell space with/without doublets

BEFORE

### Doublets viz (before removal)
EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  group_by = 'doublet_viz')

Let’s remove doublets :

sobj <- sobj[,!sobj$doublet_union]
dim(sobj)

[1] 12508  4038

AFTER removal

### Doublets viz (after removal)
EBAII.n1.SC.helper::QnD_viz(
  sobj = sobj, 
  group_by = 'doublet_viz')

8 Save

We can now save the results of our hard work !

But, how ?

?base::saveRDS

saveRDS(object = sobj, 
        file = './TD3A_01_Filtered.RDS', 
        compress = 'bzip2')

Rsession

utils::sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=fr_FR.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       

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] GenomicRanges_1.54.1        Matrix_1.6-1.1             
[3] SummarizedExperiment_1.30.2 SingleCellExperiment_1.22.0
[5] IRanges_2.36.0              S4Vectors_0.40.1           
[7] Biobase_2.62.0              BiocGenerics_0.48.0        

loaded via a namespace (and not attached):
  [1] GSEABase_1.64.0               progress_1.2.2               
  [3] nnet_7.3-19                   goftest_1.2-3                
  [5] DT_0.30                       Biostrings_2.70.1            
  [7] vctrs_0.6.4                   spatstat.random_3.2-1        
  [9] digest_0.6.33                 png_0.1-8                    
 [11] shape_1.4.6                   shinyBS_0.61.1               
 [13] registry_0.5-1                ggrepel_0.9.4                
 [15] deldir_1.0-9                  parallelly_1.36.0            
 [17] MASS_7.3-60                   reshape2_1.4.4               
 [19] httpuv_1.6.12                 foreach_1.5.2                
 [21] withr_2.5.2                   ggrastr_1.0.2                
 [23] xfun_0.41                     ellipsis_0.3.2               
 [25] survival_3.5-7                memoise_2.0.1                
 [27] ggbeeswarm_0.7.2              zoo_1.8-12                   
 [29] GlobalOptions_0.1.2           pbapply_1.7-2                
 [31] Formula_1.2-5                 prettyunits_1.2.0            
 [33] KEGGREST_1.42.0               promises_1.2.1               
 [35] httr_1.4.7                    restfulr_0.0.15              
 [37] globals_0.16.2                fitdistrplus_1.1-11          
 [39] rstudioapi_0.15.0             shinyAce_0.4.2               
 [41] miniUI_0.1.1.1                generics_0.1.3               
 [43] base64enc_0.1-3               curl_5.1.0                   
 [45] zlibbioc_1.48.0               ScaledMatrix_1.10.0          
 [47] polyclip_1.10-6               ca_0.71.1                    
 [49] GenomeInfoDbData_1.2.11       ExperimentHub_2.10.0         
 [51] SparseArray_1.2.0             RBGL_1.78.0                  
 [53] threejs_0.3.3                 interactiveDisplayBase_1.40.0
 [55] xtable_1.8-4                  stringr_1.5.0                
 [57] doParallel_1.0.17             evaluate_0.23                
 [59] S4Arrays_1.2.0                BiocFileCache_2.10.1         
 [61] hms_1.1.3                     irlba_2.3.5.1                
 [63] colorspace_2.1-0              filelock_1.0.2               
 [65] ROCR_1.0-11                   reticulate_1.34.0            
 [67] pcaExplorer_2.28.0            spatstat.data_3.0-3          
 [69] magrittr_2.0.3                lmtest_0.9-40                
 [71] Rgraphviz_2.46.0              later_1.3.1                  
 [73] viridis_0.6.4                 lattice_0.22-5               
 [75] spatstat.geom_3.2-7           NMF_0.26                     
 [77] future.apply_1.11.0           genefilter_1.84.0            
 [79] SparseM_1.81                  scattermore_1.2              
 [81] XML_3.99-0.14                 scuttle_1.12.0               
 [83] cowplot_1.1.1                 matrixStats_1.0.0            
 [85] RcppAnnoy_0.0.21              Hmisc_5.1-1                  
 [87] pillar_1.9.0                  nlme_3.1-163                 
 [89] iterators_1.0.14              gridBase_0.4-7               
 [91] compiler_4.3.1                beachmat_2.18.0              
 [93] stringi_1.7.12                Category_2.68.0              
 [95] TSP_1.2-4                     EBAII.n1.SC.helper_0.0.2     
 [97] tensor_1.5                    dendextend_1.17.1            
 [99] GenomicAlignments_1.36.0      plyr_1.8.9                   
[101] scater_1.28.0                 BiocIO_1.10.0                
[103] crayon_1.5.2                  abind_1.4-5                  
[105] SoupX_1.6.2                   locfit_1.5-9.8               
[107] sp_2.1-1                      bit_4.0.5                    
[109] dplyr_1.1.3                   codetools_0.2-19             
[111] BiocSingular_1.18.0           crosstalk_1.2.0              
[113] bslib_0.5.1                   GetoptLong_1.0.5             
[115] plotly_4.10.3                 mime_0.12                    
[117] splines_4.3.1                 circlize_0.4.15              
[119] Rcpp_1.0.11                   dbplyr_2.4.0                 
[121] sparseMatrixStats_1.14.0      knitr_1.45                   
[123] blob_1.2.4                    utf8_1.2.4                   
[125] clue_0.3-65                   BiocVersion_3.18.0           
[127] listenv_0.9.0                 checkmate_2.2.0              
[129] DelayedMatrixStats_1.24.0     tibble_3.2.1                 
[131] statmod_1.5.0                 pkgconfig_2.0.3              
[133] pheatmap_1.0.12               tools_4.3.1                  
[135] cachem_1.0.8                  RSQLite_2.3.2                
[137] viridisLite_0.4.2             DBI_1.1.3                    
[139] celldex_1.12.0                scDblFinder_1.14.0           
[141] fastmap_1.1.1                 rmarkdown_2.25               
[143] scales_1.2.1                  grid_4.3.1                   
[145] ica_1.0-3                     Seurat_4.4.0                 
[147] shinydashboard_0.7.2          Rsamtools_2.16.0             
[149] AnnotationHub_3.10.0          sass_0.4.7                   
[151] patchwork_1.1.3               BiocManager_1.30.22          
[153] dotCall64_1.1-0               graph_1.80.0                 
[155] RANN_2.6.1                    SingleR_2.4.0                
[157] rpart_4.1.21                  farver_2.1.1                 
[159] yaml_2.3.7                    AnnotationForge_1.44.0       
[161] MatrixGenerics_1.14.0         foreign_0.8-85               
[163] rtracklayer_1.60.1            cli_3.6.1                    
[165] purrr_1.0.2                   stats4_4.3.1                 
[167] webshot_0.5.5                 leiden_0.4.3                 
[169] lifecycle_1.0.3               uwot_0.1.16                  
[171] bluster_1.12.0                backports_1.4.1              
[173] BiocParallel_1.36.0           annotate_1.80.0              
[175] gtable_0.3.4                  rjson_0.2.21                 
[177] ggridges_0.5.4                progressr_0.14.0             
[179] pROC_1.18.4                   parallel_4.3.1               
[181] limma_3.58.1                  jsonlite_1.8.7               
[183] edgeR_4.0.1                   seriation_1.5.1              
[185] bitops_1.0-7                  ggplot2_3.4.4                
[187] xgboost_1.7.5.1               bit64_4.0.5                  
[189] assertthat_0.2.1              Rtsne_0.16                   
[191] spatstat.utils_3.0-4          BiocNeighbors_1.20.0         
[193] SeuratObject_5.0.0            heatmaply_1.5.0              
[195] highr_0.10                    jquerylib_0.1.4              
[197] metapod_1.10.0                dqrng_0.3.1                  
[199] scds_1.16.0                   lazyeval_0.2.2               
[201] shiny_1.7.5.1                 htmltools_0.5.6.1            
[203] GO.db_3.18.0                  sctransform_0.4.1            
[205] rappdirs_0.3.3                glue_1.6.2                   
[207] spam_2.10-0                   XVector_0.42.0               
[209] RCurl_1.98-1.12               scran_1.30.0                 
[211] gridExtra_2.3                 igraph_1.5.1                 
[213] R6_2.5.1                      tidyr_1.3.0                  
[215] DESeq2_1.42.0                 labeling_0.4.3               
[217] cluster_2.1.4                 rngtools_1.5.2               
[219] GenomeInfoDb_1.38.0           DelayedArray_0.28.0          
[221] tidyselect_1.2.0              vipor_0.4.5                  
[223] htmlTable_2.4.2               GOstats_2.68.0               
[225] xml2_1.3.5                    AnnotationDbi_1.64.0         
[227] future_1.33.0                 rsvd_1.0.5                   
[229] munsell_0.5.0                 KernSmooth_2.23-22           
[231] topGO_2.54.0                  data.table_1.14.8            
[233] htmlwidgets_1.6.2             ComplexHeatmap_2.18.0        
[235] RColorBrewer_1.1-3            biomaRt_2.58.0               
[237] rlang_1.1.1                   spatstat.sparse_3.0-3        
[239] spatstat.explore_3.2-5        fansi_1.0.5                  
[241] beeswarm_0.4.0