Infer CNV on lung cancer dataset

In this tutorial, we demonstrate how infercnvpy can be used to derive Copy Number Variation (CNV) from a single-cell RNA-seq dataset and to distinguish between tumor and normal cells.

For this tutorial, we use the dataset by [MMR+20] generated on the Smart-seq2 platform. The original dataset contains about 20,000 cells. Here, we use a downsampled version with 3,000 cells which is available through infercnvpy.datasets.maynard2020_3k().

Warning

We consider this method still experimental, but decided to already put it on GitHub because it might be useful. Treat the results with care. In particular, no validation with ground-truth data has been performed.

import scanpy as sc
import infercnvpy as cnv
import matplotlib.pyplot as plt

sc.settings.set_figure_params(figsize=(5, 5))

/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/traitlets/traitlets.py:3258: FutureWarning: --rc={'figure.dpi': 96} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
  warn(

sc.logging.print_header()

scanpy==1.9.1 anndata==0.8.0 umap==0.5.3 numpy==1.22.4 scipy==1.8.1 pandas==1.4.2 scikit-learn==1.1.1 statsmodels==0.13.2 python-igraph==0.9.10 pynndescent==0.5.7

Loading the example dataset

Note

Preprocessing data

Low-quality cells should already be filtered out and the input data must be normalized and log-transformed. For more information, see Preparing input data.

Also, the genomic positions need to be stored in adata.var. The columns chromosome, start, and end hold the chromosome and the start and end positions on that chromosome for each gene, respectively.

Infercnvpy provides the infercnvpy.io.genomic_position_from_gtf() function to read these information from a GTF file and add them to adata.var.

The example dataset is already appropriately preprocessed.

adata = cnv.datasets.maynard2020_3k()
adata.var.loc[:, ["ensg", "chromosome", "start", "end"]].head()

100%|██████████| 49.5M/49.5M [00:00<00:00, 101MB/s]

	ensg	chromosome	start	end
symbol
AL645933.5	ENSG00000288587.1	chr6	31400702	31463705
AC010184.1	ENSG00000288585.1	chr3	141449745	141456434
AC023296.1	ENSG00000288580.1	chr8	2923568	2926689
AL117334.2	ENSG00000288577.1	chr20	3406380	3410036
AC107294.4	ENSG00000288576.1	chr3	184778723	184780720

Let’s first inspect the UMAP plot based on the transcriptomics data:

sc.pl.umap(adata, color="cell_type")

Running infercnv

Let’s now run infercnvpy.tl.infercnv(). Essentially, this method sorts genes by chromosome and genomic position and compares the average gene expression over genomic region to a reference. The original inferCNV method uses a window size of 100, but larger window sizes can make sense, depending on the number of genes in your dataset.

infercnv() adds a cell x genomic_region matrix to adata.obsm[“X_cnv”].

For more information about the method check out The inferCNV method.

Note

Choosing reference cells

The most common use-case is to compare tumor against normal cells. If you have prior information about which cells are normal (e.g. from cell-type annotations based on transcriptomics data), it is recommended to provide this information to infercnv().

The more different cell-types you can provide, the better. Some cell-types physiologicallly over-express certain genomic regions (e.g. plasma cells highly express Immunoglobulin genes which are genomically adjacent). If you provide multiple cell-types, only regions are considered being subject to CNV that are different from all provided cell-types.

If you don’t provide any reference, the mean of all cells is used instead, which may work well on datasets that contain enough tumor and normal cells.

# We provide all immune cell types as "normal cells".
cnv.tl.infercnv(
    adata,
    reference_key="cell_type",
    reference_cat=[
        "B cell",
        "Macrophage",
        "Mast cell",
        "Monocyte",
        "NK cell",
        "Plasma cell",
        "T cell CD4",
        "T cell CD8",
        "T cell regulatory",
        "mDC",
        "pDC",
    ],
    window_size=250,
)

100%|██████████| 1/1 [00:31<00:00, 31.22s/it]

Now, we can plot smoothed gene expression by cell-type and chromosome. We can observe that the Epithelial cell cluster, which consists largely of tumor cells, appears to be subject to copy number variation.

cnv.pl.chromosome_heatmap(adata, groupby="cell_type")

/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/infercnvpy/pl/_chromosome_heatmap.py:59: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  tmp_adata = AnnData(X=adata.obsm[f"X_{use_rep}"], obs=adata.obs)

Clustering by CNV profiles and identifying tumor cells

To cluster and annotate cells, infercnvpy mirrors the scanpy workflow. The following functions work exactely as their scanpy counterpart, except that they use the CNV profile matrix as input. Using these functions, we can perform graph-based clustering and generate a UMAP plot based on the CNV profiles. Based on these clusters, we can annotate tumor and normal cells.

infercnvpy.tl.pca(adata[, svd_solver, …])	Compute the PCA on the result of infercnvpy.tl.infercnv().
infercnvpy.pp.neighbors(adata[, use_rep, …])	Compute the neighborhood graph based on the result from infercnvpy.tl.infercnv().
infercnvpy.tl.leiden(adata[, neighbors_key, …])	Perform leiden clustering on the CNV neighborhood graph.
infercnvpy.tl.umap(adata[, neighbors_key, …])	Compute the UMAP on the result of infercnvpy.tl.infercnv().
infercnvpy.pl.umap(adata, **kwargs)	Plot the CNV UMAP.

cnv.tl.pca(adata)
cnv.pp.neighbors(adata)
cnv.tl.leiden(adata)

After running leiden clustering, we can plot the chromosome heatmap by CNV clusters. We can observe that, as opposted to the clusters at the bottom, the clusters at the top have essentially no differentially expressed genomic regions. The differentially expressed regions are likely due to copy number variation and the respective clusters likely represent tumor cells.

cnv.pl.chromosome_heatmap(adata, groupby="cnv_leiden", dendrogram=True)

/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/infercnvpy/pl/_chromosome_heatmap.py:59: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  tmp_adata = AnnData(X=adata.obsm[f"X_{use_rep}"], obs=adata.obs)

WARNING: dendrogram data not found (using key=dendrogram_cnv_leiden). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.
WARNING: You’re trying to run this on 5314 dimensions of `.X`, if you really want this, set `use_rep='X'`.
         Falling back to preprocessing with `sc.pp.pca` and default params.
WARNING: Groups are not reordered because the `groupby` categories and the `var_group_labels` are different.
categories: 0, 1, 2, etc.
var_group_labels: chr1, chr2, chr3, etc.

UMAP plot of CNV profiles

We can visualize the same clusters as a UMAP plot. Additionally, infercnvpy.tl.cnv_score() computes a summary score that quantifies the amount of copy number variation per cluster. It is simply defined as the mean of the absolute values of the CNV matrix for each cluster.

cnv.tl.umap(adata)
cnv.tl.cnv_score(adata)

The UMAP plot consists of a large blob of “normal” cells and several smaller clusters with distinct CNV profiles. Except for cluster “12”, which consists of ciliated cells, the isolated clusters are all epithelial cells. These are likely tumor cells and each cluster represents an individual sub-clone.

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(11, 11))
ax4.axis("off")
cnv.pl.umap(
    adata,
    color="cnv_leiden",
    legend_loc="on data",
    legend_fontoutline=2,
    ax=ax1,
    show=False,
)
cnv.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
cnv.pl.umap(adata, color="cell_type", ax=ax3)

We can also visualize the CNV score and clusters on the transcriptomics-based UMAP plot. Again, we can see that there are subclusters of epithelial cells that belong to a distinct CNV cluster, and that these clusters tend to have the highest CNV score.

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    2, 2, figsize=(12, 11), gridspec_kw=dict(wspace=0.5)
)
ax4.axis("off")
sc.pl.umap(adata, color="cnv_leiden", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
sc.pl.umap(adata, color="cell_type", ax=ax3)

Classifying tumor cells

Based on these observations, we can now assign cell to either “tumor” or “normal”. To this end, we add a new column cnv_status to adata.obs.

adata.obs["cnv_status"] = "normal"
adata.obs.loc[
    adata.obs["cnv_leiden"].isin(["10", "13", "15", "5", "11", "16", "12"]), "cnv_status"
] = "tumor"

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), gridspec_kw=dict(wspace=0.5))
cnv.pl.umap(adata, color="cnv_status", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_status", ax=ax2)

... storing 'cnv_status' as categorical

Now, we can plot the CNV heatmap for tumor and normal cells separately:

cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "tumor", :])

/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/infercnvpy/pl/_chromosome_heatmap.py:59: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  tmp_adata = AnnData(X=adata.obsm[f"X_{use_rep}"], obs=adata.obs)

cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "normal", :])

/home/docs/checkouts/readthedocs.org/user_builds/holonet-tutorial/envs/latest/lib/python3.9/site-packages/infercnvpy/pl/_chromosome_heatmap.py:59: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  tmp_adata = AnnData(X=adata.obsm[f"X_{use_rep}"], obs=adata.obs)