准备环境
micromamba和conda用法一样，只是速度更快，可以把micromamba换成conda，另外conda一次装十几个包是会报错的。
micromamba create -n cpdb; micromamba activate cpdb
micromamba install -y python=3.7
pip install -U cellphonedb -i https://pypi.tuna./simple
cellphonedb database download
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 git+https://github.com/h2oai/datatable
pip install ktplotspy -i https://pypi.tuna./simple
准备输入文件
import numpy as np
import pandas as pd
import scanpy as sc
import datatable as dt
import matplotlib.pyplot as plt
import ktplotspy as kpy
标准化的count数据，行是人基因名，列是细胞名
normalized count data （row: human Symbol X col: cell name）
adata = sc.read_h5ad('output/03.inferCNV/adata.h5')
adata_raw = adata.raw.to_adata()
X=pd.DataFrame(adata_raw.X.toarray().T,columns=adata_raw.obs_names)
Symbol=dt.Frame({'Symbol':adata_raw.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版的，因此会报错...
第一列：Cell，细胞名
第二列：cell_type，细胞类型
obs = pd.DataFrame({'Cell':adata.obs.index,'cell_type':adata.obs.CellTypeS3})
obs.to_csv(OUTPUT_DIR+'/obs.csv',index=False)
运行cellphonedb
#>>>run_cellphone.sh>>>
mtx_path=output/04.CellPhoneDB/X.csv
obs_path=output/04.CellPhoneDB/obs.csv
output_path=output/04.CellPhoneDB
n_jobs=12
cellphonedb method statistical_analysis \
--counts-data hgnc_symbol \
--output-path ${output_path} \
--threshold 0.1 \
--threads ${n_jobs} \
--iterations 1000 \
--output-format csv \
${obs_path} \
${mtx_path}
#<<<run_cellphone.sh<<<
nohup bash ~/Project/SC10X/src/run_cellphone.sh &> output/04.CellPhoneDB/run_cellphone.sh.log
结果可视化
cellphonedb自带的画图太弱了，直接使用ktplotspy绘图。
ktplotspy包括三个绘图函数，plot_cpdb_heatmap画热图，plot_cpdb画点图，plot_cpdb_chord画弦图。
三个函数的API可以参考https://ktplotspy./en/latest/api.html
准备输入文件
需要用于cellphonedb的adata文件，以及cellphonedb的输出文件。
# 读取cellphonedb的输出文件
means = pd.read_csv(OUTPUT_DIR+'/means.csv')
pvals = pd.read_csv(OUTPUT_DIR+'/pvalues.csv')
decon = pd.read_csv(OUTPUT_DIR+'/deconvoluted.csv')
celltype_key="CellTypeS3"
热图
一般情况下，如果两类细胞或两组 Cluster 之间有较多的互动交流，那么它们间的基因关系也应该更为丰富。因此选用两类细胞间存在的互作数目作为评估细胞类型之间交流密切的依据。CellPhoneDB 统计 Cluster 间的蛋白互作数目作热图。横纵坐标表示分群 ID或细胞类型，颜色深蓝色到紫红色代表互作数由低到高。
kpy.plot_cpdb_heatmap(
adata=adata_raw,
pvals=pvals,
celltype_key=celltype_key,
figsize = (5,5),
title = "Sum of significant interactions"
# ,symmetrical = True
);
点图
点图可以展示蛋白互作关系在细胞类型中的平均表达水平以及显著性。如下图横坐标是细胞类型互作，纵坐标是蛋白互作，点越大表示 p 值越小，颜色代表平均表达量。
kpy.plot_cpdb(
adata=adata_raw,
cell_type1="B",
cell_type2=".", # this means all cell-types
means=means,
pvals=pvals,
max_size=8,
max_highlight_size=2,
celltype_key=celltype_key,
genes=["PTPRC", "TNFSF13"],
figsize = (10,5),
title = "interacting interactions!"
)
内置的基因家族包括，chemokines, th1, th2, th17, treg, costimulatory, coinhibitory
kpy.plot_cpdb(
adata=adata_raw,
cell_type1="B",
cell_type2="T|Mac",
means=means,
pvals=pvals,
max_size=8,
max_highlight_size=2,
celltype_key=celltype_key,
gene_family = "chemokines",
# custom_gene_family=dict(familyA=genelist),自定义基因家族
figsize = (4,4)
)
弦图
配体受体互作 Circos 图：其中最外圈代表的是配受体互作对所包含的所有 Cluster，箭头代表指向关系（x弦越宽说明表达两个基因的平均表达乘积越高）。
可以自定义face和edge的颜色，详细参考https://ktplotspy./en/latest/notebooks/tutorial.html
kpy.plot_cpdb_chord(
adata=adata,
cell_type1="B",
cell_type2=".",
means=means,
pvals=pvals,
deconvoluted=decon,
scale_lw=100,
celltype_key=celltype_key,
genes=["PTPRC", "CD40", "CLEC2D"],
# edge_cmap=plt.cm.coolwarm, # 弦的颜色
figsize=(6,6)
);
Reference
https://zhuanlan.zhihu.com/p/446055519
https://zhuanlan.zhihu.com/p/490179003
https://www.jianshu.com/p/b3f3a03afa53
