空间转录组

启动服务器端jupyter notebook的流程：

1
jupyter notebook --ip=0.0.0.0 --port=8888 --no-browser

空间转录组的流程分析

  1. 聚类
    • 不同的功能区域
      • 对比单细胞转录组分析的聚类：细胞类型
  2. 降维
    • 主要使用空间信息进行进行降维，如UMAP、t-SNE，同时考虑细胞之间的空间关系，但一般只有在多个切片联合分析才考虑降维
      • 对比单细胞转录组分析的降维：基于基因表达矩阵进行降维，通常使用PCA、t-SNE、UMAP等方法
  3. 差异分析
    • 主要发现不同空间区域之间差异表达的基因，如SpatiaIDE
      • 对比单细胞转录组分析的差异分析：用于发现不同细胞类型之间的差异表达基因，如Seurat的FindAllMarkers
  4. 富集分析
    • 主要用于解释在不同区域或细胞类型中富集的生物学过程，如GO分析、KEGG分析
      • 对比单细胞转录组分析的富集分析：用于解释在不同细胞类型中富集的生物学过程，如GO分析、KEGG分析
  5. 细胞类型
    • 相互糅合在一起，强调统一体
      • 对比单细胞转录组分析的细胞类型：批次分开，强调细胞类型的精确
  6. 肿瘤细胞区域
    • 图片可以直接识别
      • 对比单细胞转录组分析的肿瘤细胞区域：算法推断
  7. 细胞类型分布
    • 能够在组织中确定不同细胞类型的空间分布情况
      • 对比单细胞转录组分析的细胞类型分布：一个或几个cluster
  8. 共定位
    • 通常指的是组织或者细胞区域内同时表达的基因，这可以通过分析不同区域的基因表达矩阵来确定，共定位基因可能在特定的细胞类型或空间区域中表达，从而解释这些区域的生物学功能或发育过程
      • 对比单细胞转录组分析的共定位：通常指的是不同细胞类型之间共同表达的基因，这可以通过分析单细胞数据集中不同细胞簇的基因表达模式来确定，共定位基因可能多个细胞类型中表达，从而揭示这些细胞类型之间的功能联系或分析交互作用
1
2
3
4
5
6
7
8
# 获取当前路径
import os

cur_path = os.getcwd()
print(cur_path)

os.chdir("/home/oelab/xxxx")
print(os.path)
1
2
!(python --version)
# !等价于调用shell命令

DataFrame

pandas 库中的 DataFrame 是一种二维表结构，可以理解为一个类似于 Excel 表格或数据库表的对象，能够存储不同类型的数据。它的特点是带有行标签（index）和列标签（columns），便于数据的组织、操作和分析

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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
import panads as pd
import os

# 读取类似的表格文件
df = pd.read_csv('data.csv')

# 显示最大的列
# display是显示的意思，None是全部列输出
pd.set_option('display.max_columns', None)

# 展示数据框的前8行
print(df.head(8))

# 展示数据框的后8行
print(df.tail(8))

# 查看数值列的统计信息
df.describe()

# 查看DataFrame的形状（行数、列数）
df.shape

# 查看数据框所有列的名字
df.columns

# 查看所有的行名
# 如果没有指定行名，默认从0开始，每次加1
df.index

# 修改指定的第几个行名
a = df.index
# a是把所有行名提取，它的数据类型是index，接下来转换为list列表
index_list = a.tolist()
index_list[3] = 'huge'
df.index = index_list
df.index

# 修改指定行名
df = df.rename(index = {'AAACAATGTATCTCCCCA': 'huge'})

# 修改列名
df = df.rename(columns={'in_tissue': 'new_in_tissue'}, inplace=True)

# 直接用中括号提取某一列
df['clusters']

# 提取多列
df[['clusters', 'array_row', 'array_col']]

# 提取一行
df.loc['AAACAATGTATCTCCCCA-1']

# 提取一列
df.loc[:, 'array_row']

# 提取指定第一行
df.iloc[0]

# 提取指定的第一列
df.iloc[:, 0]

# 提取第一行第一列
df.iloc[0, 0]

# 删除行、列
# axis = 1：删除列
# axis = 0：删除行
df.drop('huge', axis = 0)
df.drop('array_row', axis = 1)
df.drop(['array_row', 'array_col'], axis = 1)

# 按照条件进行过滤
df_new = df[df['clusters'] > 6]

# 数据框的合并
df1 = pd.DataFrame({
    'k1': ['A', 'B', 'C', 'D'],
    'k2': ['2', '3', '4', '5']
})

df2 = pd.DataFrame({
    'm1': ['A', 'B', 'C', 'D'],
    'm2': ['2', '3', '4', '5']
})

# 取交集合并两个数据框
merge1 = pd.merge(df1, df2, left_on = 'k1', right_on = 'm1', how = 'inner')
# how 存在 inner left right outer

# concat：指定轴（行或列）将多个数据框连接在一起，通常用于堆叠数据框，适用于数据框结构相同的情况
df_con1 = pd.concat([df1, df2], axis = 0)
df_con2 = pd.concat([df1, df2], axis = 1)

# 按cluster列分组
# 类似于R中的table函数
df_group1 = df.groupby('clusters')
# 查看每个分组的大小
df_group1.size()

# for循环
for i in [0, 2, 4, 6, 8]:
    m = i + 1
    n = m * 2
    print(n)
 

# 该缩进必须要缩进，必须整齐

# 函数
def fun_name(arg1, arg2 ...):
    code1
    code2

scverse分析框架 - AnnData数据结构

等同于R中的S4对象，是数据结构

目前单细胞领域的工具分为三种：

  1. Bioconductor：R语言实现的生物信息学生态
  2. Seurat：R语言实现的单细胞分析生态
  3. Scverse：基于Python实现的单细胞分析生态
    • 单细胞主要用anndata
    • 单细胞多组学主要有mudata
    • 空间的多组学用SpatialData

SCVI-Tools(Single Cell Variational Inference Tools)：是一个开源的Python库，专为单细胞转录组数据的分析而设计

squidpy：专门用于处理和分析空间分子数据，基于scanpy和anndata，继承了它们的模块化和可拓展性，并专门处理带有空间坐标的数据提供独特的方法，包括整合组织图像进行深度分析

AnnData结构：

AnnData数据结构

  1. 核心属性 X 为一个稀疏矩阵（scipy.sparse.csr.csr_matrix），用于存储矩阵的值
  2. 矩阵的行列分别对应观察者 obs（即细胞条形码），变量 var （即基因），这两个属性分别是一个数据框(panads.core.frame.DataFrame)类型用于存储行列的注释信息，其索引就是对应的 obs_names 和 var_names
功能 类型
adata.X 矩阵 numpy矩阵
adata.obs 观测量 pandas Dataframe
adata.var 特征量 pandas Dataframe
adata.uns 非结构化数据 字典dict

X数据槽：

  • X slot存储了当前用于分析的表达矩阵，一般使用CSR的稀疏矩阵格式
  • X 的一个数据结构是numpy array，是没有行名列名消息的，行名和列名消息存储在.obs和.var里
1
2
3
4
5
6
7
8
adata.X

# 稀疏矩阵转稠密矩阵的方法一
adata.X.A
# 稀疏矩阵转稠密矩阵的方法二
df1 = adata.X.todense()
# 稀疏矩阵转数据框
df2 = adata.to_df()

Layers数据槽：

  • Layers本质上是python的字典，每个layer对应一个表达矩阵的二位数组，比如标准化、对数化等等，这些转换后的数据可以存储在layer中
1
2
3
4
5
6
7
8
9
10
11
12
# 查看layers下哪些数据
adata.layers

# 增加layer
adata.layers["raw"] = adata.X

# 将raw counts数据复制到新的layer
adata.layers["count_per_million"] = adata.layers["raw"].copy()
adata.layers["log_transformed"] = np.log1p(adata.X)

# 还原回DataFrames，直接使用.to_df()
adata.to_df(layer="log_transformed")

使用不同的Layer数据画图：

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
# 创建一个包含两个子图的画布，并指定画布大小为(10, 2.5)，ax1和ax2分别代表两个子图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 2.5))

# 定义感兴趣的基因列表，这些基因将在矩阵图中进行展示
genes_of_interest = ["COL1A2", "SYPL1"]

sc.pl.matrixplot(
    adata, # 指定要绘制的数据对象
    groupby="clusters", # 指定按照细胞类型（由louvain算法聚类得到的结果）进行分组
    var_names=genes_of_interest, # 指定要在矩阵图中展示的基因
    layer="count_per_million", # 指定要绘制的数据层，这里选择了CPM(每百万计数)
    ax=ax1, # 指定当前绘制操作在每一个子图中进行
    show=False, # 设置为不限制绘制结果，因为将在后面手动调整布局
)

# 给第一个子图设置标题，说明这个子图中展示的是CPM归一化后的数据
ax1.set_title("CPM normalization")

sc.pl.matrixplot(
    adata, # 指定要绘制的数据对象
    groupby="clusters", # 指定按照细胞类型（由louvain算法聚类得到的结果）进行分组
    var_names=genes_of_interest, # 指定要在矩阵图中展示的基因
    layer="raw", # 指定要绘制的数据层，这里选择了CPM(每百万计数)
    ax=ax2, # 指定当前绘制操作在每一个子图中进行
    show=False, # 设置为不限制绘制结果，因为将在后面手动调整布局
)

# 给第一个子图设置标题，说明这个子图中展示的是CPM归一化后的数据
ax2.set_title("row counts")

# 调整子图布局，使得子图之间的间距合适
plt.tight_layout()

导出表达矩阵（就是10X单细胞那三个文件）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 导出layers中的矩阵
mat1 = adata.to_df(layer="log_transformed")
mat1.to_csv("exprMat_X.csv")

# 保存为MTX稀疏矩阵格式
# 创建文件夹
if not os.path.isdir("sparse_matrix"):
    os.makedirs("sparse_matrix")
    
# 导出稀疏矩阵
mmwrite("sparse_matrix/matrix.mtx", adata.X)

# 另存基因名和barcode列表，通常与mtx文件配套保存
with open("sparse_matrix/features.tsv", 'w') as f:
    f.write("gene_id\n")
    for gene in adata.var_names.values:
        f.write(f"{gene}\n")
with open("sparse_matrix/barcodes.tsv", 'w') as f:
    f.write("barcodes\n")
    for gene in adata.obs_names.values:
        f.write(f"{barcode}\n")

基础信息的储存：

obs就是metadata，var就是基因注释的数据框

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
# 查看细胞注释的内容
adata.obs

# 增加新的注释信息
adata.obs["is_low_quality"] = adata.obs["pct_counts_mt"] > 0.03

# 检查更新后的信息
adata.obs

# 查看细胞名字
adata.obs_names

# 查看基因
adata.var

# 查看基因名字
adata.var_names

# 查看多维注释信息
adata.obsm

# 整理成数据框输出
pd.DataFrame(adata.obsm['X_pca'], index = adata.obs_names, 
            columns=[f'PC_{i}' for i in range(1, 51)])

# 查看所有的细胞配对信息
adata.obsp

# 查看所有非结构化数据
adata.uns_keys()

使用不同降维坐标画图：

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
# 创建一个图形，并指定大小为（12， 4）
plt.figure(figsize=(12, 4))

# 创建一个包含1行3列的子图网络，并选择第一个子图作为当前操作对象
plt.subplot(1, 3, 1)
# PCA

df6 = adata.obs
plt.scatter( # 绘制散点图
	x = adata.obsm["X_pca"][:, 0], # 指定x轴数据为PCA降维后的第一列数据
	y = adata.obsm["X_pca"][:, 1], # 指定y轴数据为PCA降维后的第一列数据
	c = adata.obs["clusters"] == "3", # 根据是否属于B细胞类型，选择颜色，True表示B细胞
	s = 3, # 设置散点的大小为3
	linewidth = 0, # 设置散点边缘宽度为0，即不显示边缘
	cmap = "coolwarm", # 选择颜色映射
)

plt.title("PCA") # 设置当前子图的标题为“PCA”
plt.axis("off") # 关闭坐标轴
plt.gca().set_aspect("equal") # 设置坐标轴的纵横比为1，保证散点图中的点不会被拉伸

sc.tl.tsne(adata)

# t-SNE
plt.subplot(1, 3, 2)
plt.scatter( # 绘制散点图
	x = adata.obsm["X_tsne"][:, 0], 
	y = adata.obsm["X_tsne"][:, 1],
	c = adata.obs["clusters"] == "3",
	s = 3,
	linewidth = 0,
	cmap = "coolwarm",
)

plt.title("t-SNE")
plt.axis("off")
plt.gca().set_aspect("equal")

# UMAP
plt.subplot(1, 3, 2)
plt.scatter( # 绘制散点图
	x = adata.obsm["X_umap"][:, 0], 
	y = adata.obsm["X_umap"][:, 1],
	c = adata.obs["clusters"] == "3",
	s = 3,
	linewidth = 0,
	cmap = "coolwarm",
)

plt.title("UMAP")
plt.axis("off")
plt.gca().set_aspect("equal")

数据合并和筛选

  1. 两个对象的合并

    1
    2
    3
    import anndata as ad
    adatas = [adata1, adata2]
    adatas = ad.concatenate(adatas)

  2. metadata更改

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # 构建control组和处理组

    adata.obs['label'] = "NA"
    c = [i for i, word in enumerate(list(adata.obs['batch1'])) if word.startswith('A')]
    obs['label'][c] = "A"
    c = [i for i, word in enumerate(list(adata.obs['batch1'])) if word.startswith('B')]
    obs['label'][c] = "B"
    c = [i for i, word in enumerate(list(adata.obs['batch1'])) if word.startswith('C')]
    obs['label'][c] = "C"
    c = [i for i, word in enumerate(list(adata.obs['batch1'])) if word.startswith('Y')]
    obs['label'][c] = "Y"

  3. 数据保存和metadata添加

    1
    2
    3
    os.getcwd()
    adata.write("xxx.h5ad", compression='gzip')
    adata.obs.to_csv("scanpy_cell_info.csv", header=True, index=True)

  4. 载入数据

    1
    2
    3
    4
    5
    6
    7
    8
    9
    # 加入数据
    anno = pd.read_csv("scanpy_cell_info.csv", index_col = 1)
    # 加入样本分组的信息
    adata.obsp['cell_groups'] = anno['label']

    # 甚至可以直接赋值 dataframe
    adata.obs['cell_groups1'] = list(anno['clusters'])
    adata.obs = anno
    df = adata.obs

  5. 数据的提取

    1
    2
    3
    # : 表示基因的过滤
    adata = adata[adata.obs.n_genes_by_counts < 2500, :]
    adata = adata[adata.obs.pct_counts_mt < 5, :]

空间转录组基本分析流程

流程：

  • 导包

    1
    2
    3
    4
    5
    6
    7
    8
    9
    import os
    import anndata
    import scanpy as sc
    import matplotlib.pyplot as plt
    import numpy as np
    import seaborn as sns
    import pooch
    from scipy.sparse import csr_matrix
    from scipy.io import mmwrite

  • 查看版本号等

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # 输出版本号
    sc.logging.print_versions()
    # 将图片的背景改为白色
    # 设置输出图像格式
    sc.set_figure_params(facecolor="white", figsize=(8, 8))
    # verbosity 的取值表示测试结果显示的详细程度，数字越大越详细
    # errors = 0
    # warnings = 1
    # info = 2
    # hints = 3
    sc.settings.verbosity = 3

  • 更改当前的路径

    1
    2
    3
    4
    5
    6
    7
    8
    # 获取当前工作路径
    print(os.getcwd())

    # 更改路径
    os.chdir('/home/oelab/RunData2/STM_AD_analyze/data/GSE220442_counts_and_images/')

    ## 更改完再检查一下有没有更改成功
    print(os.getcwd())

  • 读取数据

    1
    adata = sc.read_visium('../GSE220442_counts_and_images/counts_and_images/AD1/')

  • 查看AnnData的数据结构

    1
    adata

  • 标记线粒体基因

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    # 线粒体百分比计算
    # 进行质控

    # 将基因名去重
    adata.var_names_make_unique()

    # 标记线粒体基因
    # 在adata的数据里增加了一行mt
    """
    用于标记 adata 中的线粒体基因。这里假设线粒体基因的名称以 "MT-" 开头，adata.var["mt"] 会存储一个布尔值数组，标记每个基因是否为线粒体基因
    """
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    # pp指对空转数据的预处理
    """
    sc.pp.calculate_qc_metrics 是 Scanpy 中的一个函数，用于计算质量控制指标。
    参数 qc_vars=["mt"] 告诉函数将 adata.var["mt"] 中标记的线粒体基因作为一个单独的质量控制变量。
    这样可以计算每个细胞中线粒体基因的表达比例等指标，结果会直接添加到 adata.obs 中
    """
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

  • 两种赋值方法

    1
    2
    3
    4
    5
    # 第一种赋值方式，浅拷贝
    adata1=adata

    # 第二种赋值方式，深拷贝
    adata1=adata.copy()

  • 查看质控图

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    ## 使用 plt.subplots 创建一个包含 1 行、4 列的图形布局，每个子图的大小为 (15, 4)
    fig, axs = plt.subplots(1, 4, figsize=(15, 4))
    ## 第一个子图中绘制 adata.obs["total_counts"] 的直方图，用于显示每个细胞的总计数分布情况
    sns.histplot(adata.obs["total_counts"], kde=False, ax=axs[0])
    ## 第二个子图中绘制一个限定了范围的总计数直方图，仅包括总计数少于 10000 的细胞，以便更清晰地观察低计数细胞的分布情况
    sns.histplot(adata.obs["total_counts"][adata.obs["total_counts"] < 10000],
        kde=False, bins=40,ax=axs[1],)
    ## 第三个子图中绘制 adata.obs["n_genes_by_counts"] 的直方图，用于显示每个细胞中检测到的基因数分布
    sns.histplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])
    ## 第四个子图中绘制基因计数少于 4000 的直方图，方便观察低基因计数的细胞分布
    sns.histplot(adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000],
        kde=False, bins=60, ax=axs[3],)

  • 过滤低质量的细胞

    1
    2
    3
    4
    5
    sc.pp.filter_cells(adata, min_counts=5000)
    sc.pp.filter_cells(adata, max_counts=35000)
    adata = adata[adata.obs["pct_counts_mt"] < 20].copy()
    print(f"#cells after MT filter: {adata.n_obs}")
    sc.pp.filter_genes(adata, min_cells=10)

  • 数据归一化

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

  • pca降维，umap可视化，leiden的聚类

    1
    2
    3
    4
    5
    sc.pp.pca(adata)
    sc.pp.neighbors(adata)
    sc.tl.umap(adata)
    # 类似于R的findclusters
    sc.tl.leiden(adata, key_added="clusters" )

  • 展示umap图上基因的表达量

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    # 过滤之后的数据是什么样的
    plt.rcParams["figure.figsize"] = (4, 4)
    # 将三张图放在一张图里面
    sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

    # 和上面的结果是一样的
    # 分别展示三张图的结果
    sc.pl.umap(adata, color=["total_counts" ], wspace=0.4)
    sc.pl.umap(adata, color=[ "n_genes_by_counts"], wspace=0.4)
    sc.pl.umap(adata, color=[ "clusters"], wspace=0.4)

  • 可视化空间信息

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    # 查看空转位置上totalcounts和基因数量的分布
    plt.rcParams["figure.figsize"] = (8, 8)
    sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])

    # 查看每个clusters在空间上是如何分布的
    sc.pl.spatial(adata, img_key="hires", color="clusters", size=1.5)

    # 筛选出1和3两个cluster的分布
    sc.pl.spatial(adata,
        img_key="hires",color="clusters",
        groups=["1", "3"],crop_coord=[7000, 10000, 0, 6000],
        alpha=0.5,size=1.3,)

  • 差异分析

    1
    2
    3
    4
    # 找出来marker基因
    sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
    # 基于group2进行展示
    sc.pl.rank_genes_groups_heatmap(adata, groups="2", n_genes=10, groupby="clusters")

  • 个性化展示画图

    1
    2
    3
    4
    5
    # 查看cluster的信息
    sc.pl.spatial(adata, img_key="hires", color=["clusters" ])

    # 查看两个基因表达量的分布
    sc.pl.spatial(adata, img_key="hires", color=["COL1A2", "SYPL1"], alpha=0.7)

  • 对数据进行保存，方便下一次的读取

    1
    2
    os.chdir("/home/oelab/RunData2/STM_AD_analyze/")
    adata.write("GSE220442.data.anndata.h5ad")

  • 保存画图

    1
    fig.savefig(f'./figures/Visium-HE_{reg}.pdf')

scanpy 单个样本的单细胞分析流程

基本的分析流程：

  1. 对初始的Adata的预处理
  2. 主成分分析
  3. 计算 neighborhood graph
  4. 对 neighborhood graph 进行 embedding
  5. 对 neighborhood graph 进行聚类
  6. 找到 marker 基因（簇间比较）
  7. 根据 marker 基因注释

流程：

  1. scanpy单细胞数据读取

    • 公共数据常用的单细胞数据文件格式：
      • Matrix Market Exchange Format (MEX)：
        • .mtx 文件，该文件包含了一个稀疏矩阵，其中每一行代表一个基因，每一列代表一个细胞，而每个素则表示基因在相应细胞中的表达量。通常会有三个文件组成：barcodes.tsv（细胞标识）features.tsv（基因标识）和 matrix.mtx（表达矩阵）
      • Hierarchical Data Format (HDF5)：
        • h5 文件，HDF5 是一种用于存储和组织大规模科学数据的文件格式。在单细胞数据中，可以使用HDF5 存储表达矩阵、基因和细胞的元数据等信息
      • AnnData：
        • .h5ad 文件，是基于 HDF5 的文件格式，由 anndata 库定义。它可以包含单细胞数据的表达矩阵、基因和细胞的元数据、样品信息等
      • Comma-Separated Values（CSV）或文本文件：
        • .csv 或 .txt 文件，包含表达矩阵的文本文件。每一行代表一个基因，每一列代表一个细胞，元素是表达量

  2. 读取数据（以mtx为例，但是10X单细胞三个文件得解压）

    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
    # 导入scanpy库
    import scanpy as sc

    # 读取tab文件
    data = sc.read_10X_mtx("/GSE121121")

    # 查看数据
    data

    ----------------------------------------------------------------------------------------
    # 读取h5文件
    # 导入scanpy库
    import scanpy as sc

    # 读取h5文件
    data = sc.read_10X_h5("/xxx.h5")

    # 查看数据
    data

    # 读取h5ad文件
    adata = sc.read_h5ad('/xxx.h5ad')

    # 读取csv文件
    data4 = sc.read_csv("/xxx.csv")

    # 读取txt文件
    data5 = sc.read_text("/xxx.txt")

  3. 加载模块

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    import os
    import anndata
    import scanpy as sc
    import pandas as pd
    import matplotlib.pyplot as plt
    import numpy as np
    import seaborn as sns
    import pooch
    from scipy.sparse import csr_matrix
    from scipy.io import mmwrite

  4. 初始化设置

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

  5. 读取mtx文件

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    print(os.getcwd())
    os.chdir("/home/oelab/RunData2/STM_AD_analyze/data/")
    print(os.getcwd())

    # 声明h5ad用于存储分析结果
    result_file = 'sample1_h5ad'

    adata = sc.read_10x_mtx("../data/AD1/", # mtx文件所在的目录
                           var_names='gene_symbols', # 将 gene 作为 var
                           cache=True)  # 开启缓存读写

    # 消除重复的基因
    adata.var_names_make_unique()

  6. 对初始adata的预处理

    1
    2
    # 绘制了平均表达量前10的基因箱型图
    sc.pl.highest_expr_genes(adata, n_top=10)

  7. 质量控制

    1
    2
    3
    4
    5
    # 基本的过滤(质控)

    # pp指 pre process 预处理
    sc.pp.filter_cells(adata, min_genes=200)
    sc.pp.filter_genes(adata, min_cells=3)

  8. 线粒体百分比的计算

    1
    2
    3
    4
    5
    # 线粒体百分比的计算
    adata.var['mt'] = adata.var_names.str.startswith('MT-')
    adata.var['mt']
    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
    adata.obs

  9. 小提琴图展示

    1
    2
    3
    4
    5
    # total_counts：每个细胞的基因总计数（总表达量）
    # pct_counts_mt：每个细胞中，线粒体基因表达量占该细胞所有基因表达量的百分比

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

  10. 进行过滤

    1
    2
    3
    # 进行过滤
    adata = adata[adata.obs.n_genes_by_counts > 500, :]
    adata = adata[adata.obs.pct_counts_mt < 20, :]

  11. 进行归一化处理

    1
    2
    3
    # 进行归一化处理
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)

  12. 高变基因的展示

    1
    2
    3
    4
    5
    # 高变基因的展示
    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
    sc.pl.highly_variable_genes(adata)

    adata.raw = adata

  13. 提取高变基因的表达矩阵

    1
    2
    3
    4
    5
    6
    7
    8
    # 提取高变基因的表达矩阵
    adata = adata[:, adata.var.highly_variable]

    # 对两个注释进行回归，消除它俩的影响
    sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])

    # 按零均值单位方差标准化数据，并裁剪超过标准差10的细胞
    sc.pp.scale(adata, max_value=10)

  14. 主成分分析（PCA）

    1
    2
    3
    4
    5
    6
    sc.tl.pca(adata, svd_solver='arpack')

    sc.pl.pca_variance_ratio(adata, log=True)

    # 将结果保存到声明的文件results_file
    adata.write(result_file)

  15. 计算 neighborhood graph

    1
    2
    # 使用PCA来计算细胞的领域图，为了umap降维，也为了聚类
    sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

  16. 对 neighborhood graph 进行 embedding

    1
    2
    3
    4
    5
    6
    7
    8
    # 进行umap的可视化
    sc.tl.umap(adata)

    # 进行umap的展示
    sc.pl.umap(adata, color=['NKG7'])

    # 可以不适用原始数据
    sc.pl.umap(adata, color=['NKG7'], use_raw=False)

  17. 对 neighborhood graph 进行聚类

    1
    2
    3
    4
    5
    # scanpy推荐leiden图聚类的方法
    sc.tl.leiden(adata)

    # 用umap对leiden聚类结果的可视化
    sc.pl.umap(adata, color=['leiden', 'NKG7'])

  18. 找到marker基因

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    # 计算每个簇存在较高程度差异的基因排名
    # 在计算差异时，最简单最快的方式是t-检验
    sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')

    # 展示差异基因
    # 每个cluster的top20或者top50marker基因的展示
    sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

    # 还可以使用wilcoxon的检测方式
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
    sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

    # 差异基因以表格形式查看
    pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)

    # 显示差异基因的名称和差异度
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
    sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
    result = adata.uns['rank_genes_groups']
    groups = result['names'].dtype.names
    pd.DataFrame(
        {group + '_' + key[:1] : result[key][group]
        for group in groups for key in ['names', 'pvals']}
    ).head(5)

  19. marker基因展示

    1
    2
    3
    4
    5
    6
    7
    8
    # 用小提琴图展示前20个差异基因在该簇和剩余簇的分布差异
    from matplotlib.pyplot import rc_context

    with rc_context({'figure.figsize': (9, 1.5)}):
        sc.pl.rank_genes_groups_violin(adata, n_genes=20, jitter=False)
        
    # 直观地显示某个基因在不同簇中的表达
    sc.pl.violin(adata, ['NKG7'], groupby='leiden')

  20. 细胞类型注释

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    # 每个cluster对应的细胞类型注释
    new_cluster_names = [
        '0', '1', '2':'CD14', '3':'endo', '4', '5':'NK'
    ]

    adata.rename_categories('leiden', new_cluster_names)

    # 可视化结果展示
    pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False)

    # 用tracksplot验证基因与细胞的关系
    sc.pl.rank_genes_groups_tracksplot(adata, n_genes=3)

多个样本单细胞数据整合 - ingest和BBKNN

整合算法介绍

ingest

  • 一种基于PCA的简单方法用于集成数据（integrating data）
  • 用已知的细胞类型数据去注释未知的
  • ingest函数假设有一个带注释的参考数据集（有点类似SingleR），该数据集捕获了感兴趣的生物变异性。ingest是在参考数据上拟合模型，并使用它来投影新数据。目前，该模型是一个结合了邻居查找搜索树（neighbor lookup search tree）的PCA或者UMAP

ingest有以下特点：

  1. 由于ingest操作简单，过程清晰，因此工作流透明且快速
  2. 与BBKNN相似，ingest使数据矩阵本身保持不变。因为ingest是在embedding空间做分析
  3. 与BBKNN不同，ingest在解决了标签映射的同时，维持了可能具有特定簇或轨迹等所需属性的embedding（因为BBKNN是简单粗暴的对齐embedding）

BBKNN

  • BBKNN：Batch balanced KNN，BBKNN是一种快速直观的去除批次效应工具，可直接用于scanpy工作流。它是 scanpy.pp.neighbors() 的一种替代方法，这两个功能都是创建一个 neighbors graph，以便后续在聚类、pseudo time计算和UMAP可视化中使用
  • 标准方法首先确定整个数据集中每个细胞的k个最近邻，然后将候选细胞转换为指数相关连接，然后作为进一步分析的基础。如果数据中存在批次效应（可能是由于不同的数据采集技术、协议变更，甚至是特别严重的操作员影响），那么将很难在不同批次中链接相应的细胞类型
  • BBKNN修改了邻域构建步骤，以生成一个跨所有批次数据平衡的图。这种方法将邻居网络视为数据的主要表示形式。对于每个细胞，BBKNN图是通过独立地为每个设备提供的批次（单个数据集）中的每个细胞找到k个最近邻居来构建的，这使得每个细胞在每个批次中都有一个独立的邻居池。同类细胞的这些邻居都被加入列表联系起来
  • 注意：BBKNN会将列表中的元素用简单快速的方式对齐，即对齐不同批次的同类细胞。所以才说BBKNN实现了数据整合
  • BBKNN的主要假设是，批次间至少存在一些相同类型的细胞，并且批次效应引起的批次间相同细胞类型之间的差异小于批次内不同类型细胞之间的差异，在这种情况下，图的构造将跨批次将相似的细胞类型分组在一起，同时将不相关的细胞类型分开

使用BBKNN进行数据整合

1
2
3
4
5
6
7
8
9
10
11
from bbknn import bbknn

# BBKNN整合
sc.external.pp.bbknn(adatas, batch_key='batch')

sc.tl.umap(adatas)
sc.pl.umap(adatas, color=['batch', 'orig'])

adatas.obsm

sc.tl.leiden(adatas)

基于scanorama整合scRNA-seq和空间数据分析

官方连接：Spatial — scanpy

与单细胞数据具体整合流程见：使用 scanorama 将空间数据与 scRNA-seq 集成 — scanpy

  • Scanorama可以实现scRNA-seq数据集集成和批次矫正
  • 文章提出的方法是整合scRNA-seq的高效方案，与过去的ingest不同，不需要指定某个具体的参考数据集，其实过去方法的思想还是将数据整合到一个更大的簇中（域适应到参考数据集）
  • Scanorama利用全景生成的思想，可以实现同时多个数据集的整合，可以自动匹配不同数据集下的相同细胞类型，同时保留不同数据集下不同细胞类型的差异，是更合理的整合方法
  • 其中的异质指的是整合保留了两个数据集中真正存在差异的信息。Scanorama属于是实现了正确消除了技术，实验带来的批次效应

流程：

  1. 导包

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    import os
    import anndata
    import scanpy as sc
    import pandas as pd
    import matplotlib.pyplot as plt
    import numpy as np
    import seaborn as sns
    import pooch
    from scipy.sparse import csr_matrix
    from scipy.io import mmwrite
    import scanorama

  2. 设置一些全局参数

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

  3. 设置路径

    1
    2
    3
    print(os.getcwd())
    os.chdir("/home/oelab/RunData2/STM_AD_analyze/data/GSE220442_counts_and_images/counts_and_images/")
    print(os.getcwd())

  4. 读取数据

    1
    2
    3
    4
    5
    6
    7
    # 读取数据
    AD1 = sc.read_visium("../counts_and_images/AD1/")
    AD2 = sc.read_visium("../counts_and_images/AD2/")
    AD3 = sc.read_visium("../counts_and_images/AD3/")
    NC1 = sc.read_visium("../counts_and_images/NC1/")
    NC2 = sc.read_visium("../counts_and_images/NC2/")
    NC3 = sc.read_visium("../counts_and_images/NC3/")

  5. 质控计算

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    # 质控计算
    AD1.var_names_make_unique()
    AD2.var_names_make_unique()
    AD3.var_names_make_unique()
    NC1.var_names_make_unique()
    NC2.var_names_make_unique()
    NC3.var_names_make_unique()

    sc.pp.calculate_qc_metrics(AD1, inplace=True)
    sc.pp.calculate_qc_metrics(AD2, inplace=True)
    sc.pp.calculate_qc_metrics(AD3, inplace=True)
    sc.pp.calculate_qc_metrics(NC1, inplace=True)
    sc.pp.calculate_qc_metrics(NC2, inplace=True)
    sc.pp.calculate_qc_metrics(NC3, inplace=True)

  6. 可视化

    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
    # 可视化
    for name,adata in [
        ("AD_sample1", AD1),
        ("AD_sample2", AD2),
        ("AD_sample3", AD3),
        ("Normal_sample1", NC1),
        ("Normal_sample2", NC2),
        ("Normal_sample3", NC3)]:
        fig, axs = plt.subplots(1, 4, figsize= (12, 3))
        fig.suptitle(f"Covariates for filtering: {name}")

        sns.distplot(adata.obs["total_counts"], kde=False, ax = axs[0])
        
        sns.distplot(
            adata.obs["total_counts"][adata.obs["total_counts"] < 20000],
            kde=False, bins=40, ax=axs[1],)
        
        sns.distplot(
            adata.obs["n_genes_by_counts"],
            kde=False, bins=60, ax=axs[2],)

        sns.distplot(
            adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000],
            kde=False, bins=60, ax=axs[3],)
        plt.savefig("/home/oelab/RunData2/STM_AD_analyze/result/" + "qc1.distplot_" + name + ".png")

  7. 数据整合

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    AD1.obs["orig"] = 'AD'
    AD2.obs["orig"] = 'AD'
    AD3.obs["orig"] = 'AD'
    NC1.obs["orig"] = 'Normal'
    NC2.obs["orig"] = 'Normal'
    NC3.obs["orig"] = 'Normal'

    adatas = [AD1, AD2, AD3, NC1, NC2, NC3]

    # 数据整合
    adatas_cor = scanorama.correct_scanpy(adatas, return_dimred=True)

  8. 将多个数据整合到一起并保存在整合后的embeddings中

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    adatas_spatial = sc.concat(
        adatas_cor, label = "library_id",
        uns_merge="unique", 
        keys=[
            k
            for d in [
                adatas_cor[0].uns["spatial"],
                adatas_cor[1].uns["spatial"],
                adatas_cor[2].uns["spatial"],
                adatas_cor[3].uns["spatial"],
                adatas_cor[4].uns["spatial"],
                adatas_cor[5].uns["spatial"]
            ]
            for k, v in d.items()
        ],
        index_unique="-"
    )

  9. 降维聚类

    1
    2
    3
    4
    5
    6
    7
    8
    sc.pp.neighbors(adatas_spatial, use_rep="X_scanorama")
    sc.tl.umap(adatas_spatial)

    sc.tl.leiden(adatas_spatial, key_added="clusters")

    sc.pl.umap(
        adatas_spatial, color=["clusters", "library_id", "orig"], palette=sc.pl.palettes.default_20
    )

  10. 空间可视化

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    clusters_colors = dict(
        zip([str(i) for i in range(18)], adatas_spatial.uns["clusters_colors"])
    )

    # 开始画图
    fig, axs = plt.subplots(1, 6, figsize= (15, 10))

    for i, library in enumerate(
        ["2-8", "T4857", "2-3", "2-5", "18-64", "1-1"]
    ):
        ad = adatas_spatial[adatas_spatial.obs.library_id == library,:].copy()
        sc.pl.spatial(
            ad, img_key="hires", library_id=library,
            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()

Squidpy

教程：Squidpy - Spatial Single Cell Analysis in Python — squidpy documentation

ImageContainer对象

展示如何使用 squidpy.im.ImageContainer 与图像结构化数据交互

教程：ImageContainer object — squidpy documentation

流程：

  1. 初始化 ImageContainer

    1
    2
    3
    4
    5
    6
    7
    import squidpy as sq
    import scanpy as sc
    import anndata as ad
    import numpy as np
    import pandas as pd
    import os
    import matplotlib.pyplot as plt
    1
    2
    # 初始化ImageContainer -- 就是创建ImageContainer数据结构的对象
    img = sq.im.ImageContainer(AD1.uns["spatial"]["2-8"]["images"]["hires"], library_id = "2-8", layer="image")

  2. 添加 layers 到 ImageContainer

    1
    2
    3
    4
    5
    6
    7
    # 添加layers到ImageContainer
    ## uns中存储的是空间信息
    img.add_img(AD1.uns["spatial"]["2-8"]["images"]["hires"], 
                library_id = "2-8", layer="hires")

    # 指定属性值，优先展示哪个属性
    img["hires"].attrs["hires"] = True

  3. 可视化

    1
    2
    3
    4
    5
    6
    7
    # 可视化
    img.show(layer = "image")

    img.show("image", segmentation_layer="segmentation", segmentation_alpha=0.5)

    for i, cmap in zip(range(3), ['Reds', 'Greens', 'Blues']):
        img.show(layer='image', channel=i, cmap=cmap)

  4. 裁剪和缩放图像

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    # 图像的裁剪和放缩
    # x=0.5 代表在横坐标中间的位置
    # y=0.5 代表在纵坐标中间的位置

    # crop_center：围绕中心进行裁剪
    # 左图
    img.crop_center(y = 0.5, x = 0.5, radius=200).show(layer="image")
    # 右图
    img.crop_center(y = 0.5, x = 0.5, radius=100).show(layer="image")

    # crop_corner：是将这个点作为起始点截取向右、向上的正方形

    img.crop_corner(y = 0.2, x = 0.3, radius=400).show(layer="image")
    img.crop_corner(y = 0.2, x = 0.3, radius=200).show(layer="image")


    # 查看裁剪之后的图像
    sc.pl.spatial(
        adata_crop[adata_crop.obs["library_id"] == "2-8"], color = ["CD4"], 
        library_id="2-8", frameon=False, cmap="Spectral_r", size=1.3
    )

  5. 处理图像和提取特征

  6. 图像的内部表示

添加多个样本，对多个样本进行提取：

1
2
3
4
5
6
7
8
# 添加多个样本，对多个样本进行提取
img = []
for library_id in adata.uns["spatial"].keys():
    img = sq.im.ImageContainer(adata.uns["spatial"]["library_id"]["images"]["segmentation"])
    img.add_img(adata.uns["spatial"]["library_id"]["images"]["segmentation"])
    img["segmentation"].attrs["segmentation"] = True
    imgs.append(img)
img = sq.im.ImageContainer.concat(imgs)

使用Napari进行交互式可视化

就是展示如何在Napari中可视化ImageContainer和AnnData

使用 Napari 可视化图像并与之交互 — 使用 Napari 和 Python 进行图像分析（亥姆霍兹成像暑期学院 2022 年）

流程：

  1. 安装Napari

    1
    2
    !(conda install -c conda-forge napari pyqt -y)
    pip install deprecated

  2. 加载数据

    1
    2
    3
    4
    5
    6
    import squidpy as sq
    import os
    import scanpy as sc
    adata = sq.datasets.visium_hne_adata()
    adata = sc.read("xxxx")
    img = sq.datasets.visium_hne_image()

  3. 启动

    1
    2
    # 运行viewer
    viewer = img.interactive(adata)

Seurat 空转数据分析

链接：Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurat

代码注释见上面的链接：

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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
# 加载R包
# 部署当前环境

# 1.系统报错改为英文
Sys.setenv(LANGUAGE = "en")

# 2.禁止转化为因子
options(stringsAsFactors = FALSE)

# 3.清空环境
rm(list=ls())

library(SeuratData)
library(ggplot2)
library(patchwork)
library(dbplyr)
library(DropletUtils)
library(tibble)
library(Seurat)
library(jsonlite)
library(stringr)

setwd("/home/oelab/RunData2/STM_AD_analyze/")

# 读取10X visium数据
data <- Load10X_Spatial("../STM_AD_analyze/data/R_GSE220442/AD1/")


# 处理数据
plot1 <- VlnPlot(data, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot1

plot2 <- SpatialFeaturePlot(data, features = "nCount_Spatial") + theme(legend.position = "right")

wrap_plots(plot1, plot2)

# 进行归一化
data <- SCTransform(data, assay = "Spatial", verbose = FALSE)

# 基因可视化
SpatialFeaturePlot(data, features = c("IL1A", "IL1B", "IL6", "IL11", "IL10",  
                                     "IL18", "IL33", "IL12A", "IL12B", "IL15", 
                                     "IL23A", "IL1RN", "IL32"))

# pt.size.factor 更改点的大小
SpatialFeaturePlot(data, features = c("IL1A", "IL1B", "IL6", "IL11", "IL10",  
                                      "IL18", "IL33", "IL12A", "IL12B", "IL15", 
                                      "IL23A", "IL1RN", "IL32"), pt.size.factor = 1)

# alpha = c(0.1, 1) 就是去降低低表达点的透明度
SpatialFeaturePlot(data, features = c("IL1A", "IL1B", "IL6", "IL11", "IL10",  
                                      "IL18", "IL33", "IL12A", "IL12B", "IL15", 
                                      "IL23A", "IL1RN", "IL32"), alpha = c(0.1, 1))


# 降维、聚类、可视化
data <- RunPCA(data, assay = "SCT", verbose = FALSE)
data <- FindNeighbors(data, reduction = "pca", dims = 1:30)
data <- FindClusters(data, verbose = FALSE)
data <- RunUMAP(data, reduction = "pca", dims = 1:30)

p1 <- DimPlot(data, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(data, label = TRUE, label.size = 3)
p1 + p2

# 对于指定的cluster进行展示
SpatialDimPlot(data, cells.highlight = CellsByIdentities(object = data, idents = c(2, 1, 4)), facet.highlight = TRUE, ncol = 3)

# 进行交互式的画图
SpatialDimPlot(data, interactive = TRUE)

# 基因交互式展示，可以展示基因的大小亮度等
SpatialFeaturePlot(data, features = "IL18", interactive = TRUE)

# 识别空间的高变基因
de_markers <- FindMarkers(data, ident.1 = 5, ident.2 = 6)
SpatialFeaturePlot(object = data, features = rownames(de_markers)[1:3], alpha = c(0.1, 1), ncol = 3)

# 子集化切片区域
cortex <- subset(data, idents = c(1, 2, 3, 4, 6, 7))

p1 <- SpatialDimPlot(cortex, crop = TRUE, label = TRUE)
p2 <- SpatialDimPlot(cortex, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3)
p1 + p2


# 空转与单细胞数据整合
# 加载单细胞的数据进行预测作为新的分析
singleCelldata <- readRDS("/home/oelab/RunData2/STM_AD_analyze/scRNA_SCTrandform.rds")
# singleCelldata <- readRDS("/home/oelab/RunData2/AD_Analyze/scRNA_harmony.rds")


# 设置ncells=3000使整个数据集归一化，在3k上学习噪声模型，加快了SCTransform的速度，而性能没有损失
# Step 1: SCTransform
options(future.globals.maxSize = 10000 * 1024^2)  # 10GB
singleCelldata <- SCTransform(
  singleCelldata,
  verbose = FALSE,
  return.only.var.genes = FALSE
)

# Step 2: RunPCA
singleCelldata <- RunPCA(singleCelldata, verbose = FALSE)

# Step 3: RunUMAP
singleCelldata <- RunUMAP(singleCelldata, dims = 1:30)

# saveRDS(singleCelldata, "/home/oelab/RunData2/STM_AD_analyze/scRNA_SCTrandform.rds")

# 重新运行SCTransform，确保只有一个模型
singleCelldata <- SCTransform(singleCelldata, verbose = FALSE, return.only.var.genes = FALSE)


# 对皮质数据集进行提取

# 用单细胞来注释空转
# 查看所有模型的名称
names(singleCelldata[["SCT"]]@SCTModel.list)

# 假设你只需要第一个SCT模型，可以删除其他模型
singleCelldata@assays$SCT <- singleCelldata@assays$SCT@SCTModel.list[[1]]


# After subsetting, we renormalize cortex
data <- SCTransform(data, assay = "Spatial", verbose = FALSE) %>%
  RunPCA(verbose = FALSE)

# the annotation is stored in the 'subclass' column of object metadata
DimPlot(singleCelldata, group.by = "celltype", label = TRUE)

# 设置更大的内存限制（例如，设置为 10GB）
options(future.globals.maxSize = 10 * 1024^3)  # 10GB

# singleCelldata <- IntegrateLayers(object = singleCelldata, method = CCAIntegration, normalization.method = "SCT", verbose = F)

# 假设选取的SCT的第一个模型
singleCelldata[["SCT"]]@SCTModel.list <- list(singleCelldata[["SCT"]]@SCTModel.list[[1]])


anchors <- FindTransferAnchors(reference = singleCelldata, query = data, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = singleCelldata@meta.data[["celltype"]], prediction.assay = TRUE, weight.reduction = data[["pca"]], , dims = 1:30)
data[["predictions"]] <- predictions.assay

data[["predictions"]] <- predictions.assay

DefaultAssay(data) <- "predictions"
SpatialFeaturePlot(data, features = c("Endo.", "Olig."), pt.size.factor = 1.6, ncol = 2, crop = TRUE)

data <- FindSpatiallyVariableFeatures(data, assay = "predictions", selection.method = "moransi",
                                        features = rownames(data), r.metric = 5, slot = "data")

result <- SpatiallyVariableFeatures(data, method = "moransi")
str(result)  # 查看对象结构

top.clusters <- head(SpatiallyVariableFeatures(data, method = "moransi"), 4)

SpatialFeaturePlot(data, features = c("Astrocyte", "OPC", "Olig.", "Neuron", "Micro.", "Endo."), pt.size.factor = 1, ncol = 2, crop = FALSE, alpha = c(0.1, 1))


# 多个空转样本的整合
AD1 <- Load10X_Spatial("../STM_AD_analyze/data/R_GSE220442/AD1/")
AD1 <- SCTransform(AD1, assay = "Spatial", verbose = FALSE)
AD2 <- Load10X_Spatial("../STM_AD_analyze/data/R_GSE220442/AD2/")
AD2 <- SCTransform(AD2, assay = "Spatial", verbose = FALSE)
AD3 <- Load10X_Spatial("../STM_AD_analyze/data/R_GSE220442/AD3/")
AD3 <- SCTransform(AD3, assay = "Spatial", verbose = FALSE)

# 为每个对象添加独特的前缀
AD.merge <- merge(
  AD1, 
  y = list(AD2, AD3), 
  add.cell.ids = c("AD1", "AD2", "AD3")
)

DefaultAssay(AD.merge) <- "SCT"
VariableFeatures(AD.merge) <- c(VariableFeatures(AD1), VariableFeatures(AD2), VariableFeatures(AD3))
AD.merge <- RunPCA(AD.merge, verbose = FALSE)
AD.merge <- FindNeighbors(AD.merge, dims = 1:30)
AD.merge <- FindClusters(AD.merge, verbose = FALSE)
AD.merge <- RunUMAP(AD.merge, dims = 1:30)

DimPlot(AD.merge, reduction = "umap", group.by = c("ident", "orig.ident"))

SpatialDimPlot(AD.merge)

SpatialFeaturePlot(AD.merge, features = c("CD34"))

做图技巧

去除黑圈：

1
2
3
4
5
6
7
# 去掉黑圈

p1 <- SpatialDimPlot(AD1)
p2 <- SpatialDimPlot(AD1, stroke = NA)
p2 <- p2 + guides(color=guide_legend(override.aes = list(size = 8), ncol = 2))
p2 <- p2 + theme_gray() + xlab("") + ylab("")
p2 <- p2 + theme(axis.text = element_blank(), axis.ticks = element_blank())

利用seurat单细胞的函数实现空转spot分布：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 利用seurat单细胞的函数实现空转spot分布

spatial_corr <- AD1@images$slice1@boundaries$centroids[, c("col", "row")]
colnames(spatial_corr) <- c("s_1", "s_2")
spatial_corr <- as.matrix(spatial_corr)
AD[["spatial"]] <- CreateDimReducObject(embeddings = spatial_corr, key = )


p1 <- FeaturePlot(AD1, features = "nCount_Spatial", reduction = "spatial") +
  scale_y_reverse() +
  scale_colour_viridis(option = "inferno") + theme_void()

p1 <- FeaturePlot(AD1, features = "Hpca", reduction = "spatial") +
  scale_y_reverse() +
  scale_colour_viridis(option = "D") + theme_void()

mycolor <- c("#123456", "#121214")

p3 <- DimPlot(AD1, reduction = "spatial", cols = mycolor) +
  guides(color=guide_legend(override.aes = list(size=6), ncol = 2)) +
  theme_void() + scale_y_reverse()

p1 + p2 + p3 + plot_layout(ncol = 3, nrow = 1)

Giotto

Giotto工具箱是Giotto软件包的升级版本，提供了处理、分析和可视化各种尺度和多种分辨率的空间多组学数据的工具。它的基本框架可推广到几乎所有当前和新兴的空间组学技术上，例如：原位杂交（seqFISH+, merFISH, osmFISH, CosMx）、测序（Slide-seq, Visium, STARmap, SeqScope, Stereo-Seq）以及基于成像的多组学技术（CyCIF, MIBI, CODEX）等。这些技术在分辨率（亚细胞、单细胞或多细胞）、空间维度（2D vs 3D）、分子形态（蛋白质、RNA、DNA等）和通量（细胞和分析物的数量）方面有所不同

Seurat调用不了Python，但是Giotto可以调用Python

链接：GitHub - KPMP/Cell-State-Atlas-2022: Code for a multimodal atlas of healthy and injured cell states and niches in the human kidney

Giotto链接：Spatial Single-Cell Transcriptomics Toolbox - Giotto

配置环境：

1
conda create -n giotto_env python=3.10.2 -c https://mirror.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
1
conda activate giotto_env

各个包的版本：Giotto python 环境 — giotto_python • GiottoClass

1
2
3
4
5
6
7
8
pip3 install pandas==1.5.1 -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install networkx==2.8.8 -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install python-igraph==0.10.2 -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install leidenalg==0.9.0 -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install python-louvain==0.16 -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install scikit-learn==1.1.3 -i https://pypi.tuna.tsinghua.edu.cn/simple

pip install "numpy<2.0.0" --force-reinstall -i https://pypi.tuna.tsinghua.edu.cn/simple

在Rstudio中使用python分析空转数据：reticulate包

1
2
3
4
# 安装reticulate包
install.packages("reticulate")
# 加载R包
library(reticulate)

R中配置环境：

1
2
# 直接指定使用python版本的执行程序
use_python("/home/oelab/App/miniconda3/envs/giotto_env/bin/")

直接在Tools的global option里面选择

1
2
3
4
5
6
7
8
# 查看配置信息
py_config()

# 查看conda的版本
conda_version()

# 查看python是否可用
py_available()

启动python：

1
repl_python()

安装Giotto：

1
2
3
4
5
if(!"pak" %in% installed.packages()) {
  install.packages("pak")
}

pak::pkg_install("drieslab/Giotto")

检查是否安装成功：

1
library(Giotto)

检查环境是否能够调用成功：

1
2
3
4
5
6
7
> checkGiottoEnvironment()

Giotto can access environment found at:
 '/home/oelab/App/miniconda3/envs/giotto_env/bin/python'
 If this is the wrong environment, try specifying `envname` param
 or set option "giotto.py_path" with the desired envname or path
[1] TRUE

创建Giotto Visium对象和可视化流程

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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

library(Giotto)
library(Seurat)
library(tidyverse)
library(patchwork)
setwd("/home/oelab/RunData2/STM_AD_analyze/data/")

checkGiottoEnvironment()

## Create Giotto Instructions
python_path <- NULL
instructions <- createGiottoInstructions(save_dir = '../result/', 
                                         save_plot = TRUE,
                                         show_plot = FALSE,
                                         return_plot = FALSE,
                                         python_path = "/home/oelab/App/miniconda3/envs/giotto_env/bin/python")

## 提供visium数据所在的文件夹的路径并读取数据 
data_path <- "../data/R_GSE220442/AD1/"

## directly from visium folder
visium_brain <- createGiottoVisiumObject(visium_dir = data_path,
                                         expr_data = "filter",
                                         png_name = "tissue_lowres_image.png",
                                         gene_column_index = 2,
                                         instructions = instructions)

## show associated images with giotto object
showGiottoImageNames(visium_brain)  # "image" is the default name

## check metadata
pDataDT(visium_brain)

## 画图
spatPlot2D(gobject = visium_brain,
           cell_color = "in_tissue",
           point_size = 2,
           cell_color_code = c("0" = "lightgrey", "1" = "blue"),
           show_image = TRUE,
           image_name = "image")


## 数据处理

## 被组织覆盖的 spot 的子集
metadata <- pDataDT(visium_brain)
in_tissue_barcodes <- metadata[in_tissue == 1]$cell_ID
visium_brain <- subsetGiotto(visium_brain,
                             cell_ids = in_tissue_barcodes)

## 过滤
visium_brain <- filterGiotto(gobject = visium_brain,
                             expression_threshold = 1,
                             feat_det_in_min_cells = 50,
                             min_det_feats_per_cell = 1000,
                             expression_values = "raw",
                             verbose = TRUE)

## 归一化
visium_brain <- normalizeGiotto(gobject = visium_brain,
                                scalefactor = 6000,
                                verbose = TRUE)


## 添加基因和细胞的统计信息
visium_brain <- addStatistics(gobject = visium_brain)
## visualize 画图
## 就是测序深度的画图展示
spatPlot2D(gobject = visium_brain,
           show_image = TRUE,
           point_alpha = 0.7,
           cell_color = "nr_feats",
           color_as_factor = FALSE)

## 降维

"""

降维的顺序
1. identi fy highly variable genes (HVG)
2. perform PCA
3. identi fy number of significant prinicipal components (PCs)
4. run UMAP and/or TSNE on PCs (or directly on matrix)
5. create a shared (default) nearest network in PCA space (or directly on matrix)
cluster on nearest network with Leiden or Louvan (kmeans and hclust are alternati ves)

"""

## 找高变基因
visium_brain <- calculateHVF(gobject = visium_brain,
                             save_plot = TRUE)

## PCA降维
gene_metadata <- fDataDT(visium_brain)
featgenes <- gene_metadata[hvf == "yes" & perc_cells > 3 & mean_expr_det > 0.4]$feat_ID

## run PCA on expression values (default)
visium_brain <- runPCA(gobject = visium_brain,
                       feats_to_use = featgenes)

screePlot(visium_brain, 
          ncp = 30)

dimPlot2D(gobject = visium_brain,
          dim_reduction_to_use = "pca")


## Umap 和 tsne 可视化降维
visium_brain <- runUMAP(visium_brain,
                        dimensions_to_use = 1:10)

plotUMAP(gobject = visium_brain)

visium_brain <- runtSNE(visium_brain,
                        dimensions_to_use = 1:10)

plotTSNE(gobject = visium_brain)


# 聚类

## sNN network
visium_brain <- createNearestNetwork(gobject = visium_brain,
                                     dimensions_to_use = 1:10,
                                     k = 15)


## Leiden 聚类
visium_brain <- doLeidenCluster(gobject = visium_brain,
                                resolution = 0.4,
                                n_iterations = 1000)
## 可视化
spatDimPlot(gobject = visium_brain,
            cell_color = "leiden_clus",
            dim_point_size = 2,
            spat_point_size = 2.5)

spatDimPlot(gobject = visium_brain,
            cell_color = "nr_feats",
            color_as_factor = FALSE,
            dim_point_size = 2,
            spat_point_size = 2.5)


# dimension plots grouped by cluster
spatPlot2D(visium_brain,
           cell_color = "leiden_clus",
           coord_fix_ratio = 1)

spatPlot2D(visium_brain,
           cell_color = "leiden_clus",
           select_cell_groups = "8",
           coord_fix_ratio = 1,
           show_other_cells = TRUE,
           other_cell_color = "grey",
           other_point_size = 1.5)


## Plot with group by
## 分别展示
spatPlot2D(visium_brain,
           cell_color = "leiden_clus",
           group_by = "leiden_clus",
           coord_fix_ratio = 1,
           cow_n_col = 6,
           show_legend = FALSE)


## 数据提取

## create subset
DG_subset <- subsetGiottoLocs(visium_brain,
                              x_max = 6500, x_min = 3000,
                              y_max = -2500, y_min = -5500,
                              return_gobject = TRUE)

## show subset
spatDimPlot(gobject = DG_subset,
            cell_color = "leiden_clus",
            spat_point_size = 5)

## 找 marker 基因

## Gini markers
markers_gini <- findMarkers_one_vs_all(gobject = visium_brain,
                                       method = "gini",
                                       expression_values = "normalized",
                                       cluster_column = "leiden_clus",
                                       min_feats = 20,
                                       min_expr_gini_score = 0.5,
                                       min_det_gini_score = 0.5)

topgenes_gini <- markers_gini[, head(.SD, 2), by = "cluster"]$feats

# violinplot
violinPlot(visium_brain,
           feats = unique(topgenes_gini),
           cluster_column = "leiden_clus",
           strip_text = 8,
           strip_position = "right")

# cluster heatmap
plotMetaDataHeatmap(visium_brain,
                    selected_feats = unique(topgenes_gini),
                    metadata_cols = "leiden_clus",
                    x_text_size = 10,
                    y_text_size = 10)

## FeaturePlot 
dimFeatPlot2D(visium_brain,
              expression_values = "scaled",
              feats = markers_gini[, head(.SD, 1), by = "cluster"]$feats,
              cow_n_col = 4,
              point_size = 0.75)

## 第二种算法找marker基因---Scran Markers

markers_scran <- findMarkers_one_vs_all(gobject = visium_brain,
                                        method = "scran",
                                        expression_values = "normalized",
                                        cluster_column = "leiden_clus")

topgenes_scran <- markers_scran[, head(.SD, 2), by = "cluster"]$feats

# violinplot
violinPlot(visium_brain,
           feats = unique(topgenes_scran),
           cluster_column = "leiden_clus",
           strip_text = 10,
           strip_position = "right")

# cluster heatmap
plotMetaDataHeatmap(visium_brain,
                    selected_feats = topgenes_scran,
                    metadata_cols = "leiden_clus")

dimFeatPlot2D(visium_brain,
              expression_values = "scaled",
              feats = markers_scran[, head(.SD, 1), by = "cluster"]$feats,
              point_size = 1)


## 细胞类型注释

"""

方法分类：
1. PAGE
2. hypergeometric test
3. Rank
4. DWLS Deconvolution Corresponded Single cell dataset can be generated from here. Giotto_SC is
   processed from the downsampled Loom file and can also be downloaded from
   GiottoData::getSpatialDataset()
"""

# ## 创建单细胞数据对象
# # 假设你的Seurat对象为seurat_obj
# 
# # 提取表达矩阵
# # 使用 layer 参数提取原始计数矩阵
# 
# scRNA_harmony <- readRDS("/home/oelab/RunData2/AD_Analyze/scRNA_harmony.rds")
# 
# scRNA_harmony <- JoinLayers(scRNA_harmony)
# 
# sc_expression <- as.matrix(GetAssayData(scRNA_harmony, layer = "counts"))
# 
# # 提取元数据
# sc_metadata <- scRNA_harmony@meta.data
# 
# giotto_SC <- createGiottoObject(expression = sc_expression,
#                                 instructions = instructions)
#                             
# giotto_SC <- addCellMetadata(giotto_SC,
#                              new_metadata = data.table::fread(sc_metadata))
#                             
# giotto_SC <- normalizeGiotto(giotto_SC)


## PAGE enrichment 方法


## 1.1 用细胞signature genes创建二进制矩阵

Astrocyte_Markers <- c("GFAP", "AQP4", "SLC1A3", "S100B", "SLC1A2", 
                       "APOE", "ALDH1L1", "SOX9")

OPC_Markers <- c("ABRACL", "OLIG2", "CSPG4", "OLIG1", 
                 "VCAN", "LHFPL3", "GPR17", "APOD", "SOX6")

Olig_Markers <- c("PLP1", "MOG", "MOBP", "MAG", "CLDN11", "MBP", 
                  "OLIG1", "OLIG2", "CNP")

Neuron_Markers <- c("SNHG6", "DCX", "SNAP25", "SYT1", "GAP43",
                    "MAP2", "RBFOX3", "TUBB3", "SATB2")

Micro_Markers <- c("HAVCR2", "CSF3R", "C1QB", "C1QC", 
                   "CSF1R", "CX3CR1", "CD14", "OLR1", "CD74")

Endo_Markers <- c("TEK", "EMCN", "CLDN5", "FLT1",
                  "CLEC14A", "ADGRL4", "ADGRF5", "PECAM1", "ABCG2")


PAGE_matrix_1 <- makeSignMatrixPAGE(sign_names = c("Astrocyte",
                                                   "OPC",
                                                   "Olig.",
                                                   "Neuron", "Micro.", "Endo."),
                                    sign_list = list(Astrocyte_Markers,
                                                     OPC_Markers,
                                                     Olig_Markers, 
                                                     Neuron_Markers, 
                                                     Micro_Markers, 
                                                     Endo_Markers))


## 方法一：enrichment test with PAGE
visium_brain <- runPAGEEnrich(gobject = visium_brain,
                              sign_matrix = PAGE_matrix_1)

## heatmap of enrichment versus annotation (e.g. clustering result)

cell_types_PAGE <- colnames(PAGE_matrix_1)

plotMetaDataCellsHeatmap(gobject = visium_brain,
                         metadata_cols = "leiden_clus",
                         value_cols = cell_types_PAGE,
                         spat_enr_names = "PAGE",
                         x_text_size = 8,
                         y_text_size = 8)
## 可视化注释结果

spatCellPlot2D(gobject = visium_brain,
                  spat_enr_names = "PAGE",
                  cell_annotation_values = cell_types_PAGE[1:6],
                  cow_n_col = 2,
                  coord_fix_ratio = 1,
                  point_size = 1.25,
                  show_legend = TRUE)

spatDimCellPlot2D(gobject = visium_brain,
                  spat_enr_names = "PAGE",
                  cell_annotation_values = cell_types_PAGE[1:6],
                  cow_n_col = 1,
                  spat_point_size = 1,
                  plot_alignment = "horizontal",
                  save_param = list(base_width = 7, base_height = 10))


## 方法二：HyperGeometric test 超几何检验

visium_brain <- runHyperGeometricEnrich(gobject = visium_brain, 
                                        sign_matrix = PAGE_matrix_1)

cell_types_HyperGeometric <- colnames(PAGE_matrix_1)

spatCellPlot(gobject = visium_brain,
             spat_enr_names = "hypergeometric",
             cell_annotation_values = cell_types_HyperGeometric[1:6],
             cow_n_col = 2,
             coord_fix_ratio = NULL,
             point_size = 1.75)


## 方法三：Rank Enrichment，得用giotto_SC

# rank_matrix <- makeSignMatrixRank(sc_matrix = getExpression(giotto_SC,
#                                                             values = "normalized", 
#                                                             output = "matrix"),
#                                   sc_cluster_ids = pDataDT(giotto_SC)$Class)
# 
# colnames(rank_matrix) <- levels(factor(pDataDT(giotto_SC)$Class))
# 
# visium_brain <- runRankEnrich(gobject = visium_brain,
#                               sign_matrix = rank_matrix,
#                               expression_values = "normalized")


## 方法四：DWLS 空间反卷积

# DWLS_matrix <- makeSignMatrixDWLSfromMatrix(
#   matrix = getExpression(giotto_SC,
#                          values = "normalized",
#                          output = "matrix"),
#   cell_type = pDataDT(giotto_SC)$Class,
#   sign_gene = topgenes_scran$feats)
# 
# visium_brain <- runDWLSDeconv(gobject = visium_brain,
#                               sign_matrix = DWLS_matrix)

## Spatial Grid--空间网格

visium_brain <- createSpatialGrid(gobject = visium_brain,
                                  sdimx_stepsize = 400,
                                  sdimy_stepsize = 400,
                                  minimum_padding = 0)
showGiottoSpatGrids(visium_brain)
spatPlot2D(visium_brain,
           cell_color = "leiden_clus",
           show_grid = TRUE,
           grid_color = "red",
           spatial_grid_name = "spatial_grid")


## Spatial network---空间网络
## 创建一个空间网络，连接单个细胞基于它们彼此的物理距离

visium_brain <- createSpatialNetwork(gobject = visium_brain,
                                     method = "kNN",
                                     k = 5,
                                     maximum_distance_knn = 400,
                                     name = "spatial_network")
showGiottoSpatNetworks(visium_brain)
spatPlot2D(gobject = visium_brain,
           show_network= TRUE,
           network_color = "blue",
           spatial_network_name = "spatial_network")

## Spatial Genes---空间基因
## 如果一个基因通常被发现在近端或邻近细胞中高度表达，那么它就被认为是一个空间基因

ranktest <- binSpect(visium_brain,
                     bin_method = "rank",
                     calc_hub = TRUE,
                     hub_min_int = 5,
                     spatial_network_name = "spatial_network")

spatFeatPlot2D(visium_brain,
               expression_values = "scaled",
               feats = ranktest$feats[1:6],
               cow_n_col = 2,
               point_size = 1.5)


## 空间共表达模块
## 1. cluster the top 500 spatial genes into 20 cluster
my_spatial_genes <- ranktest[1:1500,]$feats

## 2. 这里我们使用现有的detectSpatialCorGenes函数来计算基因之间的成对距离(但将network_smoothing=0设置为使用默认聚类)

spat_cor_netw_DT <- detectSpatialCorFeats(visium_brain,
                                          method = "network",
                                          spatial_network_name = "spatial_network",
                                          subset_feats = my_spatial_genes)

## 确定一个基因最相似的空间相关基因
top10_genes <- showSpatialCorFeats(spat_cor_netw_DT, 
                                   feats = "Mbp", 
                                   show_top_feats = 10)

## 可视化
spatFeatPlot2D(visium_brain,
               expression_values = "scaled",
               feats = top10_genes$variable[1:4],
               point_size = 3)

## 空间基因的聚类
## 类似WCGNA

## cluster spatial genes
spat_cor_netw_DT <- clusterSpatialCorFeats(spat_cor_netw_DT,
                                           name = "spat_netw_clus",
                                           k = 20)

## visualize clusters
heatmSpatialCorFeats(visium_brain,
                     spatCorObject = spat_cor_netw_DT,
                     use_clus_name = "spat_netw_clus",
                     heatmap_legend_param = list(title = NULL))

## 查看每个clusters的基因
netw_ranks = rankSpatialCorGroups(visium_brain, spatCorObject = spat_cor_netw_DT,
                                  use_clus_name = 'spat_netw_clus')
cluster_genes_DT = showSpatialCorGenes(spat_cor_netw_DT,
                                       use_clus_name = 'spat_netw_clus',
                                       show_top_genes = 1)
cluster_genes = cluster_genes_DT$clus
names(cluster_genes) = cluster_genes_DT$feat_ID

## rank spatial correlated clusters and show genes for selected clusters
netw_ranks <- rankSpatialCorGroups(visium_brain,
                                   spatCorObject = spat_cor_netw_DT,
                                   use_clus_name = "spat_netw_clus")
top_netw_spat_cluster <- showSpatialCorFeats(spat_cor_netw_DT,
                                             use_clus_name = "spat_netw_clus",
                                             selected_clusters = 6,
                                             show_top_feats = 1)

## create metagene enrichment score for clusters
## 查看上述结果每个模块在空间的表达情况
cluster_genes_DT <- showSpatialCorFeats(spat_cor_netw_DT,
                                        use_clus_name = "spat_netw_clus",
                                        show_top_feats = 1)
cluster_genes <- cluster_genes_DT$clus
names(cluster_genes) <- cluster_genes_DT$feat_ID
visium_brain <- createMetafeats(visium_brain,
                                feat_clusters = cluster_genes,
                                name = "cluster_metagene")
spatCellPlot(visium_brain,
             spat_enr_names = "cluster_metagene",
             cell_annotation_values = netw_ranks$clusters,
             point_size = 1,
             cow_n_col = 5)

R中单细胞空转联合分析 - RCTD

单细胞联合空间主流软件介绍

R中常用的软件：

  • Seurat：TransferAnchors
  • RCTD：robust cell type decomposition 鲁棒细胞类型分解
  • Spotlight：非负矩阵因子分解回归
  • CellTrek：训练多元随机森林（RF）模型，使用共享降维特征预测空间坐标
  • CARD：conditional autoregressive-based deconvolution，基于条件自回归的反卷积
  • Giotto：PAGE（Parametric Analysis of Gene Set Enrichment）
  • Giotto：HyperGeometric test 超几何检验（与MIA相似）
  • Giotto：rank，它不是基于基因集的方法，而是使用单细胞所有的表达数据
  • MIA：Multimodal intersection analysis（多模态交集分析）

python中常用的软件：

  • Scanorama
  • Tangram
  • Cell2location

RCTD

RCTD：一种利用从scRNA-seq中学习细胞类型的表达谱对空间转录组 spot 内混合细胞的类型进行稳健分解的工具 RCTD（robust cell type decomposition），并考虑到了来自不同测序方法数据间的系统误差，克服了平台效应（platform effect）的影响

教程：https://github.com/dmcable/spacexr

RCTD可以将单细胞类型或细胞类型混合分配到空间转录组spots上。RCTD 有三种模式：

  1. Doublet mode：它为每个Spot分配1-2种细胞类型，推荐用于高空间分辨率的技术，如Slide-seq和MERFISH
  2. Full mode：它为每个Spot分配任意数量的细胞类型，推荐用于空间分辨率较低的技术，如100微米分辨率的Visium
  3. Multi mode：doublet mode的扩展，可以在每个点上发现两种以上的细胞类型，作为full mode的替代选项

模式一：Doublet mode

加载数据包：

1
2
3
4
# install.packages("devtools")
devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)
library(spacexr)
library(Matrix)

R中单细胞空转联合分析 - CellTrek（了解）

CellTrek：可以结合单细胞和空间转录组数据来准确地定位组织内单个细胞的位置，并构建空间细胞图谱

  • 一种可以基于scRNA-seq和ST数据将单细胞直接映射回组织切片空间坐标的计算工具包。该方法不同于ST解卷积，能够更灵活和直接地研究具有空间地形的单细胞数据
  • CellTrek工具包还提供了两个下游分析模块，包括用于空间共定位分析的SColoc和用于空间共表达分析的SCoexp