Quality Control Analysis

R Programming

Python Programming

Single-cell RNA-seq

Quality Control

This lesson guides participants through performing quality control analysis for single‑cell RNA‑seq data. Participants learn to generate key QC metrics, visualize and interpret them and apply informed filtering thresholds to remove low‑quality cells and genes. By the end of the lesson, participants will be able to assess data complexity, mitochondrial contamination and overall cell quality to prepare a clean dataset for downstream analysis.

Authors

Mary Piper

Meeta Mistry

Radhika Khetani

Lorena Pantano

Jihe Liu

Will Gammerdinger

Noor Sohail

Published

June 17, 2026

Keywords

R, Python, QC metrics, filtering, mitochondrial ratio, complexity

Approximate time: 90 minutes

Learning objectives

In this lesson, we will:

Construct quality control metrics and visually evaluate the quality of the data.
Apply appropriate filters to remove low quality cells.

Overview of lesson

At this point, we have loaded our dataset into our programming language. We saw in the metadata, that there are several columns that will help us evaluate the quality of our cells. In single-cell experiments, it is common to have some cells that are dead/dying or are generally low quality. Therefore, we want to remove these noisy cells before we run downstream analyses. This next step uses our metadata columns to filter cells using quality metrics and clean up the dataset.

Figure 1: Overview of the single-cell RNA-seq workflow.

Updating metadata

When data is loaded and the initial object is created, we have some basic metadata to work from for each of the cells in the count matrix. To take a close look at this metadata, let’s view the data frame stored in the metadata slot of our merged object:

R
Python

# Explore merged metadata
View(merged_seurat@meta.data)

Table 1: Inspecting our metadata table to get a sense for the type of information stored within it.

	orig.ident	nCount_RNA	nFeature_RNA
ctrl_AAACATACAATGCC-1	ctrl	2344	874
ctrl_AAACATACATTTCC-1	ctrl	3125	896
ctrl_AAACATACCAGAAA-1	ctrl	2578	725
ctrl_AAACATACCAGCTA-1	ctrl	3261	979
ctrl_AAACATACCATGCA-1	ctrl	746	362
ctrl_AAACATACCTCGCT-1	ctrl	3519	866

We need to calculate the basic QC metrics that we will use to filter cells later. We will accomplish this with the pp.calculate_qc_metrics() function to automatically generate scores for each cell.

# Calculate QC metrics
sc.pp.calculate_qc_metrics(merged_adata,
                           percent_top = None,
                           log1p = False,
                           inplace = True)
# Explore merged metadata
merged_adata.obs.head()

Table 2: Inspecting our metadata table to get a sense for the type of information stored within it.

	n_genes	sample	n_genes_by_counts	total_counts
AAACATACAATGCC-1_ctrl	874	ctrl	874	2344
AAACATACATTTCC-1_ctrl	896	ctrl	896	3125
AAACATACCAGAAA-1_ctrl	725	ctrl	725	2578
AAACATACCAGCTA-1_ctrl	979	ctrl	979	3261
AAACATACCATGCA-1_ctrl	362	ctrl	362	746

The columns of interest right now are:

Identity (orig.ident or sample)
Number of UMIs per cell (nCount_RNA or total_counts)
Number of genes detected per cell (nFeature_RNA or n_genes)

However, we would also like to include some additional information that would be useful to have in our metadata including cell IDs and condition information.

Storing metadata as separate dataframe

When working with a large dataset, it can be helpful to extract the metadata dataframe into a separate variable. In this way we can work with this as a separate entity from the object without the risk of affecting any other data stored inside.

R
Python

When we added columns of information to our metadata file above, we simply added it directly to the metadata slot in the Seurat object using the $ operator.

Additionally we are going to create another column titled sample which is also stored in orig.ident. We do this because this because sample is a commonly used metadata field that others will look for when they download your data.

# Add cell IDs to metadata
merged_seurat$cells <- rownames(merged_seurat@meta.data)

# Create sample column
merged_seurat$sample <- merged_seurat$orig.ident

Now you are all set up with the metrics you need to assess the quality of your data! Your final metadata table will have rows that correspond to each cell, and columns with information about those cells:

Table 3: Inspecting the first rows of our metadata table to ensure that they look formatted correctly.

	orig.ident	nCount_RNA	nFeature_RNA	cells	sample
ctrl_AAACATACAATGCC-1	ctrl	2344	874	ctrl_AAACATACAATGCC-1	ctrl
ctrl_AAACATACATTTCC-1	ctrl	3125	896	ctrl_AAACATACATTTCC-1	ctrl
ctrl_AAACATACCAGAAA-1	ctrl	2578	725	ctrl_AAACATACCAGAAA-1	ctrl
ctrl_AAACATACCAGCTA-1	ctrl	3261	979	ctrl_AAACATACCAGCTA-1	ctrl
ctrl_AAACATACCATGCA-1	ctrl	746	362	ctrl_AAACATACCATGCA-1	ctrl

# Add cell IDs to metadata
merged_adata.obs["cells"] = merged_adata.obs.index

Table 4: Inspecting the first rows of our metadata table to ensure that they look formatted correctly.

	n_genes	sample	n_genes_by_counts	total_counts	cells
AAACATACAATGCC-1_ctrl	874	ctrl	874	2344	AAACATACAATGCC-1_ctrl
AAACATACATTTCC-1_ctrl	896	ctrl	896	3125	AAACATACATTTCC-1_ctrl
AAACATACCAGAAA-1_ctrl	725	ctrl	725	2578	AAACATACCAGAAA-1_ctrl
AAACATACCAGCTA-1_ctrl	979	ctrl	979	3261	AAACATACCAGCTA-1_ctrl
AAACATACCATGCA-1_ctrl	362	ctrl	362	746	AAACATACCATGCA-1_ctrl

Assessing the quality metrics

Now that we have generated the various metrics to assess, we can explore them with visualizations. We will assess various metrics and then decide on which cells are low quality and should be removed from the analysis. Some of these values have already been calculated:

Cell counts
UMI counts per cell
Genes detected per cell

In order to create the appropriate plots for the quality control analysis, we need to calculate some additional metrics. These include:

Number of genes detected per UMI: this metric will give us an idea of the complexity of our dataset (more genes detected per UMI, more complex our data)
Mitochondrial ratio: this metric will give us a percentage of cell reads originating from the mitochondrial genes

What about doublets?

In single-cell RNA sequencing experiments, doublets are generated from two cells. They typically arise due to errors in cell sorting or capture, especially in droplet-based protocols involving thousands of cells. Doublets are obviously undesirable when the aim is to characterize populations at the single-cell level. In particular, they can incorrectly suggest the existence of intermediate populations or transitory states that do not actually exist. Thus, it is desirable to remove doublets so that they do not compromise interpretation of the results.

Many workflows use maximum thresholds for UMIs or genes, with the idea that a much higher number of reads or genes detected indicate multiple cells. While this rationale seems to be intuitive, it is not accurate. Also, many of the tools used to detect doublets tend to get rid of cells with intermediate or continuous phenotypes, although they may work well on datasets with very discrete cell types. Scrublet is a popular tool for doublet detection.

Currently, we recommend not including any thresholds at this point in time. When we have identified markers for each of the clusters, we suggest exploring the markers to determine whether the markers apply to more than one cell type.

Cell counts

The cell counts are determined by the number of unique cellular barcodes detected. For this experiment, between 12,000 -13,000 cells are expected.

In an ideal world, you would expect the number of unique cellular barcodes to correspond to the number of cells you loaded. However, this is not the case as capture rates of cells are only a proportion of what is loaded. For example, the inDrops cell capture efficiency is higher (70-80%) compared to 10X which is between 50-60%.

Automated cell counting

The capture efficiency could appear much lower if the cell concentration used for library preparation was not accurate. Cell concentration should NOT be determined by FACS machine or Bioanalyzer (these tools are not accurate for concentration determination), instead use a hemocytometer or automated cell counter for calculation of cell concentration.

The cell numbers can also vary by protocol, producing cell numbers that are much higher than what we loaded. For example, during the inDrops protocol, the cellular barcodes are present in the hydrogels, which are encapsulated in the droplets with a single cell and lysis/reaction mixture. While each hydrogel should have a single cellular barcode associated with it, occasionally a hydrogel can have more than one cellular barcode. Similarly, with the 10X protocol there is a chance of obtaining only a barcoded bead in the emulsion droplet (GEM) and no actual cell. Both of these, in addition to the presence of dying cells can lead to a higher number of cellular barcodes than cells.

R
Python

# Visualize the number of cell counts per sample
merged_seurat@meta.data %>% 
    ggplot(aes(x = sample, fill = sample)) + 
    geom_bar() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
    ggtitle("NCells")

Figure 2: Barplot of the number of cells per sample.

fig, ax = plt.subplots()
ax = sns.countplot(data = merged_adata.obs, 
                   x = "sample", hue = "sample",
                   ax = ax)
ax.set_title("Number of Cells per Sample")

# Add counts above the bars
for container in ax.containers:
    ax.bar_label(container)

plt.show()

Figure 3: Barplot of the number of cells per sample.

We see over 15,000 cells per sample, which is quite a bit more than the 12-13,000 expected. It is clear that we likely have some low-quality or empty-droplet “junk” cells present.

UMI counts (transcripts) per cell

The UMI counts per cell should generally be above 500, that is the low end of what we expect. If UMI counts are between 500-1,000 counts, it is usable but the cells probably should have been sequenced more deeply.

R
Python

# Visualize the number UMIs/transcripts per cell
merged_seurat@meta.data %>% 
    ggplot(aes(color = sample, 
                x = nCount_RNA, 
                fill = sample)) + 
    geom_density(alpha = 0.2) + 
    scale_x_log10() + 
    theme_classic() +
    ylab("Cell density") +
    geom_vline(xintercept = 500)

Figure 4: Density plot showing the distribution of UMI counts per cell for each sample

# Visualize the number UMIs/transcripts per cell
fig, ax = plt.subplots()
ax = sns.histplot(merged_adata.obs, 
                  x = "total_counts", 
                  hue = "sample",
                  log_scale = True, # apply log scale to counts
                  alpha = 0.2, ax = ax)
ax.set_title("Total Counts")

# Add line showing potential threshold
ax.axvline(500, color="black")
plt.show()

Figure 5: Histogram plot showing the distribution of UMI counts per cell for each sample

We can see that majority of our cells in both samples have 1,000 UMIs or greater, which is great.

Genes detected per cell

We have similar expectations for gene detection as for UMI detection, although it may be a bit lower than UMIs. For high quality data, the proportional histogram should contain a single large peak that represents cells that were encapsulated. If we see a small shoulder to the left of the major peak (not present in our data), or a bimodal distribution of the cells, that can indicate a couple of things. It might be that there are a set of cells that failed for some reason. It could also be that there are biologically different types of cells (i.e. quiescent cell populations, less complex cells of interest), and/or one type is much smaller than the other (i.e. cells with high counts may be cells that are larger in size). Therefore, this threshold should be assessed with other metrics that we describe in this lesson.

R
Python

# Visualize the distribution of genes detected per cell via histogram
merged_seurat@meta.data %>%
    ggplot(aes(color = sample,
               x = nFeature_RNA,
               fill = sample)) +
    geom_density(alpha = 0.2) +
    theme_classic() +
    scale_x_log10() +
    geom_vline(xintercept = 300)

Figure 6: Density plot showing the distribution of genes detected per cell for each sample

# Visualize the distribution of genes detected per cell via histogram
ax = sns.histplot(merged_adata.obs, 
                  x = "n_genes", 
                  hue = "sample",
                  log_scale = True, # apply log scale to counts
                  alpha = 0.2)
ax.set_title("Genes per Cell")

# Add line showing potential threshold
ax.axvline(300, color = "black")
plt.show()

Figure 7: Histogram plot showing the distribution of genes detected per cell for each sample

Complexity (novelty) score

We can evaluate each cell in terms of how complex the RNA species are by using a measure called the novelty score. The novelty score is computed by taking the ratio of nGenes over nUMI. If there are many captured transcripts (high nUMI) and a low number of genes detected in a cell, this likely means that you only captured a low number of genes and simply sequenced transcripts from those lower number of genes over and over again. These low complexity (low novelty) cells could represent a specific cell type (i.e. red blood cells which lack a typical transcriptome), or could be due to an artifact or contamination. Generally, we expect the novelty score to be above 0.80 for good quality cells.

The novelty score is computed as a ratio of genes to UMIs, as shown below:

\[ \text{Complexity Score} = \frac{\log_{10}(\text{Number of Genes})}{\log_{10}(\text{Number of UMIs})} \]

This value is quite easy to calculate, as we take the log10 of the number of genes detected per cell and the log10 of the number of UMIs per cell, then divide the log10 number of genes by the log10 number of UMIs. The novelty score and how it relates to the complexity of the RNA species is described in more detail later in this lesson.

R
Python

# Add number of genes per UMI for each cell to metadata
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)

After calculating the complexity score, we can now visualize the distribution of the score for each sample in our dataset:

# Visualize the overall complexity of the gene expression by visualizing 
# the genes detected per UMI (novelty score)
merged_seurat@meta.data %>%
    ggplot(aes(x = log10GenesPerUMI,
               color = sample, fill = sample)) +
    geom_density(alpha = 0.2) +
    theme_classic() +
    geom_vline(xintercept = 0.8)

Figure 8: Density plot showing the overall complexity of gene expression per cell for each sample

merged_adata.obs["log10GenesPerUMI"] = np.log10(merged_adata.obs["n_genes"]) / np.log10(merged_adata.obs["total_counts"])

After calculating the complexity score, we can now visualize the distribution of the score for each sample in our dataset:

ax = sns.histplot(merged_adata.obs, 
                  x = "log10GenesPerUMI", hue = "sample",
                  alpha = 0.2)
ax.set_title("Complexity")

# Add line showing potential threshold
ax.axvline(0.8, color = "black")
plt.show()

Figure 9: Histogram plot showing the overall complexity of gene expression per cell for each sample

Mitochondrial Ratio

This metric can identify whether there is a large amount of mitochondrial contamination from dead or dying cells. We define poor quality samples for mitochondrial counts as cells which surpass the 0.2 mitochondrial ratio mark, unless of course you are expecting this in your sample.

Think about your biological question

While using a baseline score of 0.20 is an acceptable threshold for removing high mitochondrial content cells, it is important to always go back to your original biological question. What samples are you working with? Do you expect there to be high values of mitochondrial expression due to your experimental condition?

For example, if you were studying renal oncocytomas, would you make this same choice? This disease is characterized as having aberrantly high mitochondrial expression, so would it make sense to remove cells with high mitochondrial ratio?

For our analysis, rather than using a percentage value we would prefer to work with the ratio value. As such, we will reverse that last step performed by the function by taking the output value and dividing by 100. This ratio is computed as:

\[ \text{Mitochondrial Ratio} = \frac{\text{Number of reads aligning to mitochondrial genes}} {\text{Total reads}} \]

R
Python

Seurat has a convenient function that allows us to calculate the proportion of transcripts mapping to mitochondrial genes. The PercentageFeatureSet() function takes in a pattern argument and searches through all gene identifiers in the dataset for that pattern. Since we are looking for mitochondrial genes, we are searching any gene identifiers that begin with the pattern “MT-”. For each cell, the function takes the sum of counts across all genes (features) belonging to the “Mt-” set, and then divides by the count sum for all genes (features). This value is multiplied by 100 to obtain a percentage value.

# Compute percent mito ratio
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, 
                                                pattern = "^MT-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100

We can once again visualize this metric for each sample:

# Visualize the distribution of mitochondrial gene expression detected per cell
merged_seurat@meta.data %>% 
    ggplot(aes(color = sample, x = mitoRatio, fill = sample)) + 
    geom_density(alpha = 0.2) + 
    scale_x_log10() + 
    theme_classic() +
    geom_vline(xintercept = 0.2)

Figure 10: Density plot showing the distribution of mitochondrial gene expression detected per cell for each sample.

The first step is to identify which genes are mitochondrial. We do this my identify any genes (from the .var) that starts with the characters MT- and assign True/False values to genes that fit that criteria.

# Identify mitochondrial genes
mito_genes = merged_adata.var_names.str.startswith('MT-')
merged_adata.var['mito'] = mito_genes

merged_adata.var.head()

Table 5: Inspecting our gene metadata table to see our mitochondrial information.

	gene_ids	feature_types	n_cells_by_counts	mean_counts	pct_dropout_by_counts	total_counts	mito
MIR1302-2HG	ENSG00000243485	Gene Expression	0	0	100	0	False
FAM138A	ENSG00000237613	Gene Expression	0	0	100	0	False
OR4F5	ENSG00000186092	Gene Expression	0	0	100	0	False
AL627309.1	ENSG00000238009	Gene Expression	12	0.000381631	99.9618	12	False
AL627309.3	ENSG00000239945	Gene Expression	1	3.18026e-05	99.9968	1	False

We then use the qc_vars parameter in pp.calculate_qc_metrics() to tell the function to calculate metrics (including number of counts) that belong to genes that have True set in the mito column.

# Calculate number of counts coming from mito_genes
sc.pp.calculate_qc_metrics(merged_adata, 
                           qc_vars = ["mito"],
                           inplace = True)

# We want the value as a ratio
merged_adata.obs["mito_ratio"] = merged_adata.obs["pct_counts_mito"] / 100

We can once again visualize this metric for each sample:

# Histogram of mitochondrial ratio
ax = sns.histplot(merged_adata.obs, 
                  x = "mito_ratio", hue = "sample",
                  alpha = 0.2)
ax.set_title("Mitochondrial Ratio")

# Add line showing potential threshold
ax.axvline(0.2, color = "black")
plt.show()

Figure 11: Histogram plot showing the distribution of mitochondrial gene expression detected per cell for each sample.

MT genes for organisms other than human

The pattern provided (“^MT-”) works for human gene names. You may need to adjust the pattern argument depending on your organism of interest. Additionally, if you weren’t using gene names as the gene ID then this function wouldn’t work as we have used it above as the pattern will not suffice. Since there are caveats to using this function, it is advisable to manually compute this metric. If you are interested, we have code available to compute this metric on your own.

Joint filtering effects

Considering any of these QC metrics in isolation can lead to misinterpretation of cellular signals. For example, cells with a comparatively high fraction of mitochondrial counts may be involved in respiratory processes and may be cells that you would like to keep. Likewise, other metrics can have other biological interpretations. A general rule of thumb when performing QC is to set thresholds for individual metrics to be as permissive as possible, and always consider the joint effects of these metrics. In this way, you reduce the risk of filtering out any viable cell populations.

Two metrics that are often evaluated together are the number of UMIs and the number of genes detected per cell. Here, we have plotted the number of genes versus the number of UMIs coloured by the fraction of mitochondrial reads. Jointly visualizing the count and gene thresholds and additionally overlaying the mitochondrial fraction, gives a summarized perspective of the quality per cell.

R
Python

# Visualize the correlation between genes detected and number of UMIs and 
# determine whether strong presence of cells with low numbers of genes/UMIs
merged_seurat@meta.data %>% 
    ggplot(aes(x = nCount_RNA, 
               y = nFeature_RNA, 
               color = mitoRatio)) + 
    geom_point() + 
      scale_colour_gradient(low = "gray90", high = "black") +
    stat_smooth(method = lm) +
    scale_x_log10() + 
    scale_y_log10() + 
    theme_classic() +
    geom_vline(xintercept = 500) +
    geom_hline(yintercept = 250) +
    facet_wrap(~sample)

Figure 12: Scatterplot contrasting nUMI and nGenes while coloring each cell by mitochondrial ratio.

# Order dataset so values are plotted 
# such that largest values are on top
meta = merged_adata.obs.sort_values('mito_ratio', ascending=True)

g = sns.FacetGrid(meta, col="sample", height=5)
g.map_dataframe(sns.scatterplot, 
                x = "total_counts", 
                y = "n_genes", 
                hue = "mito_ratio",
                palette = "viridis")

# Scale x,y axes by log10
for ax in g.axes.flat:
    ax.set(xscale = "log", yscale = "log")

# Add threshold lines
g.refline(x = 500, y = 300)

# Add colorbar for the legend 
sm = g.axes[0, 0].collections[0]
plt.colorbar(sm, ax = g.axes, 
             label = "mito_ratio")
plt.show()

Figure 13: Scatterplot contrasting nUMI and nGenes while coloring each cell by mitochondrial ratio.

Good cells will generally exhibit both higher number of genes per cell and higher numbers of UMIs (upper right quadrant of the plot). Cells that are poor quality are likely to have low genes and UMIs per cell, and correspond to the data points in the bottom left quadrant of the plot. With this plot we also evaluate the slope of the line, and any scatter of data points in the bottom right hand quadrant of the plot. These cells have a high number of UMIs but only a few number of genes. These could be dying cells, but also could represent a population of a low-complexity celltype (i.e. red blood cells).

Mitochondrial read fractions are only high in particularly low count cells with few detected genes (darker colored data points). This could be indicative of damaged/dying cells whose cytoplasmic mRNA has leaked out through a broken membrane, and thus, only mRNA located in the mitochondria is still conserved. We can see from the plot that these cells are filtered out by our count and gene number thresholds.

Filtering

Cell-level filtering

Now that we have visualized the various metrics, we can decide on the thresholds to apply which will result in the removal of low quality cells. Often the recommendations mentioned earlier are a rough guideline, and the specific experiment needs to inform the exact thresholds chosen. We will use the following thresholds:

nUMI > 500
nGene > 250
log10GenesPerUMI > 0.8
mitoRatio < 0.2

R
Python

# Filter out low quality cells using selected 
# thresholds - these will change with experiment
filtered_seurat <- subset(x = merged_seurat, 
                          subset = (nCount_RNA >= 500) &
                                   (nFeature_RNA >= 250) &
                                   (log10GenesPerUMI > 0.80) &
                                   (mitoRatio < 0.20))
filtered_seurat

An object of class Seurat 
33538 features across 29629 samples within 1 assay 
Active assay: RNA (33538 features, 0 variable features)
 1 layer present: counts

Now let’s filter our cells based on our metadata columns.

filtered_adata = merged_adata[merged_adata.obs["total_counts"] >= 500]
filtered_adata = filtered_adata[filtered_adata.obs["n_genes"] >= 250]
filtered_adata = filtered_adata[filtered_adata.obs["log10GenesPerUMI"] >= 0.8]
filtered_adata = filtered_adata[filtered_adata.obs["mito_ratio"] < 0.2]

filtered_adata

View of AnnData object with n_obs × n_vars = 29629 × 33538
    obs: 'n_genes', 'sample', 'n_genes_by_counts', 'total_counts', 'cells', 'log10GenesPerUMI', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mito', 'log1p_total_counts_mito', 'pct_counts_mito', 'mito_ratio'
    var: 'gene_ids', 'feature_types', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'mito', 'log1p_mean_counts', 'log1p_total_counts'
    layers: None (.X)

Gene-level filtering

Within our data we will have many genes with zero counts. These genes can dramatically reduce the average expression for a cell and so we will remove them from our data. We will start by identifying which genes have a zero count in each cell:

R
Python

# Extract counts
counts <- GetAssayData(object = filtered_seurat, layer = "counts")

# Output a logical matrix specifying for each gene on whether or not there are more than zero counts per cell
nonzero <- counts > 0

Now, we will perform some filtering by prevalence. If a gene is only expressed in a handful of cells, it is not particularly meaningful as it still brings down the averages for all other cells it is not expressed in. For our data we choose to only keep genes which are expressed in 10 or more cells. By using this filter, genes which have zero counts in all cells will effectively be removed.

# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10

# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts[keep_genes, ]

Finally, take those filtered counts and create a new Seurat object for downstream analysis.

# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, 
                                      meta.data = filtered_seurat@meta.data)
filtered_seurat

An object of class Seurat 
14065 features across 29629 samples within 1 assay 
Active assay: RNA (14065 features, 0 variable features)
 1 layer present: counts

# Remove lowly expressed genes
sc.pp.filter_genes(filtered_adata, 
                   inplace = True, 
                   min_counts = 10)

filtered_adata

AnnData object with n_obs × n_vars = 29629 × 14167
    obs: 'n_genes', 'sample', 'n_genes_by_counts', 'total_counts', 'cells', 'log10GenesPerUMI', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mito', 'log1p_total_counts_mito', 'pct_counts_mito', 'mito_ratio'
    var: 'gene_ids', 'feature_types', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'mito', 'log1p_mean_counts', 'log1p_total_counts', 'n_counts'
    layers: None (.X)

Exercise 1

After performing the filtering, it’s recommended to look back over the metrics to make sure that your data matches your expectations and is good for downstream analysis.

Perform all of the same QC plots using the filtered data.
Report the number of cells left for each sample, and comment on whether the number of cells removed is high or low. Can you give reasons why this number is still not ~12K (which is how many cells were loaded for the experiment)?
After filtering for nGene per cell, you should still observe a small shoulder to the right of the main peak. What might this shoulder represent?
When plotting the nGene against nUMI do you observe any data points in the bottom right quadrant of the plot? What can you say about these cells that have been removed?

Saving filtered object

Based on these QC metrics we would identify any failed samples and move forward with our filtered cells. Often we iterate through the QC metrics using different filtering criteria; it is not necessarily a linear process. When satisfied with the filtering criteria, we would save our filtered cell object for clustering and marker identification.

R
Python

# Create RDS object to load at any time
saveRDS(filtered_seurat, file="data/seurat_filtered.rds")

# Create h5ad object to load at any time
filtered_adata.write_h5ad("data/filtered_adata.h5ad")

Bad data

The data we are working with is pretty good quality. If you are interested in knowing what ‘bad’ data might look like when performing QC, we have some materials linked here where we explore similar QC metrics of a poor quality sample.

Next Lesson >>

Back to Schedule

Reuse

CC-BY-4.0

--- title: "Quality Control Analysis" description: | This lesson guides participants through performing quality control analysis for single‑cell RNA‑seq data. Participants learn to generate key QC metrics, visualize and interpret them and apply informed filtering thresholds to remove low‑quality cells and genes. By the end of the lesson, participants will be able to assess data complexity, mitochondrial contamination and overall cell quality to prepare a clean dataset for downstream analysis. author: - Mary Piper - Meeta Mistry - Radhika Khetani - Lorena Pantano - Jihe Liu - Will Gammerdinger - Noor Sohail date: "2026-06-17" categories: - R Programming - Python Programming - Single-cell RNA-seq - Quality Control keywords: - R - Python - QC metrics - filtering - mitochondrial ratio - complexity license: "CC-BY-4.0" editor_options: markdown: wrap: 72 --- ```{r} #| label: load_libraries_data_R #| echo: false # Load libraries and data library(Seurat) library(tidyverse) library(knitr) # Load merged object merged_seurat <- readRDS("intermediate/03_merged_seurat.RDS") ``` ```{python} #| label: load_libraries_data_py #| echo: false # Load libraries and data import scanpy as sc import numpy as np import seaborn as sns import matplotlib.pyplot as plt # Load merged object merged_adata = sc.read_h5ad("intermediate/03_merged_adata.h5ad") ``` Approximate time: 90 minutes ## Learning objectives In this lesson, we will: - Construct quality control metrics and visually evaluate the quality of the data. - Apply appropriate filters to remove low quality cells. ## Overview of lesson At this point, we have loaded our dataset into our programming language. We saw in the metadata, that there are several columns that will help us evaluate the quality of our cells. In single-cell experiments, it is common to have some cells that are dead/dying or are generally low quality. Therefore, we want to remove these noisy cells before we run downstream analyses. This next step uses our metadata columns to _filter cells using quality metrics_ and clean up the dataset. ::: {#fig-sc-workflow .figure} ![](../img/sc_workflow_2022.jpg){width=630} Overview of the single-cell RNA-seq workflow. ::: ## Updating metadata When data is loaded and the initial object is created, we have some basic metadata to work from for each of the cells in the count matrix. To take a close look at this metadata, let's view the data frame stored in the `metadata` slot of our `merged` object: ::: {.panel-tabset group="language"} ## R ```{r} #| label: inspect_metadata #| eval: false # Explore merged metadata View(merged_seurat@meta.data) ``` ```{r} #| label: tbl-table_metadata #| tbl-cap: "Inspecting our metadata table to get a sense for the type of information stored within it." #| echo: false merged_seurat@meta.data %>% head() %>% kable() ``` ## Python We need to calculate the basic QC metrics that we will use to filter cells later. We will accomplish this with the `pp.calculate_qc_metrics()` function to automatically generate scores for each cell. ```{python} #| label: calc_qc_metrics_py #| eval: false # Calculate QC metrics sc.pp.calculate_qc_metrics(merged_adata, percent_top = None, log1p = False, inplace = True) # Explore merged metadata merged_adata.obs.head() ``` ```{python} #| label: tbl-show_metadata_head_py #| tbl-cap: "Inspecting our metadata table to get a sense for the type of information stored within it." #| echo: false #| results: asis sc.pp.calculate_qc_metrics(merged_adata, percent_top = None, log1p = False, inplace = True) print(merged_adata.obs.head().to_markdown()) ``` ::: The columns of interest right now are: - Identity (`orig.ident` or `sample`) - Number of UMIs per cell (`nCount_RNA` or `total_counts`) - Number of genes detected per cell (`nFeature_RNA` or `n_genes`) However, we would also like to **include some additional information** that would be useful to have in our metadata including **cell IDs and condition information**. ::: callout-note # Storing metadata as separate dataframe When working with a large dataset, it can be helpful to extract the `metadata` dataframe into a separate variable. In this way we can work with this as a separate entity from the object without the risk of affecting any other data stored inside. ::: ::: {.panel-tabset group="language"} ## R When we added columns of information to our metadata file above, we simply added it directly to the metadata slot in the Seurat object using the `$` operator. Additionally we are going to create another column titled `sample` which is also stored in `orig.ident`. We do this because this because `sample` is a commonly used metadata field that others will look for when they download your data. ```{r} #| label: update_meta # Add cell IDs to metadata merged_seurat$cells <- rownames(merged_seurat@meta.data) # Create sample column merged_seurat$sample <- merged_seurat$orig.ident ``` Now you are **all set up with the metrics you need to assess the quality of your data**! Your final metadata table will have rows that correspond to each cell, and columns with information about those cells: ```{r} #| label: tbl-table_renamed_metadata #| tbl-cap: "Inspecting the first rows of our metadata table to ensure that they look formatted correctly." #| echo: false merged_seurat@meta.data %>% head(n = 5) %>% knitr::kable() ``` ## Python ```{python} #| label: update_meta_py # Add cell IDs to metadata merged_adata.obs["cells"] = merged_adata.obs.index ``` Now you are **all set up with the metrics you need to assess the quality of your data**! Your final metadata table will have rows that correspond to each cell, and columns with information about those cells: ```{python} #| label: tbl-table_renamed_metadata_py #| tbl-cap: "Inspecting the first rows of our metadata table to ensure that they look formatted correctly." #| echo: false #| results: asis print(merged_adata.obs.head().to_markdown()) ``` ::: ## Assessing the quality metrics Now that we have generated the various metrics to assess, we can explore them with visualizations. We will assess various metrics and then decide on which cells are low quality and should be removed from the analysis. Some of these values have already been calculated: - Cell counts - UMI counts per cell - Genes detected per cell In order to create the appropriate plots for the quality control analysis, we need to calculate some additional metrics. These include: - **Number of genes detected per UMI:** this metric will give us an idea of the complexity of our dataset (more genes detected per UMI, more complex our data) - **Mitochondrial ratio:** this metric will give us a percentage of cell reads originating from the mitochondrial genes ::: {.callout-note collapse="true"} # What about doublets? In single-cell RNA sequencing experiments, doublets are generated from two cells. They typically arise due to errors in cell sorting or capture, especially in droplet-based protocols involving thousands of cells. Doublets are obviously undesirable when the aim is to characterize populations at the single-cell level. In particular, they can incorrectly suggest the existence of intermediate populations or transitory states that do not actually exist. Thus, it is desirable to remove doublets so that they do not compromise interpretation of the results. Many workflows use maximum thresholds for UMIs or genes, with the idea that a much higher number of reads or genes detected indicate multiple cells. While this rationale seems to be intuitive, it is not accurate. Also, many of the tools used to detect doublets tend to get rid of cells with intermediate or continuous phenotypes, although they may work well on datasets with very discrete cell types. [Scrublet](https://github.com/AllonKleinLab/scrublet) is a popular tool for doublet detection. Currently, we recommend not including any thresholds at this point in time. When we have identified markers for each of the clusters, we suggest exploring the markers to determine whether the markers apply to more than one cell type. ::: ### Cell counts The cell counts are determined by the number of unique cellular barcodes detected. For this experiment, between 12,000 -13,000 cells are expected. In an ideal world, you would expect the number of unique cellular barcodes to correspond to the number of cells you loaded. However, this is not the case as capture rates of cells are only a proportion of what is loaded. For example, the inDrops cell **capture efficiency** is higher (70-80%) compared to 10X which is between 50-60%. ::: {.callout-note collapse=true} # Automated cell counting The capture efficiency could appear much lower if the cell concentration used for library preparation was not accurate. Cell concentration should NOT be determined by FACS machine or Bioanalyzer (these tools are not accurate for concentration determination), instead use a hemocytometer or automated cell counter for calculation of cell concentration. ::: The cell numbers can also vary by protocol, **producing cell numbers that are much higher than what we loaded**. For example, during the inDrops protocol, the cellular barcodes are present in the hydrogels, which are encapsulated in the droplets with a single cell and lysis/reaction mixture. While each hydrogel should have a single cellular barcode associated with it, occasionally a hydrogel can have more than one cellular barcode. Similarly, with the 10X protocol there is a chance of obtaining only a barcoded bead in the emulsion droplet (GEM) and no actual cell. Both of these, in addition to the presence of dying cells can lead to a higher number of cellular barcodes than cells. ::: {.panel-tabset group="language"} ## R ```{r} #| label: fig-cell_counts_per_sample #| fig-cap: Barplot of the number of cells per sample. # Visualize the number of cell counts per sample merged_seurat@meta.data %>% ggplot(aes(x = sample, fill = sample)) + geom_bar() + theme_classic() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + theme(plot.title = element_text(hjust = 0.5, face = "bold")) + ggtitle("NCells") ``` ## Python ```{python} #| label: fig-plot_ncells_py #| fig-cap: Barplot of the number of cells per sample. fig, ax = plt.subplots() ax = sns.countplot(data = merged_adata.obs, x = "sample", hue = "sample", ax = ax) ax.set_title("Number of Cells per Sample") # Add counts above the bars for container in ax.containers: ax.bar_label(container) plt.show() ``` ::: We see over 15,000 cells per sample, which is quite a bit more than the 12-13,000 expected. It is clear that we likely have some low-quality or empty-droplet “junk” cells present. ### UMI counts (transcripts) per cell The UMI counts per cell should generally be above 500, that is the low end of what we expect. If UMI counts are between 500-1,000 counts, it is usable but the cells probably should have been sequenced more deeply. ::: {.panel-tabset group="language"} ## R ```{r} #| label: fig-UMIs_per_cell #| fig-cap: Density plot showing the distribution of UMI counts per cell for each sample # Visualize the number UMIs/transcripts per cell merged_seurat@meta.data %>% ggplot(aes(color = sample, x = nCount_RNA, fill = sample)) + geom_density(alpha = 0.2) + scale_x_log10() + theme_classic() + ylab("Cell density") + geom_vline(xintercept = 500) ``` ## Python ```{python} #| label: fig-UMIs_per_cell_py #| fig-cap: Histogram plot showing the distribution of UMI counts per cell for each sample # Visualize the number UMIs/transcripts per cell fig, ax = plt.subplots() ax = sns.histplot(merged_adata.obs, x = "total_counts", hue = "sample", log_scale = True, # apply log scale to counts alpha = 0.2, ax = ax) ax.set_title("Total Counts") # Add line showing potential threshold ax.axvline(500, color="black") plt.show() ``` ::: We can see that majority of our cells in both samples have 1,000 UMIs or greater, which is great. ### Genes detected per cell We have similar expectations for gene detection as for UMI detection, although it may be a bit lower than UMIs. For high quality data, the proportional histogram should contain **a single large peak that represents cells that were encapsulated**. If we see a **small shoulder** to the left of the major peak (not present in our data), or a bimodal distribution of the cells, that can indicate a couple of things. It might be that there are a set of **cells that failed** for some reason. It could also be that there are **biologically different types of cells** (i.e. quiescent cell populations, less complex cells of interest), and/or one type is much smaller than the other (i.e. cells with high counts may be cells that are larger in size). Therefore, this threshold should be assessed with other metrics that we describe in this lesson. ::: {.panel-tabset group="language"} ## R ```{r} #| label: fig-genes_detected #| fig-cap: Density plot showing the distribution of genes detected per cell for each sample # Visualize the distribution of genes detected per cell via histogram merged_seurat@meta.data %>% ggplot(aes(color = sample, x = nFeature_RNA, fill = sample)) + geom_density(alpha = 0.2) + theme_classic() + scale_x_log10() + geom_vline(xintercept = 300) ``` ## Python ```{python} #| label: fig-genes_detected_py #| fig-cap: Histogram plot showing the distribution of genes detected per cell for each sample # Visualize the distribution of genes detected per cell via histogram ax = sns.histplot(merged_adata.obs, x = "n_genes", hue = "sample", log_scale = True, # apply log scale to counts alpha = 0.2) ax.set_title("Genes per Cell") # Add line showing potential threshold ax.axvline(300, color = "black") plt.show() ``` ::: ### Complexity (novelty) score We can evaluate each cell in terms of how complex the RNA species are by using a measure called the novelty score. The novelty score is computed by taking the **ratio of nGenes over nUMI**. If there are many captured transcripts (high nUMI) and a low number of genes detected in a cell, this likely means that you only captured a low number of genes and simply sequenced transcripts from those lower number of genes over and over again. These low complexity (low novelty) cells could represent a specific cell type (i.e. red blood cells which lack a typical transcriptome), or could be due to an artifact or contamination. Generally, we expect the novelty score to be above 0.80 for good quality cells. The novelty score is computed as a ratio of genes to UMIs, as shown below: $$ \text{Complexity Score} = \frac{\log_{10}(\text{Number of Genes})}{\log_{10}(\text{Number of UMIs})} $$ This value is quite easy to calculate, as we take the log10 of the number of genes detected per cell and the log10 of the number of UMIs per cell, then divide the log10 number of genes by the log10 number of UMIs. The novelty score and how it relates to the complexity of the RNA species is described in more detail later in this lesson. ::: {.panel-tabset group="language"} ## R ```{r} #| label: create_novelty_score # Add number of genes per UMI for each cell to metadata merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA) ``` After calculating the complexity score, we can now visualize the distribution of the score for each sample in our dataset: ```{r} #| label: fig-complexity_score_plot #| fig-cap: Density plot showing the overall complexity of gene expression per cell for each sample # Visualize the overall complexity of the gene expression by visualizing # the genes detected per UMI (novelty score) merged_seurat@meta.data %>% ggplot(aes(x = log10GenesPerUMI, color = sample, fill = sample)) + geom_density(alpha = 0.2) + theme_classic() + geom_vline(xintercept = 0.8) ``` ## Python ```{python} #| label: create_novelty_score_py merged_adata.obs["log10GenesPerUMI"] = np.log10(merged_adata.obs["n_genes"]) / np.log10(merged_adata.obs["total_counts"]) ``` After calculating the complexity score, we can now visualize the distribution of the score for each sample in our dataset: ```{python} #| label: fig-complexity_score_plot_py #| fig-cap: Histogram plot showing the overall complexity of gene expression per cell for each sample ax = sns.histplot(merged_adata.obs, x = "log10GenesPerUMI", hue = "sample", alpha = 0.2) ax.set_title("Complexity") # Add line showing potential threshold ax.axvline(0.8, color = "black") plt.show() ``` ::: ### Mitochondrial Ratio This metric can identify whether there is a large amount of **mitochondrial contamination from dead or dying cells**. We define poor quality samples for mitochondrial counts as cells which surpass the 0.2 mitochondrial ratio mark, unless of course you are expecting this in your sample. ::: {.callout-note collapse=true} # Think about your biological question While using a baseline score of 0.20 is an acceptable threshold for removing high mitochondrial content cells, it is important to always go back to your original biological question. What samples are you working with? Do you expect there to be high values of mitochondrial expression due to your experimental condition? For example, if you were studying renal oncocytomas, would you make this same choice? This disease is characterized as having [aberrantly high mitochondrial expression](https://www.nature.com/articles/modpathol2015101), so would it make sense to remove cells with high mitochondrial ratio? ::: For our analysis, rather than using a percentage value we would prefer to work with the ratio value. As such, we will reverse that last step performed by the function by taking the output value and dividing by 100. This ratio is computed as: $$ \text{Mitochondrial Ratio} = \frac{\text{Number of reads aligning to mitochondrial genes}} {\text{Total reads}} $$ ::: {.panel-tabset group="language"} ## R Seurat has a convenient function that allows us to calculate the **proportion of transcripts mapping to mitochondrial genes**. The `PercentageFeatureSet()` function takes in a `pattern` argument and searches through all gene identifiers in the dataset for that pattern. Since we are looking for mitochondrial genes, we are searching any gene identifiers that begin with the pattern "MT-". For each cell, the function takes the sum of counts across all genes (features) belonging to the "Mt-" set, and then divides by the count sum for all genes (features). This value is multiplied by 100 to obtain a percentage value. ```{r} #| label: compute_mito_ratio # Compute percent mito ratio merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-") merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100 ``` We can once again visualize this metric for each sample: ```{r} #| label: fig-mito_ratio_plot #| fig-cap: Density plot showing the distribution of mitochondrial gene expression detected per cell for each sample. # Visualize the distribution of mitochondrial gene expression detected per cell merged_seurat@meta.data %>% ggplot(aes(color = sample, x = mitoRatio, fill = sample)) + geom_density(alpha = 0.2) + scale_x_log10() + theme_classic() + geom_vline(xintercept = 0.2) ``` ## Python The first step is to identify which genes are mitochondrial. We do this my identify any genes (from the `.var`) that starts with the characters `MT-` and assign `True`/`False` values to genes that fit that criteria. ```{python} #| label: label_mito #| eval: false # Identify mitochondrial genes mito_genes = merged_adata.var_names.str.startswith('MT-') merged_adata.var['mito'] = mito_genes merged_adata.var.head() ``` ```{python} #| label: tbl-show_var_head_py #| tbl-cap: "Inspecting our gene metadata table to see our mitochondrial information." #| echo: false #| results: asis mito_genes = merged_adata.var_names.str.startswith('MT-') merged_adata.var['mito'] = mito_genes print(merged_adata.var.head().to_markdown()) ``` We then use the `qc_vars` parameter in `pp.calculate_qc_metrics()` to tell the function to calculate metrics (including number of counts) that belong to genes that have `True` set in the `mito` column. ```{python} #| label: calc_mito_ratio # Calculate number of counts coming from mito_genes sc.pp.calculate_qc_metrics(merged_adata, qc_vars = ["mito"], inplace = True) # We want the value as a ratio merged_adata.obs["mito_ratio"] = merged_adata.obs["pct_counts_mito"] / 100 ``` We can once again visualize this metric for each sample: ```{python} #| label: fig-mito_ratio_plot_py #| fig-cap: Histogram plot showing the distribution of mitochondrial gene expression detected per cell for each sample. # Histogram of mitochondrial ratio ax = sns.histplot(merged_adata.obs, x = "mito_ratio", hue = "sample", alpha = 0.2) ax.set_title("Mitochondrial Ratio") # Add line showing potential threshold ax.axvline(0.2, color = "black") plt.show() ``` ::: ::: {.callout-note collapse="true"} # MT genes for organisms other than human The pattern provided ("^MT-") works for human gene names. You may need to adjust the pattern argument depending on your organism of interest. Additionally, if you weren't using gene names as the gene ID then this function wouldn't work as we have used it above as the pattern will not suffice. Since there are caveats to using this function, it is advisable to manually compute this metric. If you are interested, we have [code available to compute this metric on your own](https://github.com/hbctraining/scRNA-seq/blob/master/lessons/mitoRatio.md). ::: ### Joint filtering effects Considering any of these QC metrics in isolation can lead to misinterpretation of cellular signals. For example, cells with a comparatively high fraction of mitochondrial counts may be involved in respiratory processes and may be cells that you would like to keep. Likewise, other metrics can have other biological interpretations. A general rule of thumb when performing QC is to **set thresholds for individual metrics to be as permissive as possible, and always consider the joint effects** of these metrics. In this way, you reduce the risk of filtering out any viable cell populations. Two metrics that are often evaluated together are the number of UMIs and the number of genes detected per cell. Here, we have plotted the **number of genes versus the number of UMIs coloured by the fraction of mitochondrial reads**. Jointly visualizing the count and gene thresholds and additionally overlaying the mitochondrial fraction, gives a summarized perspective of the quality per cell. ::: {.panel-tabset group="language"} ## R ```{r} #| label: fig-joint_effects_plot #| fig-cap: Scatterplot contrasting nUMI and nGenes while coloring each cell by mitochondrial ratio. # Visualize the correlation between genes detected and number of UMIs and # determine whether strong presence of cells with low numbers of genes/UMIs merged_seurat@meta.data %>% ggplot(aes(x = nCount_RNA, y = nFeature_RNA, color = mitoRatio)) + geom_point() + scale_colour_gradient(low = "gray90", high = "black") + stat_smooth(method = lm) + scale_x_log10() + scale_y_log10() + theme_classic() + geom_vline(xintercept = 500) + geom_hline(yintercept = 250) + facet_wrap(~sample) ``` ## Python ```{python} #| label: fig-joint_effects_plot_code_py #| fig-cap: Scatterplot contrasting nUMI and nGenes while coloring each cell by mitochondrial ratio. #| echo: true #| eval: false # Order dataset so values are plotted # such that largest values are on top meta = merged_adata.obs.sort_values('mito_ratio', ascending=True) g = sns.FacetGrid(meta, col="sample", height=5) g.map_dataframe(sns.scatterplot, x = "total_counts", y = "n_genes", hue = "mito_ratio", palette = "viridis") # Scale x,y axes by log10 for ax in g.axes.flat: ax.set(xscale = "log", yscale = "log") # Add threshold lines g.refline(x = 500, y = 300) # Add colorbar for the legend sm = g.axes[0, 0].collections[0] plt.colorbar(sm, ax = g.axes, label = "mito_ratio") plt.show() ``` ```{python} #| label: fig-joint_effects_plot_py #| fig-cap: Scatterplot contrasting nUMI and nGenes while coloring each cell by mitochondrial ratio. #| echo: false # Need to turn this into a function # Because intermediate plots are shown otherwise def joint_effects_plot(adata): # Order dataset so values are plotted such that largest values are on top meta = adata.obs.sort_values('mito_ratio', ascending=True) g = sns.FacetGrid(meta, col="sample", height=5) g.map_dataframe( sns.scatterplot, x="total_counts", y="n_genes", hue="mito_ratio", palette="viridis" ) # Scale x,y axes by log10 for ax in g.axes.flat: ax.set(xscale="log", yscale="log") # Add threshold lines g.refline(x=500, y=300) # Add colorbar for the legend sm = g.axes[0, 0].collections[0] plt.colorbar(sm, ax=g.axes, label="mito_ratio") return g g = joint_effects_plot(merged_adata) plt.show() ``` ::: Good cells will generally exhibit both higher number of genes per cell and higher numbers of UMIs (upper right quadrant of the plot). Cells that are **poor quality are likely to have low genes and UMIs per cell**, and correspond to the data points in the bottom left quadrant of the plot. With this plot we also evaluate the **slope of the line**, and any scatter of data points in the **bottom right hand quadrant** of the plot. These cells have a high number of UMIs but only a few number of genes. These could be dying cells, but also could represent a population of a low-complexity celltype (i.e. red blood cells). **Mitochondrial read fractions are only high in particularly low count cells with few detected genes** (darker colored data points). This could be indicative of damaged/dying cells whose cytoplasmic mRNA has leaked out through a broken membrane, and thus, only mRNA located in the mitochondria is still conserved. We can see from the plot that these cells are filtered out by our count and gene number thresholds. ## Filtering ### Cell-level filtering Now that we have visualized the various metrics, we can decide on the thresholds to apply which will result in the removal of low quality cells. Often the recommendations mentioned earlier are a rough guideline, and the specific experiment needs to inform the exact thresholds chosen. We will use the following thresholds: - nUMI > 500 - nGene > 250 - log10GenesPerUMI > 0.8 - mitoRatio < 0.2 ::: {.panel-tabset group="language"} To filter, we will go back to our Seurat object and use the `subset()` function: ## R ```{r} #| label: filtered_seurat_object # Filter out low quality cells using selected # thresholds - these will change with experiment filtered_seurat <- subset(x = merged_seurat, subset = (nCount_RNA >= 500) & (nFeature_RNA >= 250) & (log10GenesPerUMI > 0.80) & (mitoRatio < 0.20)) filtered_seurat ``` ## Python Now let's filter our cells based on our metadata columns. ```{python} #| label: filtered_adata_object filtered_adata = merged_adata[merged_adata.obs["total_counts"] >= 500] filtered_adata = filtered_adata[filtered_adata.obs["n_genes"] >= 250] filtered_adata = filtered_adata[filtered_adata.obs["log10GenesPerUMI"] >= 0.8] filtered_adata = filtered_adata[filtered_adata.obs["mito_ratio"] < 0.2] filtered_adata ``` ::: ### Gene-level filtering Within our data we will have many genes with zero counts. These genes can dramatically reduce the average expression for a cell and so we will remove them from our data. We will start by identifying which genes have a zero count in each cell: ::: {.panel-tabset group="language"} ## R ```{r} #| label: extract_non_zero_counts # Extract counts counts <- GetAssayData(object = filtered_seurat, layer = "counts") # Output a logical matrix specifying for each gene on whether or not there are more than zero counts per cell nonzero <- counts > 0 ``` Now, we will perform some filtering by prevalence. If a gene is only expressed in a handful of cells, it is not particularly meaningful as it still brings down the averages for all other cells it is not expressed in. For our data we choose to **only keep genes which are expressed in 10 or more cells.** By using this filter, genes which have zero counts in all cells will effectively be removed. ```{r} #| label: retain_expressed_genes # Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene keep_genes <- Matrix::rowSums(nonzero) >= 10 # Only keeping those genes expressed in more than 10 cells filtered_counts <- counts[keep_genes, ] ``` Finally, take those filtered counts and create a new Seurat object for downstream analysis. ```{r} #| label: reassign_filtered_counts # Reassign to filtered Seurat object filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data) filtered_seurat ``` ## Python ```{python} #| label: filter_genes # Remove lowly expressed genes sc.pp.filter_genes(filtered_adata, inplace = True, min_counts = 10) filtered_adata ``` ::: ::: callout-tip # [Exercise 1](05_quality_control-Answer_key.qmd#exercise-1) After performing the filtering, it's recommended to look back over the metrics to make sure that your data matches your expectations and is good for downstream analysis. 1. Perform all of the same QC plots using the filtered data. 2. Report the number of cells left for each sample, and comment on whether the number of cells removed is high or low. Can you give reasons why this number is still not ~12K (which is how many cells were loaded for the experiment)? 3. After filtering for nGene per cell, you should still observe a small shoulder to the right of the main peak. What might this shoulder represent? 4. When plotting the nGene against nUMI do you observe any data points in the bottom right quadrant of the plot? What can you say about these cells that have been removed? ::: ## Saving filtered object Based on these QC metrics we would identify any failed samples and move forward with our filtered cells. Often we iterate through the QC metrics using different filtering criteria; it is not necessarily a linear process. When satisfied with the filtering criteria, we would save our filtered cell object for clustering and marker identification. ::: {.panel-tabset group="language"} ## R ```{r} #| label: save_filtered_seurat_object #| eval: false # Create RDS object to load at any time saveRDS(filtered_seurat, file="data/seurat_filtered.rds") ``` ## Python ```{python} #| label: save_filtered_adata_object #| eval: false # Create h5ad object to load at any time filtered_adata.write_h5ad("data/filtered_adata.h5ad") ``` ::: ```{python} #| label: py_intructor_save #| eval: false #| echo: false filtered_adata.write_h5ad("intermediate/05_filtered_adata.h5ad") ``` ```{r} #| label: R_instructor_save #| eval: false #| echo: false saveRDS(filtered_seurat, "intermediate/05_filtered_seurat.RDS") ``` ::: callout-note # Bad data The data we are working with is pretty good quality. If you are interested in knowing what 'bad' data might look like when performing QC, we have some materials [linked here](Aside_QC_bad_data.qmd) where we explore similar QC metrics of a poor quality sample. ::: *** [Next Lesson >>](06_theory_of_PCA.qmd) [Back to Schedule](../schedule/schedule.qmd)