单细胞转录组实战05: CellphoneDB细胞通讯及可视化
- 需要为cellphonedb单独准备一个环境
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.tsinghua.edu.cn/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.tsinghua.edu.cn/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 file
标准化的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版的,因此会报错...
- Metadata file
第一列: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.readthedocs.io/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.readthedocs.io/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
相关文章
- 在 Go 里用 CGO?这 7 个问题你要关注!
- 9款优秀的去中心化通讯软件 Matrix 的客户端
- 求职数据分析,项目经验该怎么写
- 在OKR中,我看到了数据驱动业务的未来
- 火山引擎云原生大数据在金融行业的实践
- OpenHarmony富设备移植指南(二)—从postmarketOS获取移植资源
- 《数据成熟度指数》报告:64%的企业领袖认为大多数员工“不懂数据”
- OpenHarmony 小型系统兼容性测试指南
- 肯睿中国(Cloudera):2023年企业数字战略三大趋势预测
- 适用于 Linux 的十大命令行游戏
- GNOME 截图工具的新旧截图方式
- System76 即将推出的 COSMIC 桌面正在酝酿大变化
- 2GB 内存 8GB 存储即可流畅运行,Windows 11 极致精简版系统 Tiny11 发布
- 迎接 ecode:一个即将推出的具有全新图形用户界面框架的现代、轻量级代码编辑器
- loongarch架构介绍(三)—地址翻译
- Go 语言怎么解决编译器错误“err is shadowed during return”?
- 敏捷:可能被开发人员遗忘的部分
- Denodo预测2023年数据管理和分析的未来
- 利用数据推动可持续发展
- 在 Vue3 中实现 React 原生 Hooks(useState、useEffect),深入理解 React Hooks 的