CellphoneDB及可视化

生信技能树jimmy

发布于 2023-02-16 21:14:17

1.5K00

代码可运行

文章被收录于专栏：单细胞天地单细胞天地

运行总次数：0

代码可运行

CellphoneDB及可视化

CellphoneDB介绍以及结果怎么看可以参考[https://zhuanlan.zhihu.com/p/446055519]，这里就不再赘述了，直接开始跑代码。

准备环境

需要为cellphonedb单独准备一个环境

micromamba和conda用法一样，只是速度更快，可以把micromamba换成conda，另外conda一次装十几个包是会报错的。

micromamba create -n cpdb; micromamba activate cpdb
micromamba install -y python=3.7
pip install -U cellphonedb
cellphonedb database download ## V4

单细胞分析环境

micromamba create -n SC; micromamba activate SC
micromamba install -y -c conda-forge ipywidgets pandas numpy seaborn pytables matplotlib ipykernel scanpy python-igraph leidenalg scvi-tools scikit-misc python-annoy pyarrow
micromamba install -y -c bioconda harmonypy
pip install ktplotspy ## 0.1.6
pip install datatable ## 1.0.0

准备输入文件

import numpy as np
import pandas as pd
import scanpy as sc
import datatable as dt

Count file

标准化的count数据，行是人基因名，列是细胞名

## adata是之前处理好的数据
X=pd.DataFrame(adata.X.toarray().T,columns=adata.obs_names)
Symbol=dt.Frame({'Symbol':adata.var_names.values})
X = dt.cbind([Symbol,dt.Frame(X)])
X.to_csv('X.csv') ## 使用datatable加速导出

注：其实cellphonedb 3.1.0支持anndata格式的count文件，但是cellphonedb依赖的anndata是0.7版本的，而scanpy导出的anndata是0.8版的，因此会报错...

Metadata file

第一列：Cell，细胞名

第二列：cell_type，细胞类型

obs = pd.DataFrame({'Cell':adata.obs.index,'cell_type':adata.obs.CellType})
obs.to_csv('obs.csv',index=False)

运行cellphonedb

X=X.csv
obs=obs.csv
od=.
n_jobs=8

cellphonedb method statistical_analysis \
  --counts-data hgnc_symbol \
  --output-path ${od} \
  --threshold 0.1 \
  --threads ${n_jobs} \
  --iterations 1000 \
  --output-format csv \
  ${obs} \
  ${X}

结果可视化

cellphonedb自带的画图太弱了，直接使用ktplotspy绘图。

ktplotspy包括三个绘图函数:

plot_cpdb_heatmap画热图，
plot_cpdb画点图，
plot_cpdb_chord画弦图。

三个函数的API可以参考[https://ktplotspy.readthedocs.io/en/latest/api.html]

## 导库
import pandas as pd
import scanpy as sc
import datatable as dt
import matplotlib.pyplot as plt
import ktplotspy as kpy

准备输入文件

需要用于cellphonedb的adata文件，以及cellphonedb的输出文件。

## 用于cellphonedb的adata文件
adata = sc.read_h5ad('adata.h5')

## 读取cellphonedb的输出文件
means = pd.read_csv('means.csv')
pvals = pd.read_csv('pvalues.csv')
decon = pd.read_csv('deconvoluted.csv')

热图

kpy.plot_cpdb_heatmap(
        adata=adata,
        pvals=pvals,
        celltype_key="CellType", ## adata.obs中细胞类型那一列列名
        figsize = (5,5),
        title = "Sum of significant interactions"
        ## symmetrical = True ## 显得对称
    );

点图

kpy.plot_cpdb(
        adata=adata,
        cell_type1="B",
        cell_type2=".", ## .表示全部细胞类型
        means=means,
        pvals=pvals,
        celltype_key="CellType",
        genes=["PTPRC", "CD40", "CLEC2D"],
        figsize = (10,5),
        ## return_table = True, ## 返回dataframe，不返回图片
        ## default_style = False, ## squidpy风格
        highlight_size = 2, ## 显著点 外环的红色的大小
        title = "interacting interactions!"
    )

squidpy风格

基因家族

内置的基因家族包括，chemokines, th1, th2, th17, treg, costimulatory, coinhibitory

kpy.plot_cpdb(
        adata=adata,
        cell_type1="B",
        cell_type2="T|Mac",
        means=means,
        pvals=pvals,
        celltype_key="CellType",
        gene_family = "chemokines",
        ## custom_gene_family=dict(familyA=gene_list),自定义基因家族 
        figsize = (5,4)
    )

弦图

可以自定义face和edge的颜色，详细参考:

https://ktplotspy.readthedocs.io/en/latest/notebooks/tutorial.html

kpy.plot_cpdb_chord(
        adata=adata,
        cell_type1="B",
        cell_type2=".",
        means=means,
        pvals=pvals,
        deconvoluted=decon,
        celltype_key="CellType",
        genes=["PTPRC", "CD40", "CLEC2D"],
        edge_cmap=plt.cm.coolwarm, ## 弦的颜色
        figsize=(6,6)
    )