Peripheral Blood
Mononuclear Cells (PBMC)
R analytical pipeline Tutorial
• A peripheral blood mononuclear cell (PBMC) is any peripheral blood
cell having a round nucleus.
• These cells consist of lymphocytes (T cells, B cells, NK cells) and
monocytes, whereas erythrocytes and platelets have no nuclei, and
granulocytes (neutrophils, basophils, and eosinophils) have multi-
lobed nuclei.
• In humans, lymphocytes make up the majority of the PBMC
population, followed by monocytes, and only a small percentage of
dendritic cells
Reading in the data
• [Link] <- Read10X([Link] =
"C:/Users/rahul/Downloads/filtered_gene_bc_matrices/hg19")
# Initialize the Seurat object with the raw
(non-normalized data).
• pbmc <- CreateSeuratObject(counts = [Link], project = "pbmc3k",
[Link] = 3, [Link] = 200)
• pbmc is the Seurat object
• ## An object of class Seurat
• ## 13714 features across 2700 samples within 1 assay
• ## Active assay: RNA (13714 features, 0 variable features)
• ## 1 layer present: counts
contents
Standard pre-processing workflow
• This part encompass the standard pre-processing workflow for
scRNA-seq data in Seurat.
• These represent the selection and filtration of cells based on QC
metrics
• Data normalization and scaling
• Detection of highly variable features.
• QC and selecting cells for further analysis
Seurat allows you to easily explore QC metrics and
filter cells based on any user-defined criteria
• A few QC metrics commonly used by the community include
• The number of unique genes detected in each cell.
• Low-quality cells or empty droplets will often have very few genes
• Cell doublets or multiplets may exhibit an aberrantly high gene count
• Similarly, the total number of molecules detected within a cell (correlates strongly
with unique genes)
• The percentage of reads that map to the mitochondrial genome
• Low-quality / dying cells often exhibit extensive mitochondrial contamination
• We calculate mitochondrial QC metrics with the PercentageFeatureSet() function,
which calculates the percentage of counts originating from a set of features
• We use the set of all genes starting with MT- as a set of mitochondrial genes
• # The [[ operator can add columns to object metadata.
• This is a great place to put QC stats
• pbmc[["[Link]"]] <- PercentageFeatureSet(pbmc, pattern =
"^MT-")
Where are QC metrics stored in Seurat?
• In the example, we visualize QC metrics, and use these to filter cells.
• We filter cells that have unique feature counts over 2,500 or less than
200
• We filter cells that have >5% mitochondrial counts
Visualize QC metrics as a violin plot
• VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "[Link]"),
ncol = 3)
QC plots
• # FeatureScatter is typically used to
visualize feature-feature
relationships, but can be used
• # for anything calculated by the
object, i.e. columns in object
metadata, PC scores etc.
• plot1 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA", feature2 =
"[Link]")
• plot2 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA", feature2 =
"nFeature_RNA")
• plot1 + plot2
Important - subset
• pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA
< 2500 & [Link] < 5)
Normalizing the data
• After removing unwanted cells from the dataset, the next step is to normalize
the data.
• By default, seurat employ a global-scaling normalization method
“LogNormalize” that normalizes the feature expression measurements for
each cell by the total expression, multiplies this by a scale factor (10,000 by
default), and log-transforms the result.
• Normalized values are stored in pbmc[["RNA"]]$data.
pbmc <- NormalizeData(pbmc, [Link] = "LogNormalize",
[Link] = 10000)
pbmc <- NormalizeData(pbmc) (same with default value)
• While this method of normalization is standard and widely used in
scRNA-seq analysis, global-scaling relies on an assumption that each
cell originally contains the same number of RNA molecules.
• There are workflows for the single cell preprocessing that do not
make these assumptions.
Identification of highly variable features (feature
selection – part of a gene selection method)
• Calculate a subset of features that exhibit high cell-to-cell variation in the dataset
(i.e, they are highly expressed in some cells, and lowly expressed in others).
• Researchers found that focusing on these genes in downstream analysis helps to
highlight biological signal in single-cell datasets.
• The method models the mean-variance relationship inherent in single-cell data,
and is implemented in the FindVariableFeatures() function.
• By default, Seurat return 2,000 features per dataset.
• These will be used in downstream analysis, like PCA.
pbmc <- FindVariableFeatures(pbmc, [Link] = "vst", nfeatures = 2000)
variance-stabilizing transformation (VST)
?FindVariableFeatures
• # Identify the 10 most highly
variable genes
• top10 <-
head(VariableFeatures(pbmc), 10)
• # plot variable features with and
without labels
• plot1 <-
VariableFeaturePlot(pbmc)
• plot2 <- LabelPoints(plot = plot1,
points = top10, repel = TRUE)
• plot1 + plot2
Scaling the data
• Seurat apply a linear transformation (‘scaling’) that is a standard pre-
processing step prior to dimensional reduction techniques like PCA.
• The ScaleData() function:
• This step gives equal weight in downstream analyses, so that highly-
expressed genes do not dominate
• The results of this are stored in pbmc[["RNA"]]$[Link]
• By default, only variable features are scaled.
[Link] <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = [Link])
Perform linear dimensional reduction
• Seurat performs PCA on the scaled data.
• By default, only the previously determined variable features are used as input,
but can be defined using features argument if you wish to choose a different
subset (if you do want to use a custom subset of features, make sure you pass
these to ScaleData first).
• For the first principal components, Seurat outputs a list of genes with the most
positive and negative loadings, representing modules of genes that exhibit either
correlation (or anti-correlation) across single-cells in the dataset.
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
PCA Plots
• Seurat provides several useful ways of visualizing both cells and
features that define the PCA, including DimPlot(), and DimHeatmap()
• DimPlot(pbmc, reduction = "pca") + NoLegend()
• DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
Determine the ‘dimensionality’ of the dataset
• To overcome the extensive technical noise in
any single feature for scRNA-seq data, Seurat
clusters cells based on their PCA scores
• The top principal components therefore
represent a robust compression of the dataset.
However, how many components should we
choose to include? 10? 20? 100?
• An heuristic method generates an ‘Elbow plot’:
a ranking of principle components based on the
percentage of variance explained by each one
(ElbowPlot() function).
• In this example, we can observe an ‘elbow’
around PC9-10, suggesting that the majority of
true signal is captured in the first 10 PCs.
ElbowPlot(pbmc)
Knowledge of biological sample and relation
to PC
• Identifying the true dimensionality of a dataset – can be challenging/uncertain for the user.
• Seurat therefore suggests these multiple approaches for users.
• The first is more supervised, exploring PCs to determine relevant sources of heterogeneity
• The second (ElbowPlot)
• The third is a heuristic that is commonly used, and can be calculated instantly. In this example, we
might have been justified in choosing anything between PC 7-12 as a cutoff.
• Here 10 are selected, but encourage users to consider the following:
• Dendritic cell and NK aficionados may recognize that genes strongly associated with PCs 12 and
13 define rare immune subsets (i.e. MZB1 is a marker for plasmacytoid DCs). However, these
groups are so rare, they are difficult to distinguish from background noise for a dataset of this
size without prior knowledge.
• Thus users are encourage to repeat downstream analyses with a different number of PCs (10, 15,
or even 50!).
• As you will observe, the results often do not differ dramatically.
• We advise users to err on the higher side when choosing this parameter. For example, performing
downstream analyses with only 5 PCs does significantly and adversely affect results.
Cluster the cells
• Seurat applies a graph-based clustering approach
• Importantly, the distance metric which drives the
clustering analysis (based on previously identified PCs)
remains the same.
• These methods embed cells in a graph structure - for
example a K-nearest neighbor (KNN) graph, with edges
drawn between cells with similar feature expression
patterns, and then attempt to partition this graph into
highly interconnected ‘quasi-cliques’ or ‘communities’.
Clusters
• Seurat first constructs a KNN graph based on the euclidean distance in PCA space, and refine the
edge weights between any two cells based on the shared overlap in their local neighborhoods
(Jaccard similarity).
• This step is performed using the FindNeighbors() function, and takes as input the previously
defined dimensionality of the dataset (first 10 PCs).
• To cluster the cells, we next apply modularity optimization techniques such as the Louvain
algorithm, to iteratively group cells together, with the goal of optimizing the standard modularity
function.
• The FindClusters() function implements this procedure, and contains a resolution parameter that
sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater
number of clusters.
• We find that setting this parameter between 0.4-1.2 typically returns good results for single-cell
datasets of around 3K cells.
• Optimal resolution often increases for larger datasets. The clusters can be found using the
Idents() function.
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
Run non-linear dimensional reduction
(UMAP/tSNE)
• Seurat offers several non-linear dimensional reduction techniques, such as tSNE
and UMAP, to visualize and explore these datasets.
• The goal of these algorithms is to learn underlying structure in the dataset, in
order to place similar cells together in low-dimensional space.
• Therefore, cells that are grouped together within graph-based clusters
determined before should co-localize on these dimension reduction plots.
• 2D visualization techniques like tSNE and UMAP to be valuable tools for exploring
datasets, but all visualization techniques have limitations, and cannot fully
represent the complexity of the underlying data.
• Seurat encourages users to leverage techniques like UMAP for visualization, but
to avoid drawing biological conclusions solely on the basis of visualization
techniques or directly from clustering.
• Manual intervention with expert biological background knowledge is
indispensable
UMAP
• pbmc <- RunUMAP(pbmc, dims = 1:10)
• DimPlot(pbmc, reduction = "umap")
End of part 1
• #Remaining
• #Gene Selection using Omics Data (by correlation)
• #### DEG across clusters
• #### Annotation
• ### DEG across cell types
• #### Pseudotime
• ### Trajectory analysis