一种空间转录组与空间代谢组点对点匹配方法及系统

技术领域

本发明属于图像数据处理技术领域，涉及一种空间转录组与空间代谢组点对点匹配方法及系统。

背景技术

空间转录组测序是单细胞转录组技术的深入发展成果，可以同时获得细胞的空间位置信息和基因表达数据，且不需要制备细胞悬液，进一步推进了对组织原位细胞真实基因表达的研究，为组织细胞功能、微环境互作、发育过程谱系追踪、疾病病理学等多个领域提供了重要的研究手段。

空间转录组作为一项紧跟着高通量单细胞转录组技术发展起来的新兴前沿技术，其能够保留组织空间位置并同时解析组织切片转录组信息的技术。空间转录组是一种用于从空间层面上解析RNA-seq数据的技术，从而解析单个组织切片中的所有mRNA。空间条形码逆转录引物在芯片表面有序附着，使得样品处理和后续测序过程中空间位置信息的编码和获取成为可能。当组织冷冻切片附在空间转录组芯片上时，条形码引物结合并从组织中捕获邻近的mRNA。被捕获的mRNA开始逆转录，得到的cDNA中包含了空间条形码。分析测序结果中空间条形码的序列，可以将每个mRNA转录的序列映射回组织切片中的起始位置。

空间代谢组学是质谱成像技术的一种，能够直接从生物组织中获得大量已知或未知的内源性代谢物和外源性药物等分子的结构、含量和空间分布信息。无需化学或放射性标记、不需复杂样品前处理，具有高特异性、高通量和空间信息保留的突出优势。

空间代谢组学利用质谱成像结合代谢组学的方法，充分利用质谱成像的准确识别并定位多种代谢物在组织、甚至细胞间的差异性分布的能力，结合代谢组学可对目标微区域组织进行深度代谢组学分析并获得代谢物种类和含量的能力，实现了对生物学样本的代谢物空间分布的检测。空间代谢组学能够在无标记条件下直接对生物组织切片进行扫描分析，通过扫描点的扫描信息进行空间位置定位。

转录组是获得生物体内基因表达的重要方法，代谢组是生物体表型的基础和直接体现者。代谢产物是生物体在内外调控下基因转录的最终结果，是生物体表型的物质基础。在系统生物学研究时代，生物发育过程复杂多变，基因调控网络复杂，单独用一种组学对系统生物学研究，往往比较片面。利用空间转录组与空间代谢组进行关联分析，从原因和结果两个层次分析生物体内在的变化，锁定与代谢物变化重点通路，构建核心调控网络，全面探究生物的生长发育及胁迫机制，从整体上解释生物学问题。

但由于空间转录组与空间代谢组不同的空间信息标识、不同切片样本实验操作步骤及其非相同切片样本、不同空间分布的情况下，两套不同体系的数据很难进行空间位置的配准、点对点数据的匹配，更难进行空间转录组与空间代谢组的关联分析。

具体地，由于空间转录组与空间代谢组的检测平台与实验操作的不同，两套的空间分布是不同的。空间转录组以芯片作为检测平台，芯片检测区域为6.5mm×6.5mm，以spot作为检测的最小单元，单个spot直径为55μm，spot与spot中心之间的距离100μm，上下两层成交错排布。空间代谢组以质谱作为检测平台，其检测区域未作限制，最大可达10cm×10cm，以像素点作为检测的最小单元，单个像素点常见大小为100μm×100μm，像素点以正方形组合分布。空间转录组与空间代谢组的空间分布不同，导致spot与像素点之间难以进行关联，更难将spot和像素点内的底层数据进行关联。

发明内容

为了解决现有技术存在的不足，本发明的目的是提供一种空间转录组与空间代谢组点对点匹配方法及系统。

本发明提出了一种空间转录组与空间代谢组点对点匹配方法，通过将空间转录组与空间代谢组的空间信息进行匹配并对数据进行关联，实现空间转录组与空间代谢组的点对点匹配。

包括如下步骤：

步骤一、将带有标识点位的包含空间转录组信息的HE染色图与空间代谢组成像图配准，获得标识点位在空间代谢组上的空间坐标位置信息；

所述标识点位均位于空间转录组的样本检测范围内，所述标识点位的数量≥6个；所述标识点位中包含条形码barcode信息，表示标识点位的空间信息。

HE染色图和空间代谢组成像图进行配准后，HE染色图和空间代谢组成像图匹配，获得标识点位在空间代谢组成像图上的空间坐标信息。

步骤二、将空间转录组和空间代谢组的空间信息转化为统一的空间信息标识；

通过将空间转录组中条形码barcode对应的空间信息中x和y两组空间编号和空间代谢组空间信息中的x和y两组空间编号分别进行运算，转换成微米的坐标体系使空间转录组和空间代谢组的空间坐标体系统一。

步骤三、计算空间转录组与空间代谢组之间的最终缩放比和最终旋转角度；

计算获得各标志点位两两之间的缩放比，并对获得的缩放比数据进行筛选，剔除不符合要求的数据，并对剩余的缩放比数据求均值获得用于空间代谢组转换的最终缩放比。

将空间转录组中的其中任意两个标识点位的连线与空间代谢组中的两个对应的标识点位的连线计算它们之间的旋转角度；根据排列组合原理计算完全部任意两点连线的旋转角度，并将所有获得的旋转角度求均值，获得最终旋转角度。

步骤四、使用最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，使转化后的空间代谢组坐标与空间转录组坐标匹配；

步骤五、将空间代谢组的数据进行拟合，并与空间转录组的点spot进行关联；

通过将空间代谢组中的像素点扩增后，根据空间代谢组拟合后的像素点中心位置与空间转录组的点spot中心位置之间的距离进行像素点的归属。

步骤六、将每个点spot对应的像素点底层数据进行加和运算，获得空间转录组每个点spot对应的空间代谢组数据。

将归属到同一个点spot上的多个像素点的底层数据加和，获得空间转录组每个点spot对应的空间代谢组数据。

本发明还提出了实现上述空间转录组与空间代谢组点对点匹配方法的系统，及上述点对点匹配方法在空间转录组与空间代谢组的关联分析中的应用。

所述空间转录组与空间代谢组点对点匹配系统包括空间信息获取模块、空间坐标统一模块、空间参数运算模块、空间信息转化模块、空间代谢数据拟合模块、空间转录与空间代谢数据关联模块；

所述空间信息获取模块，将带有标识点位的包含空间转录组信息的HE染色图与空间代谢组成像图配准，获得标识点位在空间代谢组上的空间坐标位置信息；

所述空间坐标统一模块，将空间转录组和空间代谢组的空间信息转化为统一的空间信息标识；

所述空间参数运算模块，计算空间转录组与空间代谢组之间的最终缩放比和最终旋转角度；

所述空间信息转化模块，使用最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，使转化后的空间代谢组坐标与空间转录组坐标匹配；

所述空间代谢数据拟合模块，将空间代谢组的数据进行拟合，并与空间转录组的点spot进行关联；

所述空间转录与空间代谢数据关联模块，将每个点spot对应的像素点底层数据进行加和运算，获得空间转录组每个点spot对应的空间代谢组数据。

点对点匹配方法中，空间信息匹配通过空间转录组与空间代谢组空间信息匹配方法进行，所述空间信息匹配方法通过将空间转录组与空间代谢组的空间信息标识统一，并将空间代谢组的空间坐标体系匹配转换成空间转录组的空间坐标，实现空间转录组与空间代谢组的空间信息匹配，包括如下步骤：

步骤(1)、在空间转录组的样本检测范围内选择均匀分布(在整个样本检测范围内四周及中间进行选择分布，使选择点位扩展到整个样本检测范围)的6个及以上的点spot作为标识点位，获取标识点spot对应的条形码barcode信息并在对应HE染色图上标识选择的点spot；

所述HE染色是组织学、病理学教学与科研中最基础、使用最广泛的技术方法之一，通过HE染色呈现的病理特征与空间转录组进行关联。

所述空间转录组的表达玻片上共有4个空间转录组芯片；所述每个空间转录组芯片的检测区域大小为6.5mm×6.5mm，共4992个点spot，单个spot直径为55μm，每相邻两个spot与spot中心之间的距离为100μm；所述spot为空间转录组中检测的最小单位。

所述条形码barcode为在Visium空间转录组技术中的空间条形码，每个点spot中有其独有的barcode编码，能通过空间条形码获取对应点spot的空间位置信息。

进行HE染色后的样本依旧能够进行空间转录组检测，不会对空间转录组数据进行影响，所述HE染色图能够与空间转录组数据进行空间位置的匹配，在HE染色图上能够进行spot标识。

步骤(2)、使用步骤一中已经标识点位的HE染色图与空间代谢组成像图进行配准，获取配准后标识点位在空间代谢组上的对应像素点的空间坐标位置信息；

所述像素点是空间代谢组中检测的最小单元；所述空间代谢组中检测的最小单元的分辨率根据质谱扫描速率及移动平台速率而不同，速率越大其扫描分辨率越大。

所述空间代谢组中检测的最小单元的分辨率代表像素点与像素点之间的距离；分辨率不同，像素点与像素点之间的距离不同；在空间代谢组中空间信息的转换受到分辨率大小的影响。

所述配准操作是指以空间代谢组成像图为基准，通过MSIReader软件将HE染色图像进行旋转、缩放、翻转操作，将空间代谢组成像图和HE染色图匹配，通过软件获取HE染色图上标识位点在空间代谢组成像图上的空间坐标信息。

步骤(3)、将空间转录组的条形码barcode空间信息和空间代谢组的空间信息转化为统一的空间信息标识；

所述转化为统一的空间信息标识是指将空间转录组中条形码barcode对应的空间信息中x和y两组空间编号和空间代谢组空间信息中的x和y两组空间编号分别进行运算，以统一标准(即将空间转录组和空间代谢组的坐标均转换成微米的坐标体系)获取统一空间坐标体系。

所述空间转录组中x空间编号的转化公式如下：

UTX＝(TransX/2+0.5)*transresolution；

空间转录组中y空间编号的转化公式如下：

UTY＝TransY*sqrt(0.75)*transresolution；

其中，UTX、UTY分别为空间转录组数据统一标准后的空间坐标；TransX、TransY分别为空间转录组中x和y的两组空间编号；transresolution为空间转录组分辨率的数值；sqrt为开根号运算。

所述空间代谢组中x空间编号的转化公式如下：

UMX＝MetaX*metaresolution；

空间代谢组中y空间编号的转化公式如下：

UMY＝MetaY*metaresolution；

其中，UMX、UMY分别为空间代谢组数据统一标准后的空间坐标；MetaX、MetaY分别为空间代谢组中x和y的两组空间编号；metaresolution为空间代谢组分辨率的数值。

步骤(4)、计算空间转录组中不同标识点间距离；计算空间代谢组不同标识点间距离；并将空间转录组与空间代谢组对应的标识点间距离求缩放比；

为了进行空间转录组和空间代谢组中的标识点对应准确性评判以及校正实验仪器误差导致标志点位置差异，进行所述缩放比的计算；

所述缩放比的计算公式如下：

Ratio＝sqrt((UTXa-UTXb)

其中，UTXa、UTYa、UTXb、UTYb分别为空间转录组中标志点a和b经统一空间坐标后的X，Y坐标值；UMXa、UMYa、UMXb、UMYb分别为空间代谢组中标志点a和b经统一空间坐标后的X，Y坐标值；Ratio为空间转录组中标志点a和b的距离与空间代谢组中标志点a和b的距离的缩放比。

获得各标志点两两之间的缩放比后，还需要进行对缩放比数据的筛选，所述筛选操作包括：剔除缩放比结果大于1.1及小于0.91的数据，同时保证剩余缩放比数据的数量≥4，对剩余的缩放比数据求均值获得用于空间代谢组转换的最终缩放比；若剩余缩放比数据的数量＜4，则回到步骤(1)重新进行。

步骤(5)、通过空间转录组与空间代谢组中各标识点位两两连线的旋转角度判断空间转录组与空间代谢组样本检测过程中是否翻转并计算空间转录组与空间代谢组检测样本之间的最终旋转角度；

所述标志点位为经缩放比筛选后的标志点位；如果各标志点位两两之间连线所有旋转角度偏差在5°之内，说明样本未被翻转；如果各标志点位两两之间连线所有旋转角度偏差在5°之外，将步骤(3)中空间代谢组中y空间编号的转化公式修改为：

UMY＝-MetaY*resolution，

再进行空间转录组与空间代谢组检测样本中各标志点位两两之间的旋转角度计算，获得偏差在5°之内的所有旋转角度后对所有旋转角度求均值获得最终旋转角度；

所述旋转角度的计算公式如下：

Angle＝arctan((UMYa-UMYb)/(UMXa-UMXb))-arctan((UTYa-UTYb)/(UTXa-UTXb))；

其中，UTXa、UTYa、UTXb、UTYb分别为空间转录组中标志点a和b经统一空间坐标后的X，Y坐标值；UMXa、UMYa、UMXb、UMYb分别为空间代谢组中标志点a和b经统一空间坐标后的X，Y坐标值；arctan为反正切函数计算；Angle为空间转录组中标志点a和b的连线与空间代谢组中标志点a和b的连线间的旋转角度。

步骤(6)、通过最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，使转化后的空间代谢组坐标与空间转录组坐标匹配。

所述空间转录组中心标识点是指空间转录组检测芯片中心点位置；

所述转换空间标识方法是指通过空间坐标转换将空间代谢组中每个像素点的空间位置转换成空间转录组中的对应坐标，步骤(4)和步骤(5)中获得的最终缩放比和最终旋转角度同样满足下述计算公式：

Ratio’＝sqrt((UMX’-UTXo)

Angle’＝arctan((UMY-UMYo)/(UMX-UMXo))-arctan((UMY’-UTYo)/(UMX’-UTXo))；

其中，Ratio’为最终缩放比；Angle’为最终旋转角度；UTXo及UTYo为空间转录组检测芯片中心点位置统一坐标体系后的空间信息；UMXo及UMYo为空间转录组检测芯片中心点位置在空间代谢组对应位置的统一坐标体系后的空间信息；UMX、UMY为空间代谢组统一空间坐标后的X，Y坐标值；UMX’及UMY’为空间代谢组匹配空间转录组后的X，Y坐标值。

通过已获得的最终旋转角度、最终缩放比、空间转录组中心标识点坐标参数，即能够获取空间代谢组匹配空间转录组后的X，Y坐标值UMX’和UMY’。

由于空间转录组和空间代谢组分析体系只识别相应的空间信息模式，为了其空间信息能被相互识别，将统一并且关联的坐标体系还原转化为空间转录组和空间代谢组的空间信息体系；

统一坐标体系的空间信息转换为空间转录组空间信息体系的计算公式如下：

UMX’＝(TransX’/2+0.5)*transresolution；

UMY’＝TransY’*sqrt(0.75)*transresolution；

其中，UMX’、UMY’分别为空间代谢组匹配空间转录组后的空间信息；transresolution为空间转录组分辨率的数值；sqrt为开根号运算；TransX’、TransY’分别为空间代谢组匹配空间转录组空间信息后转化为空间转录组空间信息中的x和y的两组空间编号；

统一坐标体系的空间信息转换为空间代谢组空间信息体系的计算公式如下：

UMX’＝MetaX’*metaresolution；

UMY’＝MetaY’*metaresolution；

其中，UMX’、UMY’分别为空间代谢组匹配空间转录组后的空间信息；metaresolution为空间代谢组分辨率的数值；MetaX’、MetaY’为空间代谢组匹配空间转录组空间信息后转化为空间代谢组空间信息中x和y的两组空间编号。

本发明还提出了实现上述空间转录组与空间代谢组空间信息匹配方法的系统，及上述空间信息匹配方法在空间转录组与空间代谢组空间信息匹配中的应用。

所述空间转录组与空间代谢组空间信息匹配系统包括空间信息获取模块、空间坐标统一模块、空间参数运算模块、空间信息转化模块；

所述空间信息获取模块，将带有标识点位的包含空间转录组信息的HE染色图与空间代谢组成像图配准，获得标识点位在空间代谢组上的空间坐标位置信息；

所述空间坐标统一模块，将空间转录组和空间代谢组的空间信息转化为统一的空间信息标识；

所述空间参数运算模块，计算空间转录组与空间代谢组之间的最终缩放比和最终旋转角度；

所述空间信息转化模块，使用最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，使转化后的空间代谢组坐标与空间转录组坐标匹配。

点对点匹配方法中的数据关联是指通过对空间代谢组数据进行拟合，将拟合后数据进行空间转录组空间分布模式整合，最终实现空间代谢组与空间转录组在点与点层面上的数据关联，包括如下步骤：

步骤I、将空间代谢组的数据进行拟合，通过拟合提高空间代谢组的像素点分辨率；

所述拟合是指将空间代谢组中的像素点根据需要进行扩增，根据原始像素点位于样本区域和/或非样本区域，按照不同的拟合方法进行扩增。

若只有两个像素点，且均位于样本区域，两个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算两像素点之间扩增的像素点的底层数据：

MZ1ab＝(MZ1a+MZ1b)/2；

MZ2ab＝(MZ2a+MZ2b)/2；

…

MZnab＝(MZna+MZnb)/2；

其中，a与b像素点为相邻的像素点；MZ1a，MZ2a，…，MZna分别为a像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1b，MZ2b，…，MZnb分别为b像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1ab，MZ2ab，…，MZnab分别为像素点a与b之间扩增后的像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a和b之间扩增的像素点简写为{ab}。

若只有两个像素点，且一个像素点为样本区域像素点，另外一个像素点为非样本区域像素点，两个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算两像素点之间扩增的像素点的底层数据：

MZ1ab＝0；

MZ2ab＝0；

…

MZnab＝0；

其中，a与b像素点为相邻的像素点；MZ1ab，MZ2ab，…，MZnab分别为像素点a与b之间扩增后的像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a和b之间扩增的像素点简写为0，表示为非样本区域像素点。

若只有两个像素点，且两像素点均为非样本区域像素点，两个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算两像素点之间扩增的像素点的底层数据：

MZ1ab＝0；

MZ2ab＝0；

…

MZnab＝0；

其中，a与b像素点为相邻的像素点；MZ1ab，MZ2ab，…，MZnab分别为为像素点a与b之间扩增后的像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a和b之间扩增的像素点简写为0，表示扩增的像素点为非样本区域像素点。

步骤II、将空间代谢组拟合后的像素点与空间转录组的点spot进行关联；

所述关联是指将空间代谢组拟合后的像素点与空间转录组的点spot进行对应；通过计算空间代谢组拟合后的像素点中心位置与空间转录组的点spot中心位置之间的距离进行像素点的归属；若空间代谢组拟合后的像素点中心位置与空间转录组的某个点spot的中心位置之间的距离小于所述空间转录组的该点spot的直径，则空间代谢组拟合后的像素点对应于所述空间转录组中的点spot。

步骤III、将每个点spot对应的像素点底层数据进行加和运算；

所述加和运算是指将归属到同一个点spot上的多个像素点的底层数据整合成一个点spot的底层数据，通过如下公式计算：

MZ1spotA＝(MZ1a+MZ1b+…+MZ1x)/j；

MZ2spotA＝(MZ2a+MZ2b+…+MZ2x)/j；

…

MZnspotA＝(MZna+MZnb+…+MZnx)/j；

其中，MZ1a，MZ2a，…,MZna分别是Aspot关联的a像素点对应的MZ1，MZ2，…，MZn质核比的表达强度；MZ1b，MZ2b，…，MZnb分别是Aspot关联的b像素点对应的MZ1，MZ2，…，MZn质核比的表达强度；MZ1x，MZ2x，…，MZnx分别是Aspot关联的x像素点对应的MZ1，MZ2，…，MZn质核比的表达强度；j是Aspot关联的像素点数量；MZ1spotA，MZ2spotA，…，MZnspotA分别是对应Aspot的空间代谢组加和后的底层数据。

底层数据指的是在空间代谢组检测中在每个像素点内对应mz的表达强度数据。mz为带电粒子的质量数与电荷数之比，为物质特有数据信息，通常以mz指代物质。

步骤IV、获取得到空间转录组每个点spot对应的空间代谢组数据。

对空间代谢组数据以x和y的空间编号组合进行坐标位置处理，并在新空间坐标位置的每个mz对应下有对应表达强度。表达强度为空间代谢物在每个像素点对应mz的相对表达值，数值大小表示相对含量高低，值越大说明对应像素点的相对含量越高。

空间转录组数据以Barcode编码与上述空间代谢物坐标一一对应，并包含Barcode编码对应坐标的每个转录基因对应基因的count值。基因count值是在空间转录组中检测到对应转录基因的转录数量，转录表达是通过count值的大小判断转录组的表达高低。count值越大相应的转录组表达越高。

本发明还提出了实现上述空间转录组与空间代谢组数据关联方法的系统，及上述数据关联方法在空间转录组与空间代谢组数据关联中的应用。

所述空间转录组与空间代谢组数据关联系统包括空间代谢数据拟合模块、空间转录与空间代谢数据关联模块；

所述空间代谢数据拟合模块，将空间代谢组的数据进行拟合，并与空间转录组的点spot进行关联；

所述空间转录与空间代谢数据关联模块，将每个点spot对应的像素点底层数据进行加和运算，获得空间转录组每个点spot对应的空间代谢组数据。

本发明的有益效果包括：本发明通过一种空间转录组与空间代谢组点对点匹配方法及系统，实现了空间转录组与空间代谢组的点对点匹配，减少人工进行空间转录组与空间代谢组区域匹配的选择，点对点匹配后关联数据量大幅度提高，为后续空间转录组与空间代谢组关联分析提供更可靠的数据支持。

附图说明

图1表示本发明空间转录组与空间代谢组点对点匹配方法的流程示意图。

图2表示空间转录组芯片检测方法示意图。

图3表示空间代谢组质谱检测方法示意图。

图4表示空间转录组与空间代谢组空间信息匹配方法的流程示意图。

图5表示空间代谢组质谱检测方向及坐标标识方法示意图。

图6表示匹配方法实施例1中空间转录组中spot分布与HE染色示意图。

图7表示匹配方法实施例1中空间转录组选择spot点位标识示意图。

图8表示匹配方法实施例1中未标识点的空间代谢组成像图。

图9表示匹配方法实施例1中将空间转录组HE染色图与空间代谢组成像图匹配后将标识点标识后的示意图。

图10表示匹配方法实施例1中将空间代谢组空间转化为空间转录组空间信息的质谱成像图。

图11表示空间转录组与空间代谢组数据关联方法的流程示意图。

图12表示空间转录组检测空间分布方式示意图。

图13表示空间代谢组检测空间分布及数据存储方式示意图。

图14表示空间转录组与空间代谢组空间分布不同叠加效果图。

图15表示空间代谢组2*2像素点均为样本区域情况下，2次拟合后运算展示示意图。

图16表示空间代谢组2*2像素点中1个像素点为非样本区域情况下，2次拟合后运算展示示意图。

图17表示空间代谢组2*2像素点中2个像素点为非样本区域情况下，2次拟合后运算展示示意图。

图18表示空间代谢组2*2像素点中3个像素点为非样本区域情况下，2次拟合后运算展示示意图。

图19表示空间代谢组2*2像素点均为非样本区域情况下，2次拟合后运算展示示意图。

图20表示空间代谢组像素点与空间转录组spot之间的关联示意图。

图21表示2次拟合样本区域与非样本区域结果展示示意图。

图22表示对于样本区域及非样本区域的像素点与spot关联示意图。

图23表示实施例2空间代谢组数据转化成空间转录组分布后可视化图示意图。

图24表示实施例2空间代谢组数据拟合后成像图的变化示意图。

图25表示本发明实现上述空间转录组与空间代谢组点对点匹配方法的系统的示意图。

具体实施方式

结合以下具体实施例和附图，对本发明作进一步的详细说明。实施本发明的过程、条件、实验方法等，除以下专门提及的内容之外，均为本领域的普遍知识和公知常识，本发明没有特别限制内容。

本发明提出了一种空间转录组与空间代谢组点对点匹配方法及系统。

本发明提出的空间转录组与空间代谢组点对点匹配方法包括空间转录组和空间代谢组之间的空间信息匹配方法及数据关联方法。

其中，空间信息匹配方法用于空间转录组和空间代谢组中空间坐标位置信息的转化，使两套不同的组学体系的空间坐标位置能进行匹配。数据关联方法用于空间转录组和空间代谢组中底层数据处理，使两套不同的组学体系的数据进行关联。

具体地，通过空间转录组与空间代谢组空间信息匹配方法解决由于空间转录组与空间代谢组不同的空间信息标识、不同切片样本实验操作步骤及其非相同切片样本导致的空间位置不匹配的问题。

通过空间转录组与空间代谢组数据关联方法解决由于空间转录组上下两层成交错分布与空间代谢组正方形组合分布不同，导致无法将不同空间分布的数据进行空间位置上的关联问题。

将空间转录组与空间代谢组空间信息匹配方法和数据关联方法两个方法组合运行最终实现空间转录组与空间代谢组点对点的匹配。

具体地，所述空间信息匹配方法是指将空间转录组与空间代谢组的空间信息标识统一，并将空间代谢组的空间坐标体系匹配转换成空间转录组的空间坐标，实现空间转录组与空间代谢组的空间信息匹配。

如图4所示，包括如下步骤：

(1)在空间转录组的样本检测范围内选择6个及以上的点spot作为标识点位，获取标识spot对应的条形码Barcode空间信息并在对应HE(苏木精—伊红)染色图上标识spot。6个是最少的标识点位数量，少于6个标识点位时，匹配效果不好，且标识点位位置分散，分布在整个芯片检测区域。

(2)使用已经标识点位的HE染色图与空间代谢组成像图进行配准，获取配准后的标识点位在空间代谢组上的对应像素点的空间坐标位置信息。空间代谢组检测分辨率通常为100μm、40μm、20μm，即为像素点与像素点的距离为100μm、40μm、20μm。

(3)将空间转录组的barcode空间信息和空间代谢组空间信息转化为统一的空间信息标识。

(4)计算空间转录组不同标识点间距离；计算空间代谢组不同标识点间距离；并将空间转录组与空间代谢组对应的标识点间距离求缩放比。由于进行了空间转录组与空间代谢组的空间坐标体系的统一，理论缩放比为1，但由于空间代谢组与空间转录组的样本是相邻切片样本，在实际样本检测上是有不同的，为防止样本不同导致的选点的失误或者仪器误差导致的缩放比不为1，需要在后续运算中进行缩放比校正。

(5)剔除大于1.1及小于0.91的缩放比，剔除后剩余缩放比数量需大于等于4，(不符合要求需从步骤(1)重新开始)，计算其余缩放比均值获得最终缩放比，其作为空间转录组与空间代谢组最终缩放比，最终缩放比用于选点操作及实验操作误差导致的数据校正。

(6)使用符合缩放比条件的标识点位的空间信息判断空间转录组与空间代谢组样本检测过程中是否翻转及计算两个检测样本之间的最终旋转角度。所述最终旋转角度通过任意两个标志点位计算获得的所有旋转角度的均值，即计算在空间转录组中两点与在空间代谢组中对应两点的旋转角度，由于至少有4个以上标识点用于旋转角度的计算，最终会获得6种组合以上的旋转角度。在计算获得空间转录组和空间代谢组所有标识点组合后的旋转角度后，还需要判断旋转角度偏差是否≤5°，若是则进行后续的均值计算步骤，若否则回到步骤(3)重新开始。所述旋转角度偏差是指各个组合的旋转角度与最终旋转角度的绝对差值。

(7)通过最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，最终使转化过的空间代谢组坐标与空间转录组坐标匹配。

所述数据关联方法是指对空间代谢组数据进行拟合，将拟合后数据进行空间转录组空间分布模式整合，最终实现空间代谢组与空间转录组在点与点层面上的数据关联。

如图11所示，包括以下步骤：

(1)将空间代谢组的数据进行2次拟合，通过拟合后空间代谢组的像素点分辨率由原先的100μm×100μm提高到25μm×25μm；

(2)将空间代谢组拟合后的像素点与空间转录组的spot进行关联；

(3)将每个spot对应的像素点底层数据进行加和运算；

(4)获取得到空间转录组每个spot对应的空间代谢组数据。

实施例1

匹配方法

(1)在空间转录组的样本检测范围内选择6个spot作为标识点位，获取标识spot对应的barcode信息并在对应HE染色图上标识spot(见图6和图7)。图6为空间转录组中spot分布和HE染色的示意图；图7为空间转录组中选择的6个spot点位标识分布图。

表1空间转录组标识点的barcode及空间信息

表1中Barcode为在空间转录组技术中的空间条形码，每个spot中有其独有的barcode编码，能通过空间条形码获取对应spot的空间位置信息；TransX、TransY为空间转录组中x和y的两组空间编号；

(2)使用已经标识点位的HE染色图与空间代谢组成像图进行配准，获取配准后标识点位在空间代谢组上的对应像素点的空间坐标位置信息(见图8和图9)。图8为未标识点的空间代谢组成像图；图9为将空间转录组HE染色图和空间代谢组成像图匹配后包括标识点的示意图。

表2空间代谢组标识点的空间信息表

表2中MetaX、MetaY为空间代谢组中x和y的两组空间编号。

(3)将空间转录组的barcode空间信息和空间代谢组空间信息转化为统一的空间信息标识。

表3统一空间信息表

表3中UTX、UTY为空间转录组数据统一标准后空间坐标；UMX、UMY为空间代谢组数据统一标准后空间坐标。

(4)计算空间转录组不同标识点间距离；计算空间代谢组不同标识点间距离；并将空间转录组与空间代谢组对应的标识点间距离求缩放比。

表4两两标识点缩放比计算结果表

(5)剔除大于1.1及小于0.91的缩放比，剔除后剩余缩放比数量需大于等于4，计算其余缩放比均值获得最终缩放比，其作为空间转录组与空间代谢组最终缩放比，最终缩放比用于选点操作，选点操作是为了剔除引入过大误差的点，及实验操作误差导致的数据校正。最终缩放比为1.03。

(6)使用符合缩放比条件的标识点位的空间信息判断空间转录组与空间代谢组样本检测过程中是否翻转及计算两个检测样本之间的最终旋转角度。最终旋转角度为82.9°。各个组合的旋转角度与最终旋转角度的绝对差值称为旋转角度偏差。最终旋转角度为82.9°，各个组合的旋转角度与最终旋转角度均≤5°，均符合旋转角度偏差的阈值。

表5两两标识点旋转角度计算结果表

(7)通过最终旋转角度、最终缩放比、空间转录组中心标识点坐标三个参数使用转换空间标识方法将空间代谢组空间坐标转换，最终使转化过的空间代谢组坐标与空间转录组坐标匹配。如图10所示，图10为空间代谢组空间转化为空间转录组空间信息后的质谱成像图。

表6部分空间代谢组坐标转化为空间转录组空间信息表

表6中MetaX、MetaY为空间代谢组中x和y的两组空间编号；UMX、UMY为空间代谢组数据统一标准后空间坐标；UMX’、UMY’为空间代谢组匹配空间转录组后的空间信息；MetaX’、MetaY’为空间代谢组转匹配空间转录组空间信息后转化为空间代谢组空间信息中x和y的两组空间编号；TransX’、TransY’为空间代谢组转匹配空间转录组空间信息后转化为空间转录组空间信息中的x和y的两组空间编号。

实施例2

数据关联方法

数据关联方法中，通过以下例子对拟合提高空间代谢组的像素点分辨率进行说明：

如图15所示，若像素点按2*2排布，且均为样本区域像素点，四个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算2*2像素点中间扩增像素点的底层数据：

MZ1abcd＝(MZ1a+MZ1b+MZ1c+MZ1d)/4；

MZ2abcd＝(MZ2a+MZ2b+MZ2c+MZ2d)/4；

…

MZnabcd＝(MZna+MZnb+MZnc+MZnd)/4；

其中，像素点a、b、c、d按2*2排布；MZ1a，MZ2a，…，MZna分别为a像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1b，MZ2b，…，MZnb分别为b像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1c，MZ2c，…，MZnc分别为c像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1d，MZ2d，…，MZnd分别为d像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1abcd，MZ2abcd，…，MZnabcd分别为a、b、c、d四像素点之间扩增后的中间像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a、b、c、d之间扩增的像素点简写为{abcd}。

如图16所示，若像素点按2*2排布，且其中的1个像素点位于非样本区域，四个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算2*2像素点中间扩增像素点的底层数据，其中a像素点为非样本区域像素点：

MZ1abcd＝(MZ1b+MZ1c+MZ1d)/3；

MZ2abcd＝(MZ2b+MZ2c+MZ2d)/3；

…

MZnabcd＝(MZnb+MZnc+MZnd)/3；

其中，像素点a、b、c、d按2*2排布；MZ1b，MZ2b，…，MZnb分别为b像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1c，MZ2c，…，MZnc分别为c像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1d，MZ2d，…，MZnd分别为d像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1abcd，MZ2abcd，…，MZnabcd分别为a、b、c、d四个像素点之间扩增后的中间像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a、b、c、d之间扩增的像素点简写为{bcd}。

如图17所示，若像素点按2*2排布，且2个像素点位于非样本区域，四个像素点扩增后的中间像素点的计算方法如下：

按照下列公式计算2*2像素点中间扩增像素点的底层数据，其中a、b像素点为非样本区域像素点：

MZ1abcd＝(MZ1c+MZ1d)/2；

MZ2abcd＝(MZ2c+MZ2d)/2；

…

MZnabcd＝(MZnc+MZnd)/2；

其中，像素点a、b、c、d按2*2排布；MZ1c，MZ2c，…，MZnc分别为c像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1d，MZ2d，…，MZnd分别为d像素点对应MZ1，MZ2，…，MZn质核比的表达强度；MZ1abcd，MZ2abcd，…，MZnabcd分别为a、b、c、d四个像素点之间扩增后的中间像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a、b、c、d之间扩增的像素点简写为{cd}。

如图18所示，若像素点按2*2排布，其中3个像素点位于非样本区域，四个像素点扩增后中间像素点的计算方法如下：

按照下列公式计算2*2像素点中间扩增像素点的底层数据，其中a、b、c像素点为非样本区域像素点：

MZ1abcd＝0；

MZ2abcd＝0；

…

MZnabcd＝0；

其中，像素点a、b、c、d按2*2排布；MZ1abcd，MZ2abcd，…，MZnabcd分别为a、b、c、d四像素点之间扩增后的中间像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a、b、c、d之间扩增的像素点简写为0，表示扩增的像素点为非样本区域像素点。

如图19所示，若像素点按2*2排布，且均为非样本区域情况下，四个像素点扩增后中间像素点的计算方法如下：

按照下列公式计算2*2像素点中间扩增像素点的底层数据：

MZ1abcd＝0；

MZ2abcd＝0；

…

MZnabcd＝0；

其中，像素点a、b、c、d按2*2排布；MZ1abcd，MZ2abcd，…，MZnabcd分别为a、b、c、d四像素点之间扩增后的中间像素点对应MZ1，MZ2，…，MZn质核比的表达强度；像素点a、b、c、d之间扩增的像素点简写为0，表示扩增的像素点为非样本区域像素点。

实施例3

将小鼠肾脏空间代谢组数据与空间转录组关联

(1)将小鼠肾脏空间代谢组的数据进行2次拟合，通过拟合后空间代谢组的像素点分辨率由原先的100μm×100μm提高到25μm×25μm；(图24)

(2)将空代代谢组拟合后的像素点与空间转录组的spot进行关联；

表7部分空间代谢组坐标与空间转录组Barcode关联表

表7中Barcode为在空间转录组技术中的空间条形码，每个spot中有其独有的barcode编码，能通过空间条形码获取对应spot的空间位置信息；空间代谢组坐标为空间代谢组中x和y的两组空间编号的组合。

(3)将每个spot对应的像素点底层数据进行加和运算(图23)；底层数据指的是在空间代谢组检测中在每个像素点内对应mz的表达强度数据。mz为带电粒子的质量数与电荷数之比，为物质特有数据信息，通常以mz指代物质。

(4)获取得到空间转录组每个spot对应的空间代谢组数据。

表8部分空间代谢组数据

表8中mz为空间代谢组中质荷比数据；1-5、2-4、3-5、5-5、7-5为空间代谢组中x和y的空间编号组合对应下质荷比对应的表达强度。表达强度为空间代谢物在每个像素点对应mz的相对表达值，数值大小表示相对含量高低，值越大说明对应像素点的相对含量越高。

表9部分空间转录组数据

表9中GeneName列为空间转录组中基因名称；GCCACCCATTCCACTT、TACTCACAACGTAGTA、GTTCGGTGTGGATTTA、GGAGACATTCACGGGC、TAGAACGCCAGTAACG列的数值为空间转录组Barcode编码的对应基因的count值。

count值是在空间转录组中检测到对应GeneName的转录数量，转录表达定量通过count值的大小判断转录组的表达高低。count值越大相应的转录组表达越高。

本发明的保护内容不局限于以上实施例。在不背离本发明构思的精神和范围下，本领域技术人员能够想到的变化和优点都被包括在本发明中，并且以所附的权利要求书为保护范围。