1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

# 导入相关库
import os
import scanpy as sc
import anndata as an
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import scanorama

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

# 读取数据，即每个样本的outs/文件夹，这里原始数据已经被读入并保存为为h5ad文件。故直接读取。
adata_spatial_anterior = sc.read('adata_spatial_anterior.h5ad')
adata_spatial_posterior = sc.read('adata_spatial_posterior.h5ad')


# 使基因名唯一
adata_spatial_anterior.var_names_make_unique()
adata_spatial_posterior.var_names_make_unique()

# 首先分别进行标准化与取对数（而不是合并之后再标准化）
sc.pp.normalize_total(adata_spatial_anterior,inplace=True)
sc.pp.log1p(adata_spatial_anterior)
sc.pp.normalize_total(adata_spatial_posterior,inplace=True)
sc.pp.log1p(adata_spatial_posterior)

# 【关键步骤！】将两个样本合并。
# 两个样本的var_names（基因名）集合虽无需完全相同（程序会自动地将确实地基因用NaN来标记），
# 但是为了防止下游分析出现问题还是要预先对二者基因集取交集.
# 此外，var_names的顺序无需一致，程序会自动对齐.
# 这里两个对象的基因集已经是取过交集之后的结果，故在此不另作取交集操作。
adatas = [adata_spatial_anterior,adata_spatial_posterior]
adata = sc.concat(
    adatas,label='batch',
    uns_merge='unique',keys=[
        k
        for d in [
            adatas[0].uns['spatial'],
            adatas[1].uns['spatial'],
        ]
        for k,v in d.items()
    ],
    index_unique='-'
)

# 归一化，对数化筛选出高变基因
sc.pp.highly_variable_genes(adata,flavor='seurat',n_top_genes=2000)

# 备份，然后删去非高变基因
adata.raw = adata
adata = adata[:,adata.var.highly_variable]

# pca降维，构建邻域图，构建umap图，leiden聚类
sc.pp.pca(adata,svd_solver = 'arpack')
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata,key_added='clusters',resolution=1)

#### 可视化 ####
# umap
sc.pl.umap(adata,color=['clusters','batch'],palette=sc.pl.palettes.default_20)

# embedding_density
sc.tl.embedding_density(adata,groupby='batch')
sc.pl.embedding_density(adata,groupby='batch')

# spatial
clusters_colors = dict(
    zip([str(i) for i in range(18)], adata.uns["clusters_colors"]))

fig, axs = plt.subplots(1, 2, figsize=(12, 8))

for i, batch in enumerate(
    ["V1_Mouse_Brain_Sagittal_Anterior", "V1_Mouse_Brain_Sagittal_Posterior"]
):
    ad = adata[adata.obs.batch == batch, :].copy()
    sc.pl.spatial(
        ad,img_key="hires",
        library_id = batch,
        color="clusters",size=1.5,palette=[
            v
            for k, v in clusters_colors.items()
            if k in ad.obs.clusters.unique().tolist()],
        legend_loc=None,
        show=False,
        ax=axs[i], )

plt.tight_layout()

1
2

# 安装BBKNN
pip install bbknn

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

import bbknn

# 根据bbknn原理，该算法返回的是一个t-临近图，相当于是取代了原本的`sc.pp.neighbors(adata)`的作用。
# 当然，在此之前还是需要pca降维。
# 结果存储在`adata.uns['neighbors']`；`adata.obsp['distances']`和`adata['connectivities']`中。
bbknn.bbknn(adata)

# 根据bbknn生成的t临近图，重新计算umap图坐标、leiden聚类结果
sc.tl.umap(adata) #根据新的TNN结果计算umap图
sc.tl.leiden(adata,key_added='clusters',resolution=1) # 根据根据新的TNN结果进行leiden聚类

#### 可视化 ####
# umap
sc.pl.umap(adata,color=['clusters','batch'],palette=sc.pl.palettes.default_20)

# embedding_density
sc.tl.embedding_density(adata,groupby='batch')
sc.pl.embedding_density(adata,groupby='batch')

# spatial
clusters_colors = dict(
    zip([str(i) for i in range(18)], adata.uns["clusters_colors"]))

fig, axs = plt.subplots(1, 2, figsize=(12, 8))

for i, batch in enumerate(
    ["V1_Mouse_Brain_Sagittal_Anterior", "V1_Mouse_Brain_Sagittal_Posterior"]
):
    ad = adata[adata.obs.batch == batch, :].copy()
    sc.pl.spatial(
        ad,img_key="hires",
        library_id = batch,
        color="clusters",size=1.5,palette=[
            v
            for k, v in clusters_colors.items()
            if k in ad.obs.clusters.unique().tolist()],
        legend_loc=None,
        show=False,
        ax=axs[i], )

plt.tight_layout()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

# 合并数据
adatas = [adata_spatial_anterior,adata_spatial_posterior]
adata = sc.concat(
    adatas,label='batch',
    uns_merge='unique',keys=[
        k
        for d in [
            adatas[0].uns['spatial'],
            adatas[1].uns['spatial'],
        ]
        for k,v in d.items()
    ],
    index_unique='-'
)

# 筛选高变基因
sc.pp.highly_variable_genes(adata, flavor='seurat', n_top_genes=2000)

# 备份原始数据，并取高变基因数据子集
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

# 【关键步骤！】
sc.pp.combat(adata, key='batch')

# 后续降维与聚类
sc.pp.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added='clusters', resolution=1)

#### 可视化 ####
# umap
sc.pl.umap(adata,color=['clusters','batch'],palette=sc.pl.palettes.default_20)

# embedding_density
sc.tl.embedding_density(adata,groupby='batch')
sc.pl.embedding_density(adata,groupby='batch')

# spatial
clusters_colors = dict(
    zip([str(i) for i in range(18)], adata.uns["clusters_colors"]))

fig, axs = plt.subplots(1, 2, figsize=(12, 8))

for i, batch in enumerate(
    ["V1_Mouse_Brain_Sagittal_Anterior", "V1_Mouse_Brain_Sagittal_Posterior"]
):
    ad = adata[adata.obs.batch == batch, :].copy()
    sc.pl.spatial(
        ad,img_key="hires",
        library_id = batch,
        color="clusters",size=1.5,palette=[
            v
            for k, v in clusters_colors.items()
            if k in ad.obs.clusters.unique().tolist()],
        legend_loc=None,
        show=False,
        ax=axs[i], )

plt.tight_layout()

1
2

# 安装scanpy的扩展
pip install harmonypy

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

adatas = [adata_spatial_anterior,adata_spatial_posterior]
adata = sc.concat(
    adatas,label='batch',
    uns_merge='unique',keys=[
        k
        for d in [
            adatas[0].uns['spatial'],
            adatas[1].uns['spatial'],
        ]
        for k,v in d.items()
    ],
    index_unique='-'
)

# 归一化，对数化筛选出高变基因
## sc.pp.normalize_total(adata,inplace=True)
## sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata,flavor='seurat',n_top_genes=2000)

# 备份，然后删去非高变基因
adata.raw = adata
adata = adata[:,adata.var.highly_variable]

sc.pp.pca(adata,svd_solver = 'arpack')
# 【关键步骤！】
# 使用sc.external.pp.harmony_integrate()进行去批次，指定批次信息字段为adata.obs['batch']
sc.external.pp.harmony_integrate(adata,key='batch')
sc.pp.neighbors(adata,use_rep='X_pca_harmony') # 以X_pca_harmony为依据计算KNN
sc.tl.umap(adata) #根据新的KNN结果计算umap图
sc.tl.leiden(adata,key_added='clusters',resolution=1) # 根据根据新的KNN结果进行leiden聚类

#### 可视化 ####
# umap
sc.pl.umap(adata,color=['clusters','batch'],palette=sc.pl.palettes.default_20)

# embedding_density
sc.tl.embedding_density(adata,groupby='batch')
sc.pl.embedding_density(adata,groupby='batch')

# spatial
clusters_colors = dict(
    zip([str(i) for i in range(18)], adata.uns["clusters_colors"]))

fig, axs = plt.subplots(1, 2, figsize=(12, 8))

for i, batch in enumerate(
    ["V1_Mouse_Brain_Sagittal_Anterior", "V1_Mouse_Brain_Sagittal_Posterior"]
):
    ad = adata[adata.obs.batch == batch, :].copy()
    sc.pl.spatial(
        ad,img_key="hires",
        library_id = batch,
        color="clusters",size=1.5,palette=[
            v
            for k, v in clusters_colors.items()
            if k in ad.obs.clusters.unique().tolist()],
        legend_loc=None,
        show=False,
        ax=axs[i], )

plt.tight_layout()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

# 【关键步骤！】将所有的空间转录结果构建列表，再运行scanorma进行矫正，对批次效应进行去除
# 与BBKNN一样，也是对PCA坐标进行操作。不过程序内部就会进行主成分分析，原先的sc.pp.pca()不是必要的
adata_spatial_anterior = sc.read('data/k8/adata_spatial_anterior.h5ad')
adata_spatial_posterior = sc.read('data/k8/adata_spatial_posterior.h5ad')

# 使基因名唯一
adata_spatial_anterior.var_names_make_unique()
adata_spatial_posterior.var_names_make_unique()

###### 首先分别进行标准化与取对数？
sc.pp.normalize_total(adata_spatial_anterior,inplace=True)
sc.pp.log1p(adata_spatial_anterior)

sc.pp.normalize_total(adata_spatial_posterior,inplace=True)
sc.pp.log1p(adata_spatial_posterior)



adatas = [adata_spatial_anterior,adata_spatial_posterior,]
adatas_cor = scanorama.correct_scanpy(adatas,return_dimred=True)

# 将两个数据进行合并
adata = sc.concat(
    adatas_cor,label='batch',
    uns_merge='unique',keys=[
        k
        for d in [
            adatas_cor[0].uns['spatial'],
            adatas_cor[1].uns['spatial'],
        ]
        for k,v in d.items()
    ],
    index_unique='-'
)

# 降维聚类
sc.pp.neighbors(adata,use_rep='X_scanorama')
sc.tl.umap(adata)
sc.tl.leiden(adata,key_added='clusters',resolution = 0.6)

# 可视化
sc.pl.umap(adata,color=['clusters','batch'],palette=sc.pl.palettes.default_20)
sc.tl.embedding_density(adata,groupby='batch')
sc.pl.embedding_density(adata,groupby='batch')

# 最终在切片图片上标注上对应的细胞类型。想要达到的效果即不同切片上的相同细胞类型的颜色需要一致
clusters_colors = dict(
    zip([str(i) for i in range(18)], adata.uns["clusters_colors"]))

fig, axs = plt.subplots(1, 2, figsize=(12,8))

for i, batch in enumerate(
    ["V1_Mouse_Brain_Sagittal_Anterior", "V1_Mouse_Brain_Sagittal_Posterior"]
):
    ad = adata[adata.obs.batch == batch, :].copy()
    sc.pl.spatial(
        ad,img_key="hires",
        library_id = batch,
        color="clusters",size=1.5,palette=[
            v
            for k, v in clusters_colors.items()
            if k in ad.obs.clusters.unique().tolist()],
        legend_loc=None,
        show=False,
        ax=axs[i], )

plt.tight_layout()