单细胞空间转录组分析流程学习python版(三)

admin管理员组

文章数量:1794759

单细胞空间转录组分析流程学习python版(三)

既往推文：

单细胞空间转录组分析流程学习(一)：

单细胞空间转录组分析流程学习(二)：

分析步骤

1.配置一个环境/安装必要软件

代码语言：javascript代码运行次数：0运行复制

conda creadt -n scRNA_spatial python=3.9
conda install -c conda-forge scanpy pandas matplotlib seaborn

2.导入

代码语言：javascript代码运行次数：0运行复制

import scanpy as sc
import anndata as an
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sc.logging.print_versions()
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3

3.读入文件

代码语言：javascript代码运行次数：0运行复制

adata = sc.read_visium(path = "./RawData/GSE6716963",
                       count_file="GSM6716963_19G081_filtered_feature_bc_matrix.h5")

4.check

代码语言：javascript代码运行次数：0运行复制

adata

#AnnData object with n_obs × n_vars = 2887 × 17943
#    obs: 'in_tissue', 'array_row', 'array_col'
#    var: 'gene_ids', 'feature_types', 'genome'
#    uns: 'spatial'
#    obsm: 'spatial'

再次学习一下这里的参数含义：

n_obs × n_vars：2887 个观测点（通常是细胞或空间点）和 17943 个基因。

obs：存储每个观测点的元数据，如是否位于组织中以及在空间阵列上的坐标。

var：存储每个基因的元数据，如基因 ID、基因类型和基因组信息。

uns：存放未结构化数据，如与空间信息相关的元数据。

obsm：存储每个观测点的多维数据（如空间坐标或降维结果）。

代码语言：javascript代码运行次数：0运行复制

adata.obsm['spatial']

# array([[3426, 1363],
#       [3980, 4005],
#       [1430, 1658],
#       ...,
#       [4028, 3652],
#       [3196, 3765],
#       [1068, 3040]])

这表示每个观测点在组织切片中的位置。每一行是一个观测点的坐标，其中第一个数是 x 坐标，第二个数是 y 坐标。

代码语言：javascript代码运行次数：0运行复制

adata.obsm['spatial'].shape

# (2887, 2)

2887：这个矩阵有 2887 行，表示有 2887 个观测点（细胞或空间点）。

2：这个矩阵有 2 列，表示每个观测点有 2 个空间坐标（通常是二维的 x 和 y 坐标）。

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial'].keys()

# dict_keys(['19G081'])

adata.uns['spatial'].keys() 返回了字典中键的名称 '19G081'，这是一个样本或切片的标识符

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial']['19G081'].keys()

# dict_keys(['images', 'scalefactors', 'metadata'])

表示19G081字典下面储存着三项内容：images/scalefactors/metadata

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial']['19G081']['images'].keys()

# dict_keys(['hires', 'lowres'])

iamges下面有高低分辨率两种图片

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial']['19G081']['images']['hires'].shape
# (1878, 2000, 3)

adata.uns['spatial']['19G081']['images']['lowres'].shape
# (564, 600, 3)

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial']['19G081']['scalefactors']

#{'tissue_hires_scalef': 0.36231884,
# 'tissue_lowres_scalef': 0.10869565,
# 'fiducial_diameter_fullres': 67.115571765934,
# 'spot_diameter_fullres': 41.547734902721054}

tissue_hires_scalef 和 tissue_lowres_scalef：用于将图像的像素坐标转换为实际组织坐标，分别对应高分辨率和低分辨率图像。

fiducial_diameter_fullres：在全分辨率图像中，标识性点（fiducial marker）的直径，通常用于图像校准。

spot_diameter_fullres：在全分辨率图像中，空间点（spot）的直径，表示测量基因表达的区域的实际大小。

代码语言：javascript代码运行次数：0运行复制

adata.uns['spatial']['19G081']['metadata']

#{'chemistry_description': "Spatial 3' v1",
# 'software_version': 'spaceranger-1.3.0'}

代码语言：javascript代码运行次数：0运行复制

adata.obs.head()

这里跟seurat中的很相似，tissue代表不同的组织切片/row和col代表位置坐标

5.数据质量控制

代码语言：javascript代码运行次数：0运行复制

adata.var_names_make_unique()
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

adata.var_names_make_unique() 确保基因名称唯一。

adata.var["mt"] = adata.var_names.str.startswith("MT-") 标记线粒体基因。

sc.pp.calculate_qc_metrics 计算质量控制指标，包括线粒体基因表达的比例，用于后续数据过滤和质量评估。

代码语言：javascript代码运行次数：0运行复制

adata.obs.head()

6.增加slice列

代码语言：javascript代码运行次数：0运行复制

adata.obs['slice'] = ['slice' + str(i) for i in adata.obs['in_tissue']]
adata.obs.head()

7.绘制一下count的小提琴图

代码语言：javascript代码运行次数：0运行复制

sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], groupby="slice", jitter=0.4, multi_panel=True)

8.过滤数据

也有老师说可以不过滤hhhh

代码语言：javascript代码运行次数：0运行复制

sc.pp.filter_cells(adata, min_counts=5000)
sc.pp.filter_cells(adata, max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20]
sc.pp.filter_genes(adata, min_cells=10)

9.归一化/高变基因

代码语言：javascript代码运行次数：0运行复制

sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

10.降维聚类

代码语言：javascript代码运行次数：0运行复制

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters")

11.UMPA绘图

代码语言：javascript代码运行次数：0运行复制

plt.rcParams["figure.figsize"] = (4, 4)
sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

12.空转绘图

代码语言：javascript代码运行次数：0运行复制

plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts", "clusters"], alpha_img=1)

代码语言：javascript代码运行次数：0运行复制

sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["1","2", "3"], size = 1.3)

特定cluster绘制

13.裁剪HE图片

代码语言：javascript代码运行次数：0运行复制

adata.obsm['spatial'][:, 0].min(), adata.obsm['spatial'][:, 0].max()
# (np.int64(913), np.int64(4338))

adata.obsm['spatial'][:, 1].min(), adata.obsm['spatial'][:, 1].max()
# (np.int64(1751), np.int64(4177))

# [left, right, top, bottom]
sc.pl.spatial(adata, img_key="hires", 
              color="clusters", groups=["1","2","3"], size = 1.3, 
              crop_coord=[500, 2000, 1000, 3000])

随便剪裁了一下

14.clusters差异基因分析并绘图

计算标记基因并绘制cluster5的差异基因在不同cluster中表达情况的热图。

代码语言：javascript代码运行次数：0运行复制

sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
sc.pl.rank_genes_groups_heatmap(adata, groups="5", n_genes=10, groupby="clusters")

15.cluster和基因表达展示

代码语言：javascript代码运行次数：0运行复制

sc.pl.spatial(adata, img_key="hires", color=["clusters", "TP53"])
sc.pl.spatial(adata, img_key="hires", color=["COL1A2", "SYPL1"], alpha=0.7)

参考资料：

1. 10XGENOMICS：
2. scanpy空转：.html

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟

- END -

本文标签： 单细胞空间转录组分析流程学习python版(三)