空间转录组分析-04-单个样本分析

整个基本流程可以通过scanpy库实现。

环境准备

配置conda虚拟环境

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


conda create -n scanpy -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main python=3.10.14
conda activate scanpy
conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main scanpy==1.10.4
pip install squidpy==1.6.2
pip install dask==2024.11.2
pip install xarray==2024.10.0
pip install igraph
pip install leidenalg
pip install sopa --upgrade
pip install jupyter notebook

# 启动jupyter notebook
jupyter notebook --no-browser

导入相关python包

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21


import scanpy as sc
import squidpy as sq
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import os
sc.set_figure_params(facecolor='white',figsize=(8,8))
sc.settings.verbosity = 3
print(f"Scanpy version: {sc.__version__}")

'''
Scanpy version: 1.10.4
'''

print(sc.logging.print_header())

'''
scanpy==1.10.4 anndata==0.11.4 umap==0.5.3 numpy==1.26.4 scipy==1.15.2 pandas==2.2.3 scikit-learn==1.6.1 statsmodels==0.14.4 igraph==0.11.8 pynndescent==0.5.10
None
'''

读取数据

1
2
3
4
5


# 读取标准的10xVisium文件，即outs文件夹内的相关内容
adata = sq.read.visium('PC1/outs')

# 【重要！】基因名不唯一，应当使用.var_names_make_unique()使基因名唯一
adata.var_names_make_unique()

查看adata结构：

1
2
3
4
5
6
7
8
9


adata

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

质量控制

查看基因的表达情况

1、可视化高表达基因

1
2
3


# 可视化高表达基因，一般来说会有高表达基因中很多都是线粒体基因，需要去除
# 但是本例所使用的探针集中不含有线粒体基因，故无需去除
sc.pl.highest_expr_genes(adata,n_top=20) # 可视化前20个高表达的基因（所有spots的均值）

2、标记线粒体基因，进行质量控制统计

1
2
3
4
5
6


# 标记基因名以“MT-”开头的基因为线粒体基因，在adata.var中添加"mt"列
adata.var["mt"] = adata.var_names.str.startswith("MT-")

# 计算质量控制指标，并指定需要额外计算的质控变量为"mt"；inplace=True代表会直接在adata中进行修改
# 其实就是做描述性统计
sc.pp.calculate_qc_metrics(adata,qc_vars=["mt"],inplace=True)

3、对统计结果可视化

1
2
3


# 观察每个细胞的基因表达数、总counts数，线粒体百分比的小提琴图（分布情况）
# 画图参数jitter控制数据点在提琴图上的抖动（水平随机散布）程度。值越大点越分散。
sc.pl.violin(adata,['n_genes_by_counts','total_counts','pct_counts_mt'],jitter=0.4,multi_panel=True)

如上所属，该样本所用的探针集里不包含线粒体基因，故无法检测到线粒体基因。这不代表细胞内没有线粒体基因，只是检测不到而已。

可以看到该样本中每个细胞平均检测到大约6000个基因，12000个counts，未检测到线粒体基因。

接下来绘制散点图，可视化两个质控变量之间的关系：

1
2
3
4


fig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,4))
sc.pl.scatter(adata,x='total_counts',y='pct_counts_mt',show=False,ax=axes[0])
sc.pl.scatter(adata,x='total_counts',y='n_genes_by_counts',show=False,ax=axes[1])
plt.subplots_adjust(wspace=.3)

一般来说，total_counts和线粒体百分比没有显著的相关性，此外，n_genes_by_counts和total_counts应该呈现正相关的关系。

4、除了线粒体基因外，还可以自己定义筛选标准：

根据表达情况筛选spots和genes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20


# 筛选spots：
# 筛选总counts数在4000-35000的，以及线粒体基因占比＜20%的spots
sc.pp.filter_cells(adata,min_counts=5000)
sc.pp.filter_cells(adata,max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20]
print(f"#cells after MT filter:{adata.n_obs}")

# 筛选出至少在10个spots中出现的基因
sc.pp.filter_genes(adata,min_cells=10)

# 查看adata
adata

'''
AnnData object with n_obs × n_vars = 2939 × 15173
    obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'n_counts'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
    uns: 'spatial'
    obsm: 'spatial'
'''

寻找并保留高变基因

1、首先对每个细胞的counts进行归一化，使所有spots的总计数相同（默认 target_sum=1e4，即每个细胞的计数总和缩放至 10,000）。如果不归一化，高测序深度的spots会主导分析结果（如 PCA、聚类等）。

$$ X_{normalized}=\frac{X_i}{\sum{X_i}}\times target\_sum $$

其中$X_i$是每个spots中第i个基因的counts。

2、对归一化后的数据进行对数化（log transformation），即计算 log⁡(1+X)。

测序数据通常是高度偏态的即少数基因高表达，多数基因低表达。
对数化可以压缩数据的动态范围，使高表达基因的影响减小，同时保留低表达基因的信息。
加1可以避免出现$log\ 0$情况的出现

3、寻找方差最大的一批基因，这些基因通常是由分析意义的，因为它们在不同种类的细胞中表达有差异。只对高表达基因可以降低信噪比，减少后续数据处理的计算量。代码中计算方法一般指定未seurat法，即根据标准化的方差进行排序，选择变异度最大的基因。

1
2
3
4
5
6
7
8


# 归一化
sc.pp.normalize_total(adata,inplace=True)

# 对数化
sc.pp.log1p(adata)

# 计算前2000个高变基因，使用seurat法
sc.pp.highly_variable_genes(adata,flavor='seurat',n_top_genes=2000)

4、可视化，观察哪些基因被标记为了高变基因

1

sc.pl.highly_variable_genes(adata)

这里可以看出归一化带来的影响

5、保留高变基因，剔除其它基因

1
2
3
4
5
6


# 将adata的原始数据保存在adata的raw属性中，做一个备份。保存的内容包括adata.X和adata.var信息
# 如果想要调用原始数据，可以使用adata.raw.X和adata.raw.var来访问
adata.raw = adata

# 筛选出高变基因的列。支持额外指定参数，详情见官网
adata = adata[:,adata.var.highly_variable]

降维聚类

空间转录组的降维聚类分析与单细胞的降维聚类一致：

使用主成分分析将n维（这里是2000维）表达向量降维至低维（一般是50维，n_comps=50）。PCA 坐标（细胞 × PCs）存储在adata.obsm中，方差贡献率等统计信息存储在adata.uns[‘pca’]中。
基于 PCA 结果（adata.obsm['X_pca']）计算细胞距离。距离计算方式默认是欧式距离（metric='euclidean'），必要时可选择余弦距离（metric='cosine'）。然后构建K-近邻（KNN）图，用于后续聚类和降维。k默认值为15。
使用umap方法进行非线性降维，将高维数据从主成分空间映射到二维空间。umap坐标存储在adata.obsm[‘X_umap’]中。
基于KNN图进行降维聚类，使用Leiden法。Leiden算法是Louvain的改进版，更适合转录组数据。spots聚类结果存储在adata.obs[‘clusters’]中，聚类参数和模块度分数存储在adata.uns[’leiden’]中。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


# 主成分分析
sc.pp.pca(adata,svd_solver='arpack')

# 计算KNN图，默认k=15
sc.pp.neighbors(adata)

# umap降维
sc.tl.umap(adata)

# leiden聚类，resolution可以调节从0-1，数值越大，分出类就越多
sc.tl.leiden(adata,key_added='clusters',resolution=0.4)

进行可视化，spots被聚成了7类，并且可以可视化这六类spots在切片上的分布情况：

1
2
3
4


fig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,4))
sc.pl.umap(adata,color=["clusters"],wspace=0.4,ax=axes[0],show=False)
sc.pl.spatial(adata,img_key='hires',color=['clusters'],ax=axes[1],show=False)
plt.subplots_adjust(wspace=.3)

寻找Marker基因

使用sc.tl.rank_genes_groups()寻找Marker基因

通过t检验比较每个基因在该组与其它组之间的表达差异，筛选出每个leiden类中最具代表性的基因。最后基因按照p值大小（显著性）进行排序：

1
2
3
4
5


# 选出每个leiden cluster的Marker基因，使用student-t检验。可供选择的还包括wilcoxon或logreg。
sc.tl.rank_genes_groups(adata,'clusters',method='t-test')

# 可视化每组特异性高表达的前25个基因
sc.pl.rank_genes_groups(adata,n_genes=25,shareey=False)

例如，SPON2就是leiden1中的一个Marker基因：

Marker基因的可视化

Marker基因的可视化方法有很多，选择合适的即可：

1、基因的空间分布：

1
2
3
4
5
6


fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(17, 4))
sc.pl.umap(adata, color=["clusters"], wspace=0.4, ax=axes[0], show=False)
sc.pl.umap(adata, color=["SPON2"], wspace=0.4, ax=axes[1], show=False)
sc.pl.spatial(adata, img_key="hires", color=['clusters'], ax=axes[3], show=False)
plt.subplots_adjust(wspace=0.3)
sc.pl.spatial(adata, img_key="hires", color=['SPON2'], ax=axes[2], show=False, cmap='Spectral_r')

在umap图上和组织切片图上都可以看到SPON主要在第1类和第6类中表达

2、heatmap：

1

sc.pl.rank_genes_groups_heatmap(adata, groups="1", n_genes=10, groupby="clusters")

3、小提琴图：

观察指定分组的基因表达情况

1
2
3


from matplotlib.pyplot import rc_context
with rc_context({'figure.figsize':(9,1.5)}):
    sc.pl.rank_genes_groups_violin(adata,groups='1',n_genes=20,jitter=False)

观察指定基因在组间的表达情况

1
2


with rc_context({'figure.figsize':(9,3)}):
    sc.pl.violin(adata,['CD3D','SPON2'],groupby='cluster')

4、tracksplot:

1

sc.pl.rank_genes_groups_tracksplot(adata,n_genes=3)

指定基因集进行打分

scanpy库通过scanpy.tl.score_genes()函数，可以针对给定的基因集进行打分，即基因集中的基因在某个spots中表达量越多，则该spots的分数就越高。我们这里选择一些免疫细胞在分化过程中的标志基因作为基因集，并分别用四个基因集进行打分。如何选取基因集暂不在此介绍。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


# 选取合适的基因集
Immune_genes = [
    'CD3D','CD3E','CD2','CD4','CD8A',#T cell
    'CD79A','MZB1','MS4A1','CD79B',#B cell
    'FOXP3',"IL32",'TNFRSF18','TNFRSF4',#Treg
    'IL17A','IL17F','CD40LG',#Th17
    'S100A8','CXCL8','SOD2','NAMPT',#Neutrophil
    'SEPP1','C1QA','APOE','CD14','RNASE1',#Macrophage
    'TPSAB1','TPSB2','CPA3','HPGDS',#Mast
    'HLA-DRA','HLA-DPB1','CST3','HLA-DPA1',#mDC
    'PTGDS','SOX4','GZMB','IRF7',#pDC
    'IGHA1','IGHG1',"IGHG2",#Plasma
    'KLRF1','KLRD1','XCL2','XCL1'#NK
]
NK_genes = ['KLRF1','KLRD1','XCL2','XCL1']
T_cell_genes = ['CD3D','CD3E','CD2','CD4','CD8A']
B_cell_genes = ['CD79A','MZB1','MS4A1','CD79B']

# 对基因集进行打分
sc.tl.score_genes(adata, Immune_genes, score_name='Immune_score')
sc.tl.score_genes(adata, NK_genes, score_name='NK_score')
sc.tl.score_genes(adata, T_cell_genes, score_name='T_cell_score')
sc.tl.score_genes(adata, B_cell_genes, score_name='B_cell_score')

打完分后即可将评分可视化：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(17, 9))
sc.pl.umap(adata, color=["Immune_score"], wspace=0.4, ax=axes[0], show=False, cmap='Spectral_r', vmin=-0.5, vmax=0.5)
sc.pl.umap(adata, color=["NK_score"], wspace=0.4, ax=axes[1], show=False, cmap='Spectral_r', vmin=-0.5, vmax=0.5)
sc.pl.umap(adata, color=["T_cell_score"], wspace=0.4, ax=axes[2], show=False, cmap='Spectral_r', vmin=-0.5, vmax=0.5)
sc.pl.umap(adata, color=["B_cell_score"], wspace=0.4, ax=axes[3], show=False, cmap='Spectral_r', vmin=-0.5, vmax=0.5)
sc.pl.spatial(adata, img_key="hires", color=['Immune_score'], ax=axes[1], show=False, cmap='Spectral_r')
sc.pl.spatial(adata, img_key="hires", color=['NK_score'], ax=axes[2], show=False, cmap='Spectral_r')
sc.pl.spatial(adata, img_key="hires", color=['T_cell_score'], ax=axes[2], show=False, cmap='Spectral_r')
sc.pl.spatial(adata, img_key="hires", color=['B_cell_score'], ax=axes[2], show=False, cmap='Spectral_r')
plt.subplots_adjust(wspace=0.3)