Stereopy -- 华大空转分析

【记录使用，官方原文请参考：Quick Start (Square Bin) - Stereopy】

bin1 的纳米孔半径为 220nm；spot中心到中心的距离为 500nm；

image.png

安装

建议新建新环境，下载 requirement.txt 使用 pip 进行安装，尽量保持模块的版本与要求一致，以防后续分析中出现问题。安装完毕后使用示例数据运行下，查看是否与官网教程的结果一致。

import stereo as st
import warnings
warnings.filterwarnings('ignore')

此处展示的版本为 1.6.1，可通过st.__version__来查看

读取数据：

data_path = 'data/SS200000135TL_D1.tissue.gef'

# 查看数据信息(基因数、背景图的长宽等)
st.io.read_gef_info(data_path)

# 读取数据(查看细胞数x基因数)
data = st.io.read_gef(file_path=data_path, bin_size=50, is_sparse=True)
data

建议使用：可查看使用 gene_name_index 读取数据的spot数和基因数是否发生变化
data = st.io.read_gef(file_path=data_path, bin_size=50, is_sparse=True, gene_name_index=True)

gene_name_index=True：若 geneID 和 geneName 不一致，建议使用此参数读取 geneName，读取过程会较慢些；
is_sparse=True：设置矩阵为稀疏矩阵

质控

total_counts: 每个细胞的总数；
n_genes_by_counts：基因数；
pct_countss_mt：线粒体基因表达的百分比；.

data.tl.cal_qc()  # 统计各指标

# 绘制小提琴图
data.plt.violin(multi_panel=True)

keys='total_counts','n_genes_by_counts','pct_counts_mt']：默认展示这三个指标；
multi_panel=True：设置多个坐标轴；
out_path='test.pdf' 或 out_path='test.png'设置输出图片；

查看空间散点图：

data.plt.spatial_scatter()
# 可通过 keys 指定展示

查看指标间的关联性：

data.plt.genes_counts()

过滤：

data.tl.filter_cells：过滤细胞；
data.tl.filter_genes：过滤基因；
data.tl.filter_coordinates：过滤坐标；

但现有文献不知为何多数未进行过滤，或过滤很少部分；

data.tl.filter_cells(
  min_counts=20,  # 
  min_genes=3,  # 过滤掉基因数不足3的细胞；
  pct_counts_mt=5,  # 过滤掉线粒体基因表达比例大于5的细胞；
  inplace=True  # 替代
)

将过滤后的数据作为原始数据设置记忆点，可再后续使用时进行提取：

data.tl.raw_checkpoint()
data.tl.raw

标准化

常用的标准化方法：

normalize_total
log1p
scale
scTransform
quantile

这里运行 normalize_total 和 log1p的组合对基因表达矩阵进行归一化：

data.tl.normalize_total(target_sum=10000)
data.tl.log1p()

如果使用包含查找高度可变基因功能的data.tl.sctransform，则不需要运行data.tl.highly_variable_genes。在随后的data.tl.pca方法中，参数use_highly_genes必须设置为False。简而言之，是否使用高可变基因运行PCA取决于scTransform归一化中的filter_hvgs。

鉴定高变基因

（若是运行pca时不使用高变基因，为何要鉴定高变基因呢？）

data.tl.highly_variable_genes(
        min_mean=0.0125,
        max_mean=3,
        min_disp=0.5,
        n_top_genes=2000,
        res_key='highly_variable_genes'
        )
data.plt.highly_variable_genes(res_key='highly_variable_genes')

将基因换算成单位方差，剪辑值超过标准偏差10。zero_center=False 使用稀疏矩阵计算，减少运行所需内存。

data.tl.scale(max_value=10, zero_center=True)

PCA 主成分分析

data.tl.pca(
        use_highly_genes=False,
        n_pcs=30,
        res_key='pca'
        )

res_key 储存对应的结果，查看可使用 data.tl.key_record

选择 PCs：
根据绘图结果不应该选择 pcs=30 为何后续分析设置为 30 ？

data.plt.elbow(pca_res_key='pca')

邻近图 Neighborhood graph

和单细胞转录组 Seurat 分析中相同的道理

data.tl.neighbors(
        pca_res_key='pca',
        n_pcs=30,
        res_key='neighbors'
        )

# compute spatial neighbors
data.tl.spatial_neighbors(
        neighbors_res_key='neighbors',
        res_key='spatial_neighbors'
        )

UMAP

data.tl.umap(pca_res_key='pca', neighbors_res_key='neighbors', res_key='umap')

data.plt.umap(gene_names=['Atpif1', 'Tmsb4x'], res_key='umap')

建议设置两三个值看看其在 UMAP 上的分布，本人当初的 UMAP 图展示就出问题了，还在傻傻地往后分析呢。（注意 numba 与 umap 版本的问题）

分群 Clustering

常用的聚类方法，包括Leiden、Louvain和Phenograph

Leiden

data.tl.leiden(neighbors_res_key='neighbors', res_key='leiden')
data.plt.cluster_scatter(res_key='leiden')

data.tl.leiden(neighbors_res_key='spatial_neighbors', res_key='spatial_leiden')
data.plt.cluster_scatter(res_key='spatial_leiden')

Louvain

data.tl.louvain(neighbors_res_key='neighbors', res_key='louvain')
data.plt.cluster_scatter(res_key='louvain')

Phenograph

data.tl.phenograph(phenograph_k=30, pca_res_key='pca', res_key='phenograph')
data.plt.cluster_scatter(res_key='phenograph')

对比不同算法分析的结果：

以 leiden 为例，展示 UMAP图、单个cluster空间散点图、可交互式图，辅助后续注释：

data.plt.umap(res_key='umap', cluster_key='leiden')
data.plt.cluster_scatter(res_key='leiden', groups=['1', '2'])
data.plt.interact_cluster(res_key='leiden')

鉴定差异基因

若数据的差异基因具有特异性，可根据差异基因进行注释；若组织轮廓较明显的，可根据组织轮廓对应的分群进行注释。

注释：

annotation_dict = {
    '1':'a', '2':'b',
    '3':'c', '4':'d',
    '5':'e', '6':'f',
    '7':'g', '8':'h',
    '9':'i', '10':'j',
    '11':'k', '12': 'l',
    '13': 'm', '14': 'n',
    '15': 'o', '16': 'p',
    '17': 'q', '18': 'r',
    '19': 's', '20': 't',
    '21': 'u', '22': 'v',
    '23': 'w', '24': 'x',
    '25': 'y', '26': 'z'
    }
data.tl.annotation(
        annotation_information=annotation_dict,
        cluster_res_key='leiden',
        res_key='anno_leiden'
        )

保存数据

adata = st.io.stereo_to_anndata(data, output='SS200000135TL_D1_seurat.h5ad')

坑

自己输出的UMAP图和官方流程中输出的 UMAP 图超级不一致：

自己输出的 UMAP 可视化图.png

官方流程中输出的 UMAP 可视化图.png

初步怀疑是 umap 模块问题，然后怀疑 numba 模块问题，igraph 模块问题，后续怀疑 plotly 模块问题。依次按照安装中的 requirement.txt 进行检查。。。
最后解决：重新创建新环境，重新安装。

Stereopy -- 华大空转分析