EBAII n1 2022 : Single Cell Analysis Training PREPROCESSING - QC, Filtering and normalization

Show TOC

1
2
3

This file describes the different steps to perform the training course for the EBAII n1 2022 Single Cell analysis, covering the preprocessing steps : * Raw matrix loading * Empty droplets filtering * Ambient RNA contamination estimation * QC metrics (%mito, %ribo, %stress, cell cycle status) * Barcodes filtering * Normalization

1 CHEAT SHEET

Here is a list of different guidelines you should refer to, as they correspond to actions you would perform frequently.

1.1 JupyterHub

1.1.1 Connecting, and opening a session

• Connect to the IFB JupyterHub using your credentials :
• Select the reservation dedicated to our training session : form_2022_32
• Reservation :
  • For a Rstudio session, request 2 CPUs and 15 GB of memory (let the other options to their default value)
  • For a shell session, request the default minimum (1 CPU, 1 GB)

1.1.2 Inside a session

• To open an app, click on its icon.
• To get multiple apps open simultaneously
  • Click on [File > New Launcher]
  • or just click on the big blue [ + ] button

1.1.3 Closing a session

• Close your browser/tab, won’t close your session. It will still be running in background despite being invisible : this is a normal behavior, so that you can get back later to your current active session.
• To actually close and finish your session, you have to :
  • Select from the JupyterHub menu [File > Hub Control Panel]
  • Click on the red button [Stop my server]
  • This will release the resources used by your server

1.2 singleCellTK

1.2.1 Starting scTK

• In the Rstudio console :

scTKcomp::scTK_launch()

• Let the interface start, it should take a few dozen of seconds
• singleCellTK opens a new tab in your favorite internet browser, opened by default on the dataset loading page.

1.2.2 Loading a dataset (count matrix)

• [Data > Import Single Cell data]
  • For this course, the data source will always be a SingleCellExperiment as RDS file
  • For each dataset you want to load
    • [Browse] to select the ‘*.rds’ file containing the SingleCellExperiment you want to load
    • Wait for loading bar to be fullfilled (should be immediate for early steps ; may take a few dozen seconds for later steps, when the SCE object contains multiple matrices)
    • Then, [Add to sample list]
  • If you are importing a new dataset but another one is already loaded, please remember to Overwrite existing SCE object (radio button) before importing !
  • Once you added all the datasets you need, [Import]
  • Let the features for display as Rownames (which, for our course, are symbols) and [Set]
  • You may switch to the next analysis steps

1.2.3 Exporting a dataset

For this course, we will systematically export objects from singleCellTK in the SCE (SingleCellExperiment) format, saved as a RDS archive (a compressed format which contained a single, unnamed R object)

• [Data > Export Single Cell data]
  • Select the output directory (your R working directory should already be well positioned)
  • Let the export type as RDS
  • You can use whatever file name, but it is heavily recommended to follow this course’s suggestions.

1.2.4 Deleting cell annotations

Sometimes, it will be recommended to remove some entries in the working SCE object, mainly to : * Avoid being confused among multiple matrices in the object, discarding those that won’t be required anymore at further steps * Slim down the object for faster import / export.

• [Data > Delete Single Cell data]
  • Just check the data to remove.
  • Confirm with [Delete]
  • Please remember that any suppression is definitive for the currently loaded object.

1.2.5 Stopping singleCellTK

• Close your scTK browser tab
• In Rstudio :
  • Clic in the console and press [Esc].
  • or clic on the red [Stop] button.

1.3 scTKcomp

This package hosts a collection of functions built around singleCellTK to :

• Fill some gaps existing in its course
  • Create a SingleCellExperiment object from multiple sources of data (Cell Ranger, UMI-tools, BUStools, Alevin) (scTK_load)
  • Filter empty droplets (scTK_edf)
  • Evaluate cell cycle phase status / scores, using Seurat (scTK_cc_seurat)
  • Regress of cell covariates on any type of matrix, using Seurat too (scTK_regress_covariate)
• Bypass some bugs
  • Perform a UMAP reduction (still with Seurat) on a reduced dimension object which is not a PCA/ICA, without creating a SCE object that can’t be properly imported anymore (scTK_SeuratUMAP)
• Add new functionalities
  • Quickly describe the content of a SingleCellExperiment object (scTK_descriptor)
  • Perform a quick and dirty UMAP for early stage matrices (for early visualization) (scTK_QnDuMAP)
  • Use of “sensitive genes” instead of “highly variable genes” for the feature-based dimension reduction, thanks to scSensitiveGeneDefine (scTK_scSG)
  • Assess the weight of cell covariates on a matrix, with a heatmap visualization (scTK_assess_covar)
  • Perform trajectory analysis with TinGa (scTK_TinGa)
• Misc
  • Recompress a RDS with a better compression ratio than default (ie, from gzip to bzip2 compression) (scTK_recompress)

We are now done with the blah blah, let’s practice for real !

2 WARMUP

2.1 Preamble

2.1.1 The analysis framework/interface : singleCellTK

singleCellTK : * Allows to use some of the R tools we bioinformaticians use * Adds a graphical interface * Allows visualization of some of the results * Obviously there are some limits in what can be done, compared to the command line * Relies on a format to store and transform the data :

2.1.2 The SingleCellExperiment object

We will now see together the scructure of the R data object you will manipulate for the next fewdays : the SingleCellExperiment object.

2.2 Getting the source data

• Create an IFB JupyterHub session for Rstudio
• Open a terminal
• Create a new “TD” directory in your project folder

MYPROJECTNAME=my_project_name
mkdir -p /shared/projects/${MYPROJECTNAME}/TD/DATA

• Copy the source data

cp -r /shared/projects/form_2022_32/SingleCellRNASeq/TD/DATA/START_OBJECTS/* /shared/projects/${MYPROJECTNAME}/TD/DATA/

• Copy the required resource files

cp -r /shared/projects/form_2022_32/SingleCellRNASeq/TD/RESOURCES /shared/projects/${MYPROJECTNAME}/TD/

• Exit your terminal

exit

2.3 Starting Rstudio

• Just click on the Rstudio logo
• Optionally : load the scTKcomp package

library(scTKcomp)

3 FROM COUNTS TO CELLS

3.1 Purpose

Our typical input data (once the reads mapping has been done, and the count matrix generated) is a features (row) by barcodes (columns) matrix.

Our aim here is to go from these experimental barcodes to actual cells (or, at least, most probable cells)

We will now : * Filter barcodes that should correspond to empty droplets * Perform QC, evaluating : * The amount of UMIs and number of expressed features per barcode * Which barcodes correspond to cell multiplets (mainly doublets) * The level of expression per barcode of : * Mitochondrial genes (proxy for over-lysis / cell mortality) * Ribosomal protein genes (proxy for metabolic level status of cells) * Dissociation / mechanical stress response genes (from Macchiado PhD report) * We will then filter out barcodes according to these generated QC metrics, using a series of criteria.

3.2 Generating an SCE object from raw data

For this very first step, we will use a rather small public dataset provided by 10X Genomics, which consists into an assay containing ~1000 human PBMCs (Peripheral Blood Mononucleus Cells). The source data for this sample consists in a raw count matrix in the 10X Cell Ranger format, and are located in the [TD/DATA/PBMC1K] folder

We will first generate a SingleCellExperiment object from the raw Cell Ranger matrix :

• Position your R working directory t o where the data are :

setwd('/shared/projects/<yourprojectname>/TD/DATA/PBMC1K')

• Call the help for the scTK_load() function

?scTKcomp::scTK_load

• Run it on the 10X PBMC 1K dataset

sname <- 'PBMC1K'
scTKcomp::scTK_load(
  data_path = '.',
  sample_name = sname,
  exp_name = sname
)

• This should have created the [PBMC1K_00a_Raw.rds] file that contains a SingleCellExperiment
• Read the console outputs while they pop
  • Q: Is the observed sparsity level alarming ?

3.3 Filtering empty droplets

We will now filter out all these empty barcodes/droplets from this huge matrix :

• Call the help for the scTK_edf() function

?scTKcomp::scTK_edf

• Run it on the freshly created [PBMC1K_00a_Raw.rds] SCE object

scTKcomp::scTK_edf(in_rds = './PBMC1K_00a_Raw.rds')

• Read the console outputs while they pop
  • Q1 : How many droplets were retained ?
  • Q2 : Compare this value from emtyDrops with the one obtained with Cell Ranger (in the HTML report). Are they concordant ?
• Observe the generated plots (kneeplot, saturation plot)
  • Q3 : What could you say from the kneeplot ?
• We will now take a look inside the guts of the generated SCE object :
  • Call the help for the function scTK_descriptor()
  • • ?scTKcomp::scTK_descriptor
  • Use this function on the newly created object [PBMC1K_01a_EDf.rds]
• Observe the console output (and note the sparsity level)

3.4 Estimating contamination

We will now work with the ‘TDCT’ sample from Paiva et al. It consists in a cell dissociation of a grafted murine thymus from a juvenile mouse to another juvenile mouse, that serves as a control in their study.

• Start Rstudio (if not already started)
• Get to the data directory corresponding to the TDCT sample

setwd('/shared/projects/<myprojectname>/TD/DATA/TDCT')

• Start scTK

3.4.1 Importing data

• Import the [TDCT_01a_EDf.rds] SCE object
• Get to [QC & Filtering > QC]
• Run SoupX with its default parameters
• Observe the generated uMAPs (and example soup markers)
• Save the SCE as [TDCT_01X_SoupX.rds]
• Stop scTK
• Use scTK_descriptor on [TDCT_01X_SoupX.rds]

scTKcomp::scTK_descriptor('TDCT_01X_SoupX.rds')

• Observe the console output
  • Q1 : What did SoupX actually do ?
• (You may now delete [TDCT_01X_SoupX.rds], we won’t use it anymore)

Just for the show : Run scTK_descriptor on a post-SoupX RDS on an other sample : N705, from Bacon et al, 2018) :

scTKcomp::scTK_descriptor('/shared/projects/form_2022_32/SingleCellRNASeq/TD/DATA/ALL_STEPS_OBJECTS/N705/N705_01X_SoupX.rds')

3.5 QC Metrics : Genesets expression & doublets

We will now generate some QC metrics for our cells

• Start scTK
• Set your working directory on the sample data folder

setwd('/shared/projects/<yourprojectname>/TD/DATA/TDCT')

• Import [TDCT_01a_EDf.rds] object again
• Import the gene lists required for some QCs :
  • (Mus musculus mitochondrial genes are bundled with singleCellTK)
  • Ribosomal proteins geneset :
    • [Data > Import Gene Sets]
      • Let the list selection to ‘Upload a GMT file’
      • [Browse] to select the ‘../../RESOURCES/GENELISTS/mus_musculus_ribo_symbols_20221105.gmt’
      • Set ‘Ribo’ as Collection Name
      • [Upload].
  • Dissociation / mechanical stress geneset :
    • [Data > Import Gene Sets]
      • In the same way, let the list selection to ‘Upload a GMT file’
      • [Browse] to select the ‘../../RESOURCES/GENELISTS/mus_musculus_stress_symbols_20221107.gmt’
      • Set ‘Stress’ as Collection Name
      • [Upload].
• Perform doublets identification with method :
  • Team A : scDblFinder
  • Team B : scds (cxds bcds hybrid)
  • Team C : both
  • Team D : DoubletFinder
• Observe generated uMAP
• Perform QC (and observe QC plots each time) with
  • Mito
  • Ribo
  • Stress
• Get to the [Cell Viewer] to observe QC metrics on the [QC_UMAP]
• Delete all metadata but the first 5 annotations, *_percent, percent.top.* and *_call
• Save the SCE object as [TDCT_01b_QC.rds]
• Stop scTK

3.6 Estimate the cell cycle phase

We will now estimate the cell cyle phase status of our cells thanks to Seurat

• Call the help for the scTK_cc_seurat() function

?scTKcomp::scTK_cc_seurat

• This function will also need another RDS object (available in the [/shared/projects//TD/RESOURCES/GENELISTS] direcory) that contains 2 genesets :
  • S phase signature
  • G2M phase signature
• Run it on the freshly created [TDCT_01b_QC.rds] SCE object.

scTKcomp::scTK_cc_seurat(
  in_rds = 'TDCT_01b_QC.rds',
  assay = 'counts',
  cc_seurat_file = '/shared/projects/<yourprojectname>/TD/RESOURCES/GENELISTS/mus_musculus_Seurat_cc.genes_20191031.rds')

• Read the contingency table for estimated phases in the console output
• Start scTK
• Set your working directory on the sample data folder

setwd('/shared/projects/<yourprojectname>/TD/DATA/TDCT')

• Import the freshly created [TDCT_02a_counts_CC.seurat.rds] SCE object
• in the [Cell Viewer], on the QC_UMAP, observe cell annotations :
  • CC phase
  • SmG2M score, with a split by CC Phase

3.7 Filter cells and features

• At the barcode level : Apply these filtering criteria (step by step, to measure their impact):
  • nFeature > 199
  • nFeature < 2999
  • nCount >999
  • Mito_percent < 5 (%)
  • doublets_call : Singlet
• At the feature level : keep features with :
  • 1 count in at least 5 cells
• Delete unneeded QC_UMAP reduction
• Save object as [TDCT_02b_QCf.rds]

4 NORMALIZATION

4.1 Merge datasets

The study from which came the TDCT sample had another sample : TD3A. This second sample consists into another graft of a juvenile mouse thymus to another juvenile one, except that in this case the host is mutated so that it cannot produce the expected T cell progenitors from its spine bone marrow.

In normal circumstances, when there is a short period during which the spine bone marrow is unable to produce such progenitors, it is already known that a mechanism (called autonomy) activates, such that some unknown immature cell type from the thymus behaves as a progenitor. The aim of this study was to identify which cell type in the T cell maturation stages has this ability.

It is also important to note that such mechanism, when lasting to long, is very prone to start leukemia…

We will know merge the two samples to analyze them, in a comparison purpose.

• Stop scTK
• Set your working directory on the sample data folder

setwd('/shared/projects/<yourprojectname>/TD/DATA/INTEGRATED')

• Start scTK
• Import [../TDCT/TDCT_02b_QCf.rds] and [../TD3A/TD3A_02b_QCf.rds].
  • WARNING : This time you import the 2 objects at the same time. Be careful to add them both to the sample list before importing.
• Delete “rownames” (bypasses a bug for QC plots)
  • WARNING : Do not skip this step or you’ll be in trouble soon !
• Perform basic QCs and observe boxplots
  • If you are lost in an infinite waiting loop, you did not read the last warning :D
  • Q1 : How similar are these two samples when considering :
    • Their globally measured expression level ?
    • The expression level of their top 50 expressed genes ?
• Delete unneeded barcode annotations : sum, detected, total
• Save the SCE object as [TDCT.TD31_02b_QCf.rds]
• Stop scTK
• Run scTK_descriptor (only describing assays) on these 3 SCE objects :
  • [TDCT_02b_QCf.RDS] (Control)
  • [TD3A_02b_QCf.RDS] (Autonomy)
  • [TDCT.TD31_02b_QCf.RDS] (Merged)
    • Q2 : Do things add up (features / cells) ?

4.2 Normalization

• Start scTK
• Set your working directory on the sample data folder

setwd('/shared/projects/<yourprojectname>/TD/DATA/INTEGRATED')

• Import again the [TDCT.TD3A_02b_QCf.rds] SCE object
• Get to [Normalization & Batch Correction > Normalization]
• Use :
  • Method : Seurat - LogNormalize
  • Input assay : counts
  • Output assay name = ‘SLNst’
  • Use scaling
  • Use trimming, let default limits
• Save the SCE object as [TDCT.TD3A_02c_SLNst.rds]
  • This should take more time than you are used to from earlier steps. This is due to the fact that while the raw ‘counts’ matrix consists in a very sparse, thus very compressed matrix, the normalized ‘SLNst’ matrix contains numeric values, which compress with a mediocre ratio
• Stop scTK
• Describe the freshly saved [TDCT.TD3A_02c_SLNst.rds]

scTKcomp::scTK_descriptor('TDCT.TD3A_02c_SLNst.rds')

• Q1 : What could you say for our new assay ?