Danger

Please access this document in its canonical location as the currently accessed page may not be rendered correctly: Core plotting functions

Core plotting functions#

Author: Fidel Ramírez

This tutorial explores the visualization possibilities of scanpy and is divided into three sections:

Scatter plots for embeddings (eg. UMAP, t-SNE)
Identification of clusters using known marker genes
Visualization of differentially expressed genes

In this tutorial, we will use a dataset from 10x containing 68k cells from PBMC. Scanpy, includes in its distribution a reduced sample of this dataset consisting of only 700 cells and 765 highly variable genes. This dataset has been already preprocessed and UMAP computed.

In this tutorial, we will also use the following literature markers:

B-cell: CD79A, MS4A1
Plasma: IGJ (JCHAIN)
T-cell: CD3D
NK: GNLY, NKG7
Myeloid: CST3, LYZ
Monocytes: FCGR3A
Dendritic: FCER1A

Scatter plots for embeddings#

With scanpy, scatter plots for tSNE, UMAP and several other embeddings are readily available using the sc.pl.tsne, sc.pl.umap etc. functions. See here the list of options.

Those functions access the data stored in adata.obsm. For example sc.pl.umap uses the information stored in adata.obsm['X_umap']. For more flexibility, any key stored in adata.obsm can be used with the generic function sc.pl.embedding.

import scanpy as sc
from matplotlib.pyplot import rc_context

sc.set_figure_params(dpi=100, color_map="viridis_r")
sc.settings.verbosity = 0
sc.logging.print_header()

scanpy==1.10.0rc2.dev6+g14555ba4.d20240226 anndata==0.11.0.dev78+g64ab900 umap==0.5.5 numpy==1.26.3 scipy==1.11.4 pandas==2.2.0 scikit-learn==1.3.2 statsmodels==0.14.1 igraph==0.10.8 pynndescent==0.5.11

Load pbmc dataset#

pbmc = sc.datasets.pbmc68k_reduced()

# inspect pbmc contents
pbmc

AnnData object with n_obs × n_vars = 700 × 765
    obs: 'bulk_labels', 'n_genes', 'percent_mito', 'n_counts', 'S_score', 'G2M_score', 'phase', 'louvain'
    var: 'n_counts', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
    uns: 'bulk_labels_colors', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

Visualization of gene expression and other variables#

For the scatter plots, the value to plot is given as the color argument. This can be any gene or any column in .obs, where .obs is a DataFrame containing the annotations per observation/cell, see AnnData for more information.

# rc_context is used for the figure size, in this case 4x4
with rc_context({"figure.figsize": (4, 4)}):
    sc.pl.umap(pbmc, color="CD79A")

../_images/57be388755512bf1c4bece588c87eeb6fc5c9961ca87791671a11abb9472dcf6.png

Multiple values can be given to color. In the following example we will plot 6 genes: ‘CD79A’, ‘MS4A1’, ‘IGJ’, CD3D’, ‘FCER1A’, and ‘FCGR3A’ to get an idea on where those marker genes are being expressed.

Also, we will plot two other values: n_counts which is the number of UMI counts per cell (stored in .obs), and bulk_labels which is a categorical value containing the original labelling of the cells from 10X.

The number of plots per row is controlled using the ncols parameter. The maximum value plotted can be adjusted using vmax (similarly vmin can be used for the minimum value). In this case we use p99, which means to use as max value the 99 percentile. The max value can be a number or a list of numbers if the vmax wants to be set for multiple plots individually.

Also, we are using frameon=False to remove the boxes around the plots and s=50 to set the dot size.

color_vars = [
    "CD79A",
    "MS4A1",
    "IGJ",
    "CD3D",
    "FCER1A",
    "FCGR3A",
    "n_counts",
    "bulk_labels",
]
with rc_context({"figure.figsize": (3, 3)}):
    sc.pl.umap(pbmc, color=color_vars, s=50, frameon=False, ncols=4, vmax="p99")

../_images/f5fb8a678ef83c003c35c7b020573f9e2aab56cddf2a5a701e28f2735f83422b.png

In this plot we can see the groups of cells that express the marker genes and the agreement with the original cell labels.

The functions for scatterplots have many options that allow fine tuning of the images. For example, we can look at the clustering as follows:

# compute clusters using the leiden method and store the results with the name `clusters`
sc.tl.leiden(
    pbmc,
    key_added="clusters",
    resolution=0.5,
    n_iterations=2,
    flavor="igraph",
    directed=False,
)

with rc_context({"figure.figsize": (5, 5)}):
    sc.pl.umap(
        pbmc,
        color="clusters",
        add_outline=True,
        legend_loc="on data",
        legend_fontsize=12,
        legend_fontoutline=2,
        frameon=False,
        title="clustering of cells",
        palette="Set1",
    )

../_images/d21985bcdf27832fccc189399bdaef7bee22ef093bdd40da4061dc10b9b1baed.png

Identification of clusters based on known marker genes#

Frequently, clusters need to be labelled using well known marker genes. Using scatter plots we can see the expression of a gene and perhaps associate it with a cluster. Here, we will show other visual ways to associate marker genes to clusters using dotplots, violin plots, heatmaps and something that we call ‘tracksplot’. All of these visualizations summarize the same information, expression split by cluster, and the selection of the best results is left to the investigator do decide.

First, we set up a dictionary with the marker genes, as this will allow scanpy to automatically label the groups of genes:

marker_genes_dict = {
    "B-cell": ["CD79A", "MS4A1"],
    "Dendritic": ["FCER1A", "CST3"],
    "Monocytes": ["FCGR3A"],
    "NK": ["GNLY", "NKG7"],
    "Other": ["IGLL1"],
    "Plasma": ["IGJ"],
    "T-cell": ["CD3D"],
}

dotplot#

A quick way to check the expression of these genes per cluster is to using a dotplot. This type of plot summarizes two types of information: the color represents the mean expression within each of the categories (in this case in each cluster) and the dot size indicates the fraction of cells in the categories expressing a gene.

Also, it is also useful to add a dendrogram to the graph to bring together similar clusters. The hierarchical clustering is computed automatically using the correlation of the PCA components between the clusters.

sc.pl.dotplot(pbmc, marker_genes_dict, "clusters", dendrogram=True)

../_images/2aa081589d93ba1ae17c10722df228186adf87c3c43905fc2b887185d1dde828.png

Using this plot, we can see that cluster 4 correspond to B-cells, cluster 2 is T-cells etc. This information can be used to manually annotate the cells as follows:

# create a dictionary to map cluster to annotation label
cluster2annotation = {
    "0": "Monocytes",
    "1": "NK",
    "2": "T-cell",
    "3": "Dendritic",
    "4": "Dendritic",
    "5": "Plasma",
    "6": "B-cell",
    "7": "Dendritic",
    "8": "Other",
}

# add a new `.obs` column called `cell type` by mapping clusters to annotation using pandas `map` function
pbmc.obs["cell type"] = pbmc.obs["clusters"].map(cluster2annotation).astype("category")

sc.pl.dotplot(pbmc, marker_genes_dict, "cell type", dendrogram=True)

../_images/95d9bc0fcfc65573d6ce802c62532cbda2ce8ecbc4ad06e1f1a731afa50527b3.png

sc.pl.umap(
    pbmc,
    color="cell type",
    legend_loc="on data",
    frameon=False,
    legend_fontsize=10,
    legend_fontoutline=2,
)

../_images/02d019ffed26debe94c69e5adbb17bb52797c8b6b9a7a18512cbf8c7d7fb59ac.png

violin plot#

A different way to explore the markers is with violin plots. Here we can see the expression of CD79A in clusters 5 and 8, and MS4A1 in cluster 5.Compared to a dotplot, the violin plot gives us and idea of the distribution of gene expression values across cells.

with rc_context({"figure.figsize": (4.5, 3)}):
    sc.pl.violin(pbmc, ["CD79A", "MS4A1"], groupby="clusters")

../_images/4d533bc71a9c43e23187261bcaf2bf4886d883f16a387ccd810c160d2568eeb4.png

Note Violin plots can also be used to plot any numerical value stored in .obs. For example, here violin plots are used to compare the number of genes and the percentage of mitochondrial genes between the different clusters.

with rc_context({"figure.figsize": (4.5, 3)}):
    sc.pl.violin(
        pbmc,
        ["n_genes", "percent_mito"],
        groupby="clusters",
        stripplot=False,  # remove the internal dots
        inner="box",  # adds a boxplot inside violins
    )

../_images/cbc036ce11fe42ce6583bcbd1e43942af679d7892b61bca04050a791bacbf83b.png

stacked-violin plot#

To simultaneously look at the violin plots for all marker genes we use sc.pl.stacked_violin. As previously, a dendrogram was added to group similar clusters

ax = sc.pl.stacked_violin(
    pbmc, marker_genes_dict, groupby="clusters", swap_axes=False, dendrogram=True
)

../_images/477a99eb112ecbc4b0e1c526f14ec5a27a1fb25bb6eae3c49e5ec57142142d0d.png

matrixplot#

A simple way to visualize the expression of genes is with a matrix plot. This is a heatmap of the mean expression values per gene grouped by categories. This type plot basically shows the same information as the color in the dotplots.

Here, scale the expression of the genes from 0 to 1, being the maximum mean expression and 0 the minimum.

sc.pl.matrixplot(
    pbmc,
    marker_genes_dict,
    "clusters",
    dendrogram=True,
    cmap="Blues",
    standard_scale="var",
    colorbar_title="column scaled\nexpression",
)

../_images/8bde7de95a2351bd40bce07c97f561d37524ee1a236956626ee6b2b73d5637bd.png

Other useful option is to normalize the gene expression using sc.pp.scale. Here we store this information under the scale layer. Afterwards we adjust the plot min and max and use a diverging color map (in this case RdBu_r where _r means reversed).

# scale and store results in layer
pbmc.layers["scaled"] = sc.pp.scale(pbmc, copy=True).X

sc.pl.matrixplot(
    pbmc,
    marker_genes_dict,
    "clusters",
    dendrogram=True,
    colorbar_title="mean z-score",
    layer="scaled",
    vmin=-2,
    vmax=2,
    cmap="RdBu_r",
)

../_images/6b10b32141a9383cc2d3ec112957ac1cfbee8154c449683ff9be2f612420b6f0.png

Combining plots in subplots#

An axis can be passed to a plot to combine multiple outputs as in the following example

import matplotlib.pyplot as plt

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 4), gridspec_kw={"wspace": 0.9})

ax1_dict = sc.pl.dotplot(
    pbmc, marker_genes_dict, groupby="bulk_labels", ax=ax1, show=False
)
ax2_dict = sc.pl.stacked_violin(
    pbmc, marker_genes_dict, groupby="bulk_labels", ax=ax2, show=False
)
ax3_dict = sc.pl.matrixplot(
    pbmc, marker_genes_dict, groupby="bulk_labels", ax=ax3, show=False, cmap="viridis"
)

../_images/de3758d79d14fbdb767d6a7bd7e4168274b76a9681e24a4e9fd8ecbb3a422b66.png

Heatmaps#

Heatmaps do not collapse cells as in previous plots. Instead, each cells is shown in a row (or column if swap_axes=True). The groupby information can be added and is shown using the same color code found for sc.pl.umap or any other embedding.

ax = sc.pl.heatmap(
    pbmc, marker_genes_dict, groupby="clusters", cmap="viridis", dendrogram=True
)

../_images/b092a195a2e9e7d9ca1d9584a76b5b83ff3ad24c61bf77ffc0b61536bb8a4265.png

The heatmap can also be plotted on scaled data. In the next image, similar to the previus matrixplot the min and max had been adjusted and a divergent color map is used.

ax = sc.pl.heatmap(
    pbmc,
    marker_genes_dict,
    groupby="clusters",
    layer="scaled",
    vmin=-2,
    vmax=2,
    cmap="RdBu_r",
    dendrogram=True,
    swap_axes=True,
    figsize=(11, 4),
)

../_images/72018721bf8762bbb23bdd5d118ce8c53e823e86f38104920ff825660e8251a2.png

Tracksplot#

The track plot shows the same information as the heatmap, but, instead of a color scale, the gene expression is represented by height.

ax = sc.pl.tracksplot(pbmc, marker_genes_dict, groupby="clusters", dendrogram=True)

../_images/94d6e876d08841a146bf45d65240823eb8a13b2445b54ed983aea66c19f19d74.png

Visualization of marker genes#

Instead of characterizing clusters by known gene markers as previously, we can identify genes that are differentially expressed in the clusters or groups.

To identify differentially expressed genes we run sc.tl.rank_genes_groups. This function will take each group of cells and compare the distribution of each gene in a group against the distribution in all other cells not in the group. Here, we will use the original cell labels given by 10x to identify marker genes for those cell types.

sc.tl.rank_genes_groups(pbmc, groupby="clusters", method="wilcoxon")

Visualize marker genes using dotplot#

The dotplot visualization is useful to get an overview of the genes that show differential expression. To make the resulting image more compact we will use n_genes=4 to show only the top 4 scoring genes.

sc.pl.rank_genes_groups_dotplot(pbmc, n_genes=4)

../_images/eca426ba0c175565e5e34a5fb0e16eb75ea1933a4008b87634d6251931bab373.png

In order to get a better representation we can plot log fold changes instead of gene expression. Also, we want to focus on genes that have a log fold change >= 3 between the cell type expression and the rest of cells.

In this case we set values_to_plot='logfoldchanges' and min_logfoldchange=3.

Because log fold change is a divergent scale we also adjust the min and max to be plotted and use a divergent color map. Notice in the following plot that is rather difficult to distinguish between T-cell populations.

sc.pl.rank_genes_groups_dotplot(
    pbmc,
    n_genes=4,
    values_to_plot="logfoldchanges",
    min_logfoldchange=3,
    vmax=7,
    vmin=-7,
    cmap="bwr",
)

../_images/8a5b9a3a844761620a81333f1f9aa834f42f58dbe1c07556eede789fa73d3821.png

Focusing on particular groups#

Next, we use a dotplot focusing only on two groups (the groups option is also available for violin, heatmap and matrix plots). Here, we set n_genes=30 as in this case it will show all the genes that have a min_logfoldchange=4 up to 30.

sc.pl.rank_genes_groups_dotplot(
    pbmc,
    n_genes=30,
    values_to_plot="logfoldchanges",
    min_logfoldchange=4,
    vmax=7,
    vmin=-7,
    cmap="bwr",
    groups=["3", "7"],
)

../_images/5cba17b1487a5b9f8e4f82747a7b5402b238afd0f8b9771329cc2cafe762323f.png

Visualize marker genes using matrixplot#

For the following plot the we use the previously computed ‘scaled’ values (stored in layer scaled) and use a divergent color map.

sc.pl.rank_genes_groups_matrixplot(
    pbmc, n_genes=3, use_raw=False, vmin=-3, vmax=3, cmap="bwr", layer="scaled"
)

../_images/d92ac8cff005d3d675db226476fcf16ce95c0103d3028eb2770e7aa8c5ae5f78.png

Visualize marker genes using stacked violin plots#

sc.pl.rank_genes_groups_stacked_violin(pbmc, n_genes=3, cmap="viridis_r")

../_images/33587a56c76354a9006e5432138fc644f0098910297dce64d387667ac6ad53f6.png

Visualize marker genes using heatmap#

sc.pl.rank_genes_groups_heatmap(
    pbmc,
    n_genes=3,
    use_raw=False,
    swap_axes=True,
    vmin=-3,
    vmax=3,
    cmap="bwr",
    layer="scaled",
    figsize=(10, 7),
    show=False,
);

../_images/a8f32a9fad1c3bca6c889411fa6340403d93cff435aa8365331db6d409594068.png

Showing 10 genes per category, turning the gene labels off and swapping the axes. Notice that when the image is swapped, a color code for the categories appear instead of the ‘brackets’.

sc.pl.rank_genes_groups_heatmap(
    pbmc,
    n_genes=10,
    use_raw=False,
    swap_axes=True,
    show_gene_labels=False,
    vmin=-3,
    vmax=3,
    cmap="bwr",
)

../_images/26e16bbaa5daced14abe7dc8619f966a7bc30f364dba20e461e1a8b3f272ba98.png

Visualize marker genes using tracksplot#

sc.pl.rank_genes_groups_tracksplot(pbmc, n_genes=3)

../_images/e1d9a1ff1adb3487a37c60f6580f69b4d4c0da2a2d01025db06e83cca852729b.png

Comparison of marker genes using split violin plots#

In scanpy, is very easy to compare marker genes using split violin plots for all groups at once.

with rc_context({"figure.figsize": (9, 1.5)}):
    sc.pl.rank_genes_groups_violin(pbmc, n_genes=20, jitter=False)

../_images/100284609ffeda2f93a6b1a417142135c4b7efe38bdf3c2883769701a7b4299d.png

../_images/9e13c66be5f38fec2706fb5171b52fed4f99fce233114d52b1dec726fbb68961.png

../_images/4d34fb0bcd64d9de622e9376144625d1cee775569df40d6c2953ef9c5c7eb989.png

../_images/23de59d1185f3141aad01dd923dd8751ace35d6223b6dd3283d195b4b5a35bda.png

../_images/d392ef9bf7bccb0cfa92bd7cfd20f15708623f461850d487d268ace8881bacbe.png

../_images/9be6cbb75c066ca431aa1697c02b3bbc6387d91d32376e5cb2fff533a28d0dcf.png

../_images/0f430485100fa1e738391c6c573147391e3b05235873271ca76308ff27873f8b.png

../_images/07218b82400958fd118f851eba7a1c593d3d083f790cbb9b46c10ea4dffa173c.png

../_images/3d1d0aaab00e772893316173eaca085db98615b901a1ac7ade4d25dff68a672e.png

Dendrogram options#

Most of the visualizations can arrange the categories using a dendrogram. However, the dendrogram can also be plotted independently as follows:

# compute hierarchical clustering using PCs (several distance metrics and linkage methods are available).
sc.tl.dendrogram(pbmc, "bulk_labels")

ax = sc.pl.dendrogram(pbmc, "bulk_labels")

../_images/d468b69af1c7f0af35798c70c6312a99e836c25fd2134bb9edbdbbe4b7e9cded.png

Plot correlation#

Together with the dendrogram it is possible to plot the correlation (by default ‘pearson’) of the categories.

ax = sc.pl.correlation_matrix(pbmc, "bulk_labels", figsize=(5, 3.5))

../_images/b12877b4e058b5310108e83e7b492a837cf49fb4203e4b68e48c678b612e34d8.png

Core plotting functions

Contents

Core plotting functions#

Scatter plots for embeddings#

Load pbmc dataset#

Visualization of gene expression and other variables#

Identification of clusters based on known marker genes#

dotplot#

violin plot#

stacked-violin plot#

matrixplot#

Combining plots in subplots#

Heatmaps#

Tracksplot#

Visualization of marker genes#

Visualize marker genes using dotplot#

Focusing on particular groups#

Visualize marker genes using matrixplot#

Visualize marker genes using stacked violin plots#

Visualize marker genes using heatmap#

Visualize marker genes using tracksplot#

Comparison of marker genes using split violin plots#

Dendrogram options#

Plot correlation#