具体实施方式

下面结合具体实施例对本发明作进一步解释，但是应当注意的是，以下实施例仅用以解释本发明，而不能用来限制本发明，所有与本发明相同或相近的技术方案均在本发明的保护范围之内。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段，所用原料为市售商品。

所述组织透化液、荧光标签均选自试剂盒，试剂盒名称为Visium Spatial GeneExpression Reagent Kits-Tissue Optimization，来自10X genomics，可购自上海仁科生物科技有限公司，货号PN-100 0193；所述文库构建时的PCR扩增过程采用试剂盒扩增，试剂盒名称为Visium Spatial Gene Expression Reagent Kits，来自10X genomics，可购自上海仁科生物科技有限公司，货号PN-1000184，扩增程序为98℃3min，(98℃5s，63℃30s)25个循环。

实施例一种空间转录组测序数据的分析方法

所述空间转录组数据生物信息分析方法，包括以下步骤：

S1、样本检测建库：

(1)样本准备

将组织速冻后，使用OCT包埋剂对组织进行包埋。在低温环境下，使用冷冻切片技术从组织块上切取一张长、宽各10mm，10μm厚的切片分别贴附到组织优化芯片或文库制备芯片上。然后使用预冷甲醇固定切片，苏木精-伊红染色法(HE染色)对切片进行染色。组织染色结果通过明场显微镜记录组织形态和spot位置。

样品准备是整个技术的起始，其难点就在于冷冻切片。空间转录组的芯片单spot直径是55μm，并没有达到真正的单细胞级，单点捕获的真实细胞数量是由切片厚度、切片平整度和细胞直径共同决定的。

不同组织的细胞大小是不一样的。在普通组织切片中，虽可以追求极致的单层细胞而选择2-3μm的切片厚度。但是在空间转录组，2-3μm意味着一层破碎的细胞，是不适合进行下游实验的。所以，切片的厚度也应随着细胞大小而做出调整。

不同组织的脂质、蛋白质含量是不一样的，导致在相同的冷冻温度下其硬度和延展性也存在差异。对于一块组织来说，温度过低，组织脆性会增加，获得的切片会有或大或小的裂隙；温度过高，组织延展性会降低，切片容易出现褶皱。所以，切片环境温度会成为组织相关的个性化点。

(2)组织优化

使用组织优化芯片对组织透化时间进行优化。首先，设置梯度透化时间(一般为5min，10min，15min，20min，25min，30min的时间梯度)向完成了明场成像的芯片添加组织透化液；然后在透化完成的芯片上进行荧光反转录，将成功捕获的mRNA反转录为带荧光标签的cDNA；反转录完成后，去除掉芯片上的组织以避免切片对荧光信号的影响；最后，使用荧光显微镜对芯片进行成像，依据荧光成像结果的亮度、弥散程度和组织透化时间选取最佳组织透化时间用于文库构建。

优选的，最佳组织透化时间通过三个条件决定：(1)荧光信号最强，(2)荧光信号弥散最弱，(3)前两个情况一致的情况下选择较长的组织透化时间。

(3)文库制备

使用文库制备芯片进行空间转录组文库制备。将完成了明场成像的芯片按照最佳组织透化时间进行组织透化，释放组织mRNA与芯片上的捕获序列结合，逆转录合成cDNA片段并标记空间位置，PCR扩增合成cDNA二链。通过碱性环境使cDNA变性将cDNA二链释放到液体游离环境，以cDNA二链为模板进一步扩增合成大量cDNA。将cDNA酶切打断成200-300bp的片段，加上P5接头、i5样本索引、read 2、i7样本索引和P7接头构建标准测序文库。

(4)文库测序

利用Illumina测序平台的双端测序模式对建好的文库进行高通量测序。在Read 1端，包含16bp的barcode信息和12bp的UMI信息用于确定转录本位置和表达量定量；在Read2端，包含了cDNA片段，用于参考基因组比对确定mRNA所对应的基因。

S2、测序数据生物信息分析

(1)测序数据质控及基因表达量定量：测序获得的原始数据、组织染色照片和芯片序列号一起输入space ranger进行数据质控和序列比对，以获得高质量的测序数据、spot在组织中的空间分布位置和基因在各个spot的表达量情况。具体包括如下步骤：

①数据测序的基本质控

将测序原始数据输入space ranger中，对spot位置进行比对获得覆盖在组织切片下的spot作为有效spot，结合芯片的序列号获得有效spot对应的barcode信息。基于有效barcode信息，对测序数据的reads质量、测序饱和度和有效barcode进行评估。

Space Ranger依靠图像处理算法来解决与载玻片图像有关的两个关键问题：决定组织被放置在哪里，以及对准校准点图案。通过组织检测来确定哪些捕获点被组织覆盖，从而确定哪些barcode将用于分析。校准点的对齐是必需的，Space Ranger就能知道每个spot在像中的位置。校准点的对齐首先提取"看起来"像校准点的特征，然后尝试将这些候选校准点与已知的校准点模式对齐。从图像中提取出来的spot必然会包含一些遗漏的地方，例如，在组织覆盖的地方，以及一些假阳性的地方，例如载玻片上的碎片或组织特征可能看起来像spot。这个过程的输出是一个spot坐标值表格，用于将Visium spot位置与组织图像相关联。组织检测时使用灰度的组织图像。计算并比较组织切片位置的多个估计值。这些估计值用于训练一个统计分类器，将捕获区域内的每个像素分别标记为组织或背景，从而判定哪些区域属于组织。

②参考基因组比对和基因组表达定量

将reads在space ranger中与参考基因组进行比对，对所有spot的基因数、有效reads数等参数进行统计评估。根据基因组比对结果、UMI计数和barcode序列信息，获得不同基因在不同的spot中的表达量丰度信息。该矩阵可用于后续分析的进行。

Space Ranger调用STAR(Spliced Transcripts Alignment to a Reference)比对软件将Read2比对到参考基因组。依据GTF注释的结果，将比对上的reads分类为外显子、内含子或基因间区的reads：如一条reads有至少50％位于外显子区则属于外显子reads，以此类推。如果reads比对到单一的外显子区域和一个/多个非外显子区，则优先归类为外显子且它的MAPQ被置为255。Space Ranger进一步将外显子reads比对到已注释的转录本上：1.如果reads比对到某条转录本的外显子且两者方向相同，则认为它比对到了转录本上，记为转录本reads；2.在转录组reads中，如果reads只比对到一个基因，则认为它们是唯一比对reads。只有唯一比对reads才用于UMI计数(表达量定量)。

Space Ranger接着会校正UMI测序错误。具有相同的barcode、UMI和基因注释的reads会被归为同一组。如果两组reads具有相同的barcode序列和基因注释，但是UMI具有一个碱基的差异，那么其中一条UMI序列可能存在测序错误。那么，reads数较低的UMI会被修正为reads数较高的一组UMI的序列。基于修正后的UMI序列，再次将具有相同的barcode、UMI和基因注释的reads归为同一组。如果两组或两组以上的reads具有相同的barcode序列和UMI序列但是不同的基因注释，将reads数最高的一组reads用于UMI计数，其余组抛弃；当出现多组reads的reads数都很高时，将所有reads组都排出UMI计数范围外。在完成以上矫正和过滤后，将每个barcode的每个基因对应的UMI去重，计算unique UMI的数量做为该spot中该基因的表达量。

(2)spots聚类分析

空间转录组中每个spot包含1-10个细胞，利用spot分子特征的相似性对spots进行聚类分群，获得的spots分群结果可与spots对应组织类型形成对应，辅助空间异质性解析。后续，将space ranger从得到的表达量信息、spot空间位置信息转入seurat进行后续分析。

对表达量数据进行SCTransform均一化处理；然后，利用均一化的表达量进行PCA(Principal components analysis)分析，减少背景噪音对聚类和降维的影响；最后，采用基于图论的聚类算法对spot进行聚类和分群。基于spot亚群分类结果，进一步使用tSNE(t-Distributed Stochastic Neighbor Embedding)非线性降维的方法将所有spot投影至二维平面中可视化。具体过程如下：

①数据归一化

在进行spots聚类之前，首先对基因表达量进行均一化处理。使用SCT ransform算法对基因表达量进行均一化处理，该算法将基因表达量的对数值放入一个符合负二项分布的回归模型中。在该回归模型中，基因i的表达量应符合如下公式：log(E(x_i))＝β₀+β₁lg₁₀m，其中，x_i是基因i在各个spot中的表达量，m是各个spot中的UMI总和。在该模型中，也相应计算出每个基因表达量的离散系数θ和表达量均值μ。使用负二项分布对模型参数进行矫正后，使用表达量的数学期望值来代替原有的表达量。表达量的数学期望值公式为：μ_ij＝exp(β_0i+β_1ilog₁₀m_j),表达量归一化的计算公式为：

在这两个公式中，x_ij为基因i在spot j中的实际表达量，m_j为spot j中所有UMI数的总和，μ_ij为基因i在spot j中的表达量数学期望，θ_ij为基因i在spot j中的离散系数。

数据均一化的目的，主要是降低测序深度对基因表达量的影响，若未进行归一化的数据，基因表达量随着测序深度提高而提高，测序深度高的细胞对表达量变异程度贡献高。传统常规的log均一化对低丰度基因校正能力强，对高丰度基因较正能力差，所以log均一化的数据高丰度基因表达量随着测序深度的提高而提高，测序深度低的细胞对高丰度基因的变异度贡献被高估。本发明采用，SCTransform均一化，SCTransform归一化的数据，除了丰度最高的一组基因，其余基因的表达丰度都不会随着测序深度变化，不同测序深度的细胞对基因丰度变异度贡献基本一致。所以SCTransform归一化相对于log均一化表现出更真实的差异基因检验结果。

②主成分分析(Principal component analysis，PCA)

为了克服单个基因表达量对空间转录组数据的噪音干扰，Seurat通过PCA将spot进行聚类，每个主成分代表表达模式相近的一个基因集。确定下游分析使用多少个主成分关系到后续分析的背景噪音和精确度。为了确定下游分析使用多少个主成分，使用重复取样检验法检验各主成分对空间转录组数据的影响程度。将检验结果中富集最多的低P值基因的主成分作为显著影响空间转录组数据集的主成分，用于下游分析。

③spot聚类

Seurat使用基于图论的方法对spot进行聚类。首先，基于之前确定的主成分，计算spot之间的欧氏距离；然后，基于spot间的欧式距离将spot嵌入KNN(K-nearest neighborgraph)图中；接着，根据两个spot在局部相邻的Jaccard距离将KNN图分割成高度互连的准群体；最后，使用Louvain算法将spot聚类分群。以此，将空间转录组中的所有spot聚类划分为若干spot亚群，便于后续分析的进行。

④tSNE可视化

Seurat使用非线性降维法t-SNE将高维的空间转录组数据映射到二维空间中，将表达模式相似的spot聚集在一起，将表达模式不同的spot分隔得更远，更直观地展现spot间差异性。

(3)差异基因分析

①亚群差异表达分析

使用Seurat的MAST检验(Model-based Analysis of Single cellTranscriptomics)分别对不同spot亚群进行亚群上调基因分析，并按照以下阈值来筛选显著上调基因：

a.基因上调倍数的自然对数logFC≥0.25，即基因上调倍数≥e ^0.25≈1.28；

b.P值≤0.01；

c.在目标亚群中，有25％以上的spot表达目标基因；

MAST检验具有更适用于复杂实验设计，以及对表达基因的细胞比例做拟合，降低测序误差影响的优点。

②亚群差异基因功能分析

对筛选的亚群上调基因进行KEGG富集分析，DO富集分析，GO富集分析，Reactome富集分析等分析。

③亚群特征基因筛选和分析

进一步挑选每个亚群表达量上调倍数最高的20个基因做为亚群特征基因，并以热图、组织映射图和气泡图来展示基因分布特征。

热图：使用SCTransform法对基因表达量进行归一化处理，然后将归一化后的基因表达量用于热图绘制。如此，避免热图中极端值的出现，使得基因在不同亚群的表达特征更加明显。

组织映射图：使用SCTransform均一化法对基因表达量进行处理，减小基因在不同spot间的表达差异性；然后将均一化后的基因表达量用于组织映射图的绘制。如此，在保留一定原始表达量特征的情况下减小了极值对数据分布的影响，便于观察基因在不同组织区域间的表达特征。

气泡图：对每个基因在每个细胞亚群的未取log值的均一化表达量取均值，然后基于表达量均值在不同亚群间对基因平均表达量做SCTransform归一化，使用归一化的基因表达量绘制气泡图中气泡的颜色深浅。另外，统计表达目标基因的细胞在亚群中占比，用于绘制气泡图中气泡的大小。气泡图同时展示了基因在细胞亚群中的平均表达量和表达目标基因细胞的比例。

④亚群特征基因蛋白质互作网络分析

通过string数据库(https://string-db.org/)，获得蛋白质间互作信息。通过这些信息，建立亚群特征基因间的蛋白质互作网络，进一步了解不同spot亚群之间存在的潜在的调控关系。

如背景技术所述，空间转录组的数据信息是很庞大，每个spot点除了基因丰度信息外，还带有空间位置信息。在上述基础分析中，将spot进行分群的分析方式仍并未将空间位置信息充分利用起来。所以，需要进一步进行分析方法，尽可能多地利用空间位置信息，充分挖掘特征基因信息。

(4)空间特征基因分析

①空间差异基因分析

为了深入探索表达量对空间位置依赖的基因-空间差异基因，首先通过Seurat的markvariogram算法对高变基因进行检验，再通过trendsceek的mark-segregation假设检验对基因显著性进行检验，最后通过BH法对显著性P值进行多重检验矫正。通过以下条件筛选空间差异基因：

a.r.metric.5≤0.8

b.P值≤0.01

本发明采用的markvariogram算法，具有不依赖spot分群结果和组织区域鉴定结果，可以更充分的挖掘与空间位置相关的变异基因的优点，同时也可进行显著性检验。

②空间差异基因功能分析

对筛选的空间差异基因，进行KEGG富集分析，DO富集分析，GO富集分析，Reactome富集分析等分析。

③空间特征基因筛选和分析

进一步挑选P值最小的200个基因作为空间特征基因，并以热图、组织映射图和气泡图来展示基因分布特征。具体的：

热图：使用SCTransform法对基因表达量进行归一化处理，然后将归一化后的基因表达量用于热图绘制。如此，避免热图中极端值的出现，使得基因在不同亚群的表达特征更加明显。但是，对于空间特征基因来说，其分布并不完全依赖于spot亚群，所以热图只能显示出部分基因的分布特征。

组织映射图：使用SCTransform均一化法对基因表达量进行处理，减小基因在不同spot间的表达差异性；然后将均一化后的基因表达量用于组织映射图的绘制。如此，在保留一定原始表达量特征的情况下减小了极值对数据分布的影响，便于观察基因在不同组织区域间的表达特征。空间特征基因因为是在空间分布中分析得到的基因，所以组织映射图可以最直观得体现空间特征基因的分布特征。

气泡图：气泡图同时展示了基因在细胞亚群中的平均表达量和表达目标基因细胞的比例。为了绘制气泡图，对每个基因在每个细胞亚群的未取log值的均一化表达量取均值，然后基于表达量均值在不同亚群间对基因平均表达量做SCTransform归一化，使用归一化的基因表达量绘制气泡图中气泡的颜色深浅。另外，统计表达目标基因的细胞在亚群中占比，用于绘制气泡图中气泡的大小。与热图类似，气泡图也只能显示出部分空间特征基因的表达特征。

④空间特征基因蛋白质互作网络分析

通过string数据库，获得蛋白质间互作信息。通过这些信息，建立空间特征基因间的蛋白质互作网络，进一步了解不同组织区域之间存在的潜在的调控关系。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明做了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。