ixxmu
commented
7 months ago Scanpy生态教程12：推断CNV by 生信会客厅
设置分析环境
import scanpy as sc
import infercnvpy as cnv
import matplotlib.pyplot as plt
import os
import warningswarnings.simplefilter(action="ignore")
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=80, figsize=(5, 4.5))projpath = 'D:\\公众号\\2024scanpy'
chapter = "10_inferCNV"os.chdir(projpath)
if not os.path.exists(chapter):
os.makedirs(chapter)
os.chdir(chapter)推断CNV
infercnvpy.io.genomic_position_from_gtf可以帮助我们从gtf文件中提取信息并保存到adata对象。adata = cnv.datasets.maynard2020_3k()
adata
AnnData object with n_obs × n_vars = 3000 × 55556obs: 'age', 'sex', 'sample', 'patient', 'cell_type'var: 'ensg', 'mito', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_counts', 'chromosome', 'start', 'end', 'gene_id', 'gene_name'uns: 'cell_type_colors', 'neighbors', 'umap'obsm: 'X_scVI', 'X_umap'obsp: 'connectivities', 'distances'
# 检查数据类型，此数据集是log1p数据
print(f"{adata.X.min()}\n{adata.X.max()}")
adata.to_df()0.013.5234375
sc.pl.umap(adata, color="cell_type")# We provide all immune cell types as "normal cells".
cnv.tl.infercnv(
adata, n_jobs=6,
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,
)100%|██████████| 1/1 [00:42<00:00, 42.87s/it]# cnv矩阵保存在obsm['X_cnv']，行为细胞列为基因组分箱窗口（与基因平滑处理时的窗口不同），每个窗口的宽度大约5M
adata.obsm['X_cnv'].shape(3000, 5314)cnv.pl.chromosome_heatmap(adata, groupby="cell_type")CNV聚类
cnv.tl.pca(adata)
cnv.pp.neighbors(adata)
cnv.tl.leiden(adata)computing PCAwith n_comps=50finished (0:00:01)computing neighborsfinished: added to `.uns['cnv_neighbors']``.obsp['cnv_neighbors_distances']`, distances for each pair of neighbors`.obsp['cnv_neighbors_connectivities']`, weighted adjacency matrix (0:00:13)running Leiden clusteringfinished: found 18 clusters and added'cnv_leiden', the cluster labels (adata.obs, categorical) (0:00:00)
cnv.pl.chromosome_heatmap(adata, groupby="cnv_leiden", dendrogram=True)CNV降维与评分
cnv.tl.umap(adata)
cnv.tl.cnv_score(adata)computing UMAPfinished: added'X_umap', UMAP coordinates (adata.obsm) (0:00:07)
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图上展示数据
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 11), gridspec_kw={"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)肿瘤细胞分类
adata.obs["cnv_status"] = "normal"
adata.obs.loc[
adata.obs["cnv_leiden"].isin(["10", "16", "13", "8", "12", "17", "1", "14", "11"]), "cnv_status"] = "tumor"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), gridspec_kw={"wspace": 0.5})
cnv.pl.umap(adata, color="cnv_status", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_status", ax=ax2)... storing 'cnv_status' as categorical# 单独展示恶心细胞的cnv热图
cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "tumor", :])# 单独展示正常细胞的cnv热图
cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "normal", :])交流合作 Kinesin专业从事单细胞与空转数据的分析与培训，有此需求的朋友欢迎加微信洽谈合作。本公众号交流群向科研院校师生和公司生信人员开放，希望入群的朋友也可加我微信申请。
https://mp.weixin.qq.com/s/byA6MmhvEs8c7oDU-DMzVg