社区所有版块导航
Python
python开源   Django   Python   DjangoApp   pycharm  
DATA
docker   Elasticsearch  
aigc
aigc   chatgpt  
WEB开发
linux   MongoDB   Redis   DATABASE   NGINX   其他Web框架   web工具   zookeeper   tornado   NoSql   Bootstrap   js   peewee   Git   bottle   IE   MQ   Jquery  
机器学习
机器学习算法  
Python88.com
反馈   公告   社区推广  
产品
短视频  
印度
印度  
Py学习  »  Python

用python版InferCNVpy加速运算

生信技能树 • 2 年前 • 907 次点击  

本质上,inferCNVpy这个包是InferCNV的python版重现。主要还是遵循R包版本的计算步骤,进行了少量修改。inferCNVpy通过使用numpy、scipy和稀疏矩阵,使其计算效率大大提高。inferCNVpy可以在Linux,Mac环境下运行。Windows下可参考:

inferCNVpy github见:https://github.com/icbi-lab/infercnvpy

官网认为这个方法仍处于实验阶段,但他们还是决定将其放到GitHub上,因为它可能会有用。但是需要大家小心看待结果:

首先是inferCNV软件的环境配置:

conda create -n infercnv
conda activate infercnv
pip install infercnvpy

umap-learn 0.5.2版本有一点bug,如果下面的UMAP图像显示异常,可以考虑用pip做一个降级

在Linux环境下利用Jupyter Notebook编译器运行下述python代码:

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

sc.settings.set_figure_params(dpi=80,figsize=(55))

Step1.载入示例数据

input的数据应该是已过滤后的单细胞数据,例如我们可以从Seurat流程进行QC过滤,然后导出数据。

示例数据是Smart-seq2生成的3000个单细胞:

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

可视化:

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

Step2.Running infercnv

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

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

3000个Smart-seq2单细胞约运行32.39秒。

与inferCNV R版教程一样,inferCNVpy可设置参考细胞:

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

Step3.可视化

绘制热图

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

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

CNV聚类和肿瘤细胞鉴定

为了对细胞进行聚类和注释,inferCNVpy镜像了scanpy工作流。以下函数与它们的scanpy对应函数完全一样,只是它们使用CNV剖面矩阵作为输入。使用这些函数,我们可以执行基于UMAP的聚类,并根据CNV数据生成UMAP图。基于这些clusters,我们可以注释肿瘤细胞和正常细胞:

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(22, figsize=(1111))
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(
    22, figsize=(1211), 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)

Step5.Classifying tumor cells

基于这些观察,我们现在可以将细胞分配给“tumor”或“normal”。为此,我们在 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(12, figsize=(125), 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", :])
- END -


Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/130647
 
907 次点击