10X单细胞（10X空间转录组）CNV分析之inferCNVpy

阳了个阳~~~?

文章在10X单细胞（10X空间转录组）CNV分析回顾之CopyKAT详细回顾了copycat，还有分享的文章copyKAT推断单细胞转录组肿瘤细胞CNV（自动识别肿瘤normal和tumor），但是copyCAT的引用率还是比较低，相比于inferCNV 2000多的引用率，还是差很多，所以，inferCNV还是主流，我们这一篇就来回顾一下python版本的inferCNVpy。

其实inferCNV已经是非常老的一个软件了，14年就发表了文章，所以我们来分享其最新的python版本。

安装

pip install infercnvpy

加载,与scanpy无缝衔接

import scanpy as scimport infercnvpy as cnvimport matplotlib.pyplot as plt
sc.settings.set_figure_params(figsize=(5, 5))

Loading the example dataset

前处理

• 应该已经过滤掉低质量的细胞，并且必须对输入数据进行归一化和对数转换
• 此外，基因组位置需要存储在 adata.var 中。 "chromosome", "start", "end"分别保存，每个基因的染色体以及该染色体上的开始和结束位置。
• Infercnvpy provides the infercnvpy.io.genomic_position_from_gtf() function to read these information from a GTF file and add them to adata.var.

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

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

Running infercnv

现在运行 infercnvpy.tl.infercnv()。本质上，该方法通过染色体和基因组位置对基因进行分类，并将基因组区域的平均基因表达与参考进行比较。原始的 inferCNV 方法使用上下游50个基因作为窗口，但更大的窗口大小可能有意义，具体取决于数据集中的基因数量。

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

Choosing reference cells(最为关键的一步)

• 最常见的用例是将肿瘤与正常细胞进行比较。如果有关于哪些细胞正常的先验信息（例如，来自基于转录组学数据的细胞类型注释），建议将此信息提供给 infercnv()。
• 可以提供的不同细胞类型越多越好。一些细胞类型在生理上过度表达某些基因组区域（例如浆细胞高度表达基因组相邻的免疫球蛋白基因）。如果提供多种细胞类型，则仅考虑与所有提供的细胞类型不同的区域受CNV影响。
• 如果不提供任何参考，则使用所有细胞的平均值，这可能适用于包含足够肿瘤和正常细胞的数据集。

# 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,)

现在，可以按细胞类型和染色体绘制平滑的基因表达。 可以观察到主要由肿瘤细胞组成的上皮细胞cluster似乎受到拷贝数变化的影响。

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

Clustering by CNV profiles and identifying tumor cells

为了对细胞进行聚类和注释，infercnvpy 结合了 scanpy 工作流程。以下函数与它们的 scanpy 对应函数完全一样，除了使用 CNV 分析矩阵作为输入。使用这些函数，可以执行基于图形的聚类并根据CNV分析矩阵生成 UMAP 图。基于这些聚类，可以注释肿瘤和正常细胞。

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

运行 leiden 聚类后，可以通过 CNV 聚类绘制染色体热图。可以观察到，与底部的clusters相反，顶部的clusters基本上没有差异表达的基因组区域。差异表达的区域可能是由于拷贝数变异，并且各自的clusters可能代表肿瘤细胞。

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

UMAP plot of CNV profiles

我们可以将相同的clusters可视化为UMAP图。此外，infercnvpy.tl.cnv_score() 计算一个汇总分数，量化每个cluster的拷贝数变异量。它被简单地定义为每个cluster的 CNV 矩阵的绝对值的平均值。

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

UMAP 图由一大团“正常”细胞和几个具有不同 CNV 分布的较小cluster组成。除了由纤毛细胞组成的cluster“12”外，孤立的cluster都是上皮细胞。这些可能是肿瘤细胞，每个cluster代表一个单独的亚克隆。

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)

还可以在基于转录组学的 UMAP 图上可视化 CNV 分数和cluster。同样，可以看到存在属于不同 CNV cluster的上皮细胞subcluster，并且这些cluster往往具有最高的 CNV 分数。

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

基于这些观察，现在可以将细胞分配给“肿瘤”或“正常”。为此，在 adata.obs 中添加一个新列 cnv_status。

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)

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

The inferCNV method

本质上，这个包是 infercnv 的 Python 重新实现。通过使用 numpy、scipy 和稀疏矩阵，它的计算效率要高得多。

Computation steps

1、从所有细胞中减去参考基因表达。由于数据在对数空间中，这有效地计算了对数倍数变化。如果有多个类别的引用可用（即为 reference_cat 指定了多个细胞类型），则log fold change是“bounded”：

• 分别计算每个类别的平均基因表达。
• 在所有参考平均值的最小值和最大值范围内的值会收到 0 的对数倍数变化，因为它们不被视为与背景不同。
• 从小于所有参考平均值的最小值的值中减去该最小值。
• 从大于所有参考平均值的最大值的值中减去该最大值。

此过程避免了由于聚集基因区域（例如不同免疫细胞类型中的免疫球蛋白或 HLA 基因）的细胞类型特异性表达而调用假阳性 CNV 区域。

2、Clip the fold changes at -lfc_cap and +lfc_cap.

3、通过基因组位置平滑基因表达。 计算长度为 window_size 的运行窗口的平均值。仅计算每第 n 个窗口以节省时间和空间，其中 n = step。

4、通过从每个细胞中减去每个细胞的中位数，按细胞将平滑的基因表达居中。

5、执行噪声过滤。值 < dynamic_theshold * STDDEV 设置为 0，其中 STDDEV 是平滑基因表达的标准偏差

6、Smooth the final result using a median filter.

Preparing input data

• annadata.AnnData 对象应该已经被过滤为低质量的细胞。adata.X 需要进行规范化和对数转换。该方法应该对不同的归一化方法（scanpy.pp.normalize_total()、scran 等）相当稳健。
• 该方法需要一个“参考”值，与基因组区域的表达进行比较。如果数据集包含不同的细胞类型并且包括肿瘤细胞和正常细胞，则可以使用所有细胞的平均值作为参考。这是默认设置。
• 如果已经知道哪些细胞是“正常的”，可以提供从 adata.obs 到 reference_key 的列，其中包含注释。reference_cat 在reference_key 中指定一个或多个引用正常细胞的值。

生活很好，有你更好