CN116842996A - 一种基于深度压缩感知的空间转录组方法及装置 - Google Patents

一种基于深度压缩感知的空间转录组方法及装置 Download PDF

Info

Publication number
CN116842996A
CN116842996A CN202310514516.7A CN202310514516A CN116842996A CN 116842996 A CN116842996 A CN 116842996A CN 202310514516 A CN202310514516 A CN 202310514516A CN 116842996 A CN116842996 A CN 116842996A
Authority
CN
China
Prior art keywords
model
data
matrix
encoder
compressed sensing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202310514516.7A
Other languages
English (en)
Inventor
瞿昆
林俊
蒋辰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Science and Technology of China USTC
Original Assignee
University of Science and Technology of China USTC
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Science and Technology of China USTC filed Critical University of Science and Technology of China USTC
Priority to CN202310514516.7A priority Critical patent/CN116842996A/zh
Publication of CN116842996A publication Critical patent/CN116842996A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G06N3/0455—Auto-encoder networks; Encoder-decoder networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/0464—Convolutional networks [CNN, ConvNet]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Linguistics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明属于分子生物学和生物信息学领域，提供了一种自编码器，包含所述自编码器的压缩感知模型CSNet，基于深度学习的压缩感知方法，以及基于深度压缩感知的空间转录组方法及装置。基于本发明的基于深度压缩感知的空间转录组方法csFISH重构的数据不仅与smFISH技术相近，而且本发明成像所用时间仅仅是smFISH技术的1/28。并且，本发明的csFISH方法不仅可以只进行少量的采样即可恢复高维数据，还克服了光学拥挤的难题，能够用低倍镜更快速地成像。该方法在保证高质量数据的同时，极大提升了成像速度。

Description

一种基于深度压缩感知的空间转录组方法及装置
技术领域
本发明涉及一种空间转录组方法和装置，特别是涉及一种基于深度压缩感知的空间转录组方法及装置，属于分子生物学和生物信息学领域。
背景技术
空间转录组技术是一种用于研究组织中基因表达的新兴组学技术，能够捕获细胞层面的基因表达图谱并保留组织的位置信息。它是对特定性质的细胞类型、组织或器官内的转录本表达进行全面和系统研究的方法，并且不会失去细胞在组织中的空间位置信息。相比于传统的转录组技术，空间转录组技术可以获得细胞在组织生理环境下真实的基因表达特征，并且此技术能够探明细胞与微环境之间的关系，有助于加深对正常和病理状态下细胞特性的理解。因此，它成为了分子生物学和生物医学领域重要的研究方法之一。
目前，主流的空间转录组技术主要分为两类：基于下一代测序的空间转录组技术和基于原位成像的空间转录组技术。其中，原位成像技术又可分为基于原位杂交的空间转录组技术和基于原位测序的空间转录组技术。
基于下一代测序的空间转录组技术利用芯片给组织切片中的RNA添加空间位置信息，之后再进行测序并重构其空间分布。这一类技术目前已有多个商业化产品，包括美国10X公司的Visium技术、Slide-seq技术和华大的Stereo-seq技术等。
虽然，Visium技术、Slide-seq技术和Stereo-seq技术具有各自的优点，但是，基于下一代测序的空间转录组技术不能同时具有高空间分辨率与高测序深度的特点，难以获取高质量的单细胞精度的数据。同时，其成本随样品数目线性增加，不利于大规模样品的实验。
主流的基于原位成像的空间转录组技术通过结合荧光原位杂交(FISH)技术和显微成像技术，能够实现对单个细胞中不同mRNA转录本的定位和统计，从而得到单细胞水平的空间转录组信息。其中，FISH技术是利用核酸探针对目标mRNA转录本进行特异性的标记并搭配荧光物质进行检测；而显微成像技术则是通过高分辨率的成像设备对这些荧光信号进行捕获；然后使用图像分析算法进行图像处理，得到每个细胞中mRNA的准确位置，即“亚细胞精度”信息；最后使用细胞分割算法识别和统计每一个细胞中的mRNA数量，从而获得空间转录组信息。
smFISH技术、seqFISH技术和MERFISH技术是目前使用较为普遍的基于原位成像的空间转录组技术。但是，在基于原位成像的空间转录组技术中，基因通量与成像速度不能兼得。smFISH技术基因通量低；seqFISH技术采用编码方案，能够测量的基因数最多，但是由于受到光学衍射极限的限制，存在光学拥挤的问题；而MERFISH技术为了避免光学拥挤的问题，降低了编码压缩率，采用了更高标准的光学仪器来获得高分辨率图像，所需成像时间较长。
因此，本领域亟需一种新的空间转录组技术，其既能兼顾基因通量和成像速度，也能避免光学拥挤的问题。
发明内容
本发明针对现有技术中存在的上述缺陷，提供一种自编码器，其是将输出恢复成输入的神经网络结构，包含编码器和解码器，其中：
编码器是由一个线性压缩层组成的神经网络，将高维数据转换为低维表示，实现数据的压缩；
解码器由多个全连接层组成，其采用残差网络ResNet作为网络结构，将压缩后的数据重构成原始输入数据，实现数据的重构。
在本发明优选的实施方案中，所述自编码器通过最小化输入x和输出之间的误差同时完成压缩和重构两个步骤，其中，
编码器的线性压缩层可以将基因表达矩阵X通过与采样矩阵A进行矩阵乘法进行线性变换得到复合测量矩阵Y，其中A是m×g的二维矩阵，X是g×n的二维矩阵，Y是m×n的二维矩阵，编码器的输出可以表示为如下公式，
Y＝AX
通过所述公式，编码器将高维的输入数据转换为低维的线性表示，在数据采集时实现压缩；
解码器将编码器的输出结果Y作为解码器的输入，通过如下公式，重构基因表达矩阵
。
本发明另一方面提供了一种基于深度学习的压缩感知模型CSNet，其包含本发明所述的自编码器。
本发明另一方面提供了一种基于深度学习的压缩感知方法，其将深度学习技术应用到压缩感知中，包括以下步骤：
1、构建本发明所述的自编码器，或者本发明所述的模型CSNet；
2、采用步骤1构建的自编码器或模型CSNet进行模型训练和推理，在进行模型训练中，通过反向传播算法更新网络权重以最小化损失函数。
在本发明优选的实施方案中，步骤2包含编码器和解码器联合训练、解码器单独训练和解码器推理三个过程。
在本发明优选的实施方案中，步骤2中的模型训练包括训练步骤、验证步骤以及测试步骤。
在本发明优选的实施方案中，使用基于PyTorch开发的Pytorch-Lightning高级框架来实现该压缩感知算法。
在本发明优选的实施方案中，步骤2在进行模型训练之前，还包括对数据进行处理的步骤，更加优选对复合测量矩阵进行归一化和对数转换。
本发明另一方面提供了本发明所述的自编码器或者本发明所述的模型CSNet在基于深度学习的压缩感知方法中的应用。
本发明另一方面提供了一种基于深度压缩感知的空间转录组方法，包括以下步骤：
1、模型构建，构建本发明所述的自编码器，或者本发明所述的模型CSNet；
2、模型训练，获得训练好的模型和采样矩阵；
3、组织成像，进行荧光原位杂交成像实验得到数据图像；
4、数据重构，从图像中提取复合测量矩阵并输入压缩感知模型中重构基因表达矩阵。
在本发明优选的实施方案中，步骤2通过单细胞数据集训练模型，得到采样方案。
在本发明优选的实施方案中，步骤2还包括根据皮尔逊相关系数指标和聚类指标ARI选择最优的模型参数和采样矩阵。
在本发明优选的实施方案中，步骤3还包括对图像进行滤波去噪以及图像处理步骤，所述图像处理优选为对图像进行归一化、拼接、对齐和分割。
在本发明优选的实施方案中，步骤4还包括对对经过图像处理得到的复合测量矩阵进行预处理的步骤，所述预处理优选为进行归一化和对数转换。
本发明另一方面提供了本发明所述的自编码器或者本发明所述的模型在基于深度压缩感知的空间转录方法中的应用。
本发明再一方面提供了一种基于深度压缩感知的空间转录组装置，包括以下模块：
1、模型构建模块，用于构建本发明所述的自编码器，或者本发明所述的模型CSNet；
2、模型训练模块，用于获得训练好的模型和采样矩阵；
3、组织成像模块，用于进行荧光原位杂交成像实验得到数据图像；
4、数据重构模块，用于从图像中提取复合测量矩阵并输入压缩感知模型中重构基因表达矩阵。
在本发明优选的实施方案中，模型训练模块通过单细胞数据集训练模型，得到采样方案。
在本发明优选的实施方案中，模型训练模块根据皮尔逊相关系数指标和聚类指标ARI选择最优的模型参数和采样矩阵。
在本发明优选的实施方案中，组织成像模块还包括滤波去噪和图像处理模块，滤波去噪用于对图像进行滤波将噪声去除，图像处理用于对图像进行处理，所述图像处理更加优选为对图像进行归一化、拼接、对齐和分割。
在本发明优选的实施方案中，数据重构模块还包括复合测量矩阵预处理模块，用于对经过图像处理得到的复合测量矩阵进行预处理，所述预处理更加优选为进行归一化和对数转换。
本发明通过结合深度学习方法和压缩感知方法开发了一种重构空间转录组的新算法，并验证了在同等条件下该算法能够比其他工具更好地重构空间转录组。然后，通过结合深度压缩感知算法和荧光原位杂交成像技术，本发明开发出一种新的空间转录组技术，该技术可以在低倍镜成像条件下对基因探针进行少量混合采样，将采样的复合测量矩阵输入到已经从单细胞数据集中学习到先验知识的模型中重构出高维的空间转录组数据，在提升了成像速度的同时提高了基因通量。本发明不仅在保证高基因通量的同时加快了空间转录组技术的成像速度，还解决由于mRNA转录本过多和图像分辨率不足而导致的光学拥挤问题。
本发明通过深度压缩感知算法将低倍镜和欠采样条件下获取的数据重构成高质量的基因表达矩阵，能够大大减少图像采集时间和图像文件大小，从而加快了数据的采样、传输和处理速度。训练的模型压缩率越高，成像所用的镜头倍数越低，则成像速度越快。与smFI SH技术相比，本发明的csFI SH技术有望在10倍镜、5倍镜的成像条件下通过几十轮测量成像数千个基因，能够将成像速度提升上千倍。本发明的应用将有助于大规模空间转录组测序，从而更全面地解析大组织切片中的基因表达分布，为疾病诊断和机制研究提供了有力的工具。
综上，基于本发明的基于深度压缩感知的空间转录方法csFISH重构的数据不仅与smFISH技术相近，而且本发明成像所用时间仅仅是smFISH技术的1/28。并且，本发明的csFISH方法不仅可以只进行少量的采样即可恢复高维数据，还克服了光学拥挤的难题，能够用低倍镜更快速地成像。该方法在保证高质量数据的同时，极大提升了成像速度。
附图说明
图1为本发明基于深度学习的压缩感知方法的整体流程图；
图2为本发明基于深度压缩感知的空间转录组方法的整体流程图；
图3为本发明基于深度压缩感知的空间转录组装置的整体结构图；
图4为本发明基于深度压缩感知的空间转录组装置的优选结构图；
图5为本发明自编码器网络结构图；
图6为本发明残差网络结构图；
图7为本发明5种算法在10个空间数据集中的评估结果图；
图8为本发明5种算法的性能随细胞数变化的结果图；
A：箱型图根据细胞占比分组展示；B：箱型图根据模型分组展示；
图9为本发明5种算法的性能随测量数变化的结果图；
A：箱型图根据multiple分组展示；B：箱型图根据模型分组展示；
图10为本发明csFISH方法流程图；
图11为本发明5种算法工具重构的113个基因的聚类效果图；
图12为本发明去噪前后图像对比图；
A：去噪前图像；B：去噪后图像；
图13为本发明对细胞核进行DAPI染色和细胞掩膜对比图；
A：细胞核DAPI染色通道；B：细胞掩膜；
图14为本发明基于csFISH方法的小鼠海马体空间转录组聚类结果图；
A：csFISH数据的聚类结果UMAP分布；B：8个标记基因Gad1，Slc17a7，Slc17a6，Ctgf，Cldn11，Gja1，Cd14，Esam的UMAP分布；
图15为本发明csFISH方法与smFISH技术的小鼠海马体空间转录组图谱；
A：基于csFISH方法测得的数据的聚类结果；
B：基于csFISH方法测得的数据的8个标记基因(Gad1，Slc17a7，Slc17a6，Ctgf，Cldn11，Gja1，Cd14，Esam)的空间分布；
C：基于smFISH技术测得的数据的聚类结果；
D：基于smFISH技术测得的数据的8个标记基因(Gad1，Slc17a7，Slc17a6，Ctgf，Cldn11，Gja1，Cd14，Esam)的空间分布；
E：基于csFISH与smFISH技术获得的平均表达量间的相关性，其中点代表基因，皮尔逊相关性为0.875；
F：基于csFISH与smFISH技术测得的数据的10个亚群的平均表达量间的相关性热图。
具体实施方式
下面通过附图和实施例对本申请进一步详细说明。通过这些说明，本申请的特点和优点将变得更为清楚明确。
实施例1：转录组数据收集与分析方法
(一)数据收集与处理
从GEO(Barrett T,Wilhite S E,Ledoux P,et al.NCBI GEO:archive forfunctional genomics data sets—update[J].Nucleic acids research,2012,41(D1):D991-D995)和NEMO(Iannuccelli E,Mompart F,Gellin J,et al.NEMO:a tool foranalyzing gene and chromosome territory distributions from 3D-FISHexperiments[J].Bioinformatics,2010,26(5):696-697.)数据库中获取单细胞转录组和空间转录组数据集。
使用scanpy(Wolf F A,Angerer P,Theis F J.SCANPY:large-scale single-cell gene expression data analysis[J].Genome biology,2018,19:1-5.)，这是一个为单细胞转录组数据处理和分析设计的Python(Python W.Python[J].Python ReleasesWind,2021,24.)库。为了方便操作，将所有数据集处理为h5ad格式，并进行预处理，去除低质量的细胞。
(二)数据预处理
整理好数据集后，首先需要对数据进行质量控制，使用sc.pp.filter_cells函数过滤低质量的细胞和使用sc.pp.filter_genes函数过滤低质量的基因，再去除掉线粒体基因表达占比过高的细胞。单细胞转录组数据有一个标准的文库大小(library size)，即总基因表达量，位于标准的预期范围之外的细胞(过低或过高)可能代表无需进行下游分析的低质量细胞，将这些细胞过滤掉。由于不同细胞中所捕获的mRNA转录本数量差异较大，使用sc.pp.normalize_total函数对每个细胞的基因表达量进行标准化处理，以消除这些差异对后续分析的影响。并且由于基因表达服从幂律分布(skewed distribution)，需要使用sc.pp.log1p函数对特征矩阵进行log转换，转换后的结果更符合高斯分布(normaldistribution)，以使单细胞数据更加可靠地进行下游分析和可视化。
(三)聚类分析
细胞聚类是一种常用的数据分析技术，可以将具有相似基因表达模式的细胞归为同一类别，以便更好地理解细胞功能和细胞分化，并揭示不同疾病中细胞的变化模式。在处理高维数据时，需要先将其映射成低维稠密矩阵才能进行聚类分析。主成分分析(PCA)(Ratajczak W.Principal components analysis(PCA)[J].Computers&Geosciences,1993,19(3):303-342)是一种通用的降维方法，在单细胞转录组数据处理中可以使用sc.tl.pca函数来实现降维。接下来，可以在PCA结果的基础上使用sc.pp.neighbors函数构建K近邻图(KNN)(Eppstein D,Paterson M S,Yao F F.Onnearest-neighbor graphs[J].Discrete&Computational Geometry,1997,17:263-282.)，并用sc.tl.leiden函数调用聚类算法对数据进行聚类(Traag V A,Waltman L,Van Eck NJ.From Louvain to Leiden:guaranteeing well-connected communities[J].Scientific reports,2019,9(1):5233.)。
(四)数据可视化
使用PCA方法得到的低维数据仍然具有较多的维度，不易于可视化展示。采用流形学习方法(如：t-SNE和UMAP)(Van Der Maaten L,Hinton G.Visualizing data using t-SNE[J].Journal of machine learning research,2008,9(11).；Mcinnes L,Healy J,Melville J.Umap:Uniform manifold approximation and projection for dimensionreduction[J].arXiv preprint arXiv:1802.03426,2018.)能够将高维数据投影到二维空间并保留数据间的局部线性关系，是一种常用的降维方法。在单细胞转录组数据处理中，可以使用sc.tl.umap函数进行UMAP降维，并使用sc.pl.umap函数进行可视化展示。这种可视化方式能够更好地展示单细胞转录组数据之间的关系，有助于研究者观察基因表达水平和聚类信息来进一步探索数据的生物学意义。
(五)基因差异分析
基因差异分析是基于单细胞转录组数据，用来识别在不同条件下具有显著差异的基因。该分析可以帮助研究者更好地理解生物学过程和疾病机制，同时也有助于筛选出与特定生物学过程或疾病相关的候选基因。sc.tl.rank_genes_groups函数使用Wilcoxon秩和检验方法来计算每个基因在两个不同群体之间的表达差异，并生成一个基因表达差异评分(Liao C,Li S,Luo Z.Gene selection using wilcoxon rank sum test and supportvector machine for cancer classification[J].Lecture Notes in ComputerScience,2007,4456:57.)。分数越高表示该基因在不同群体之间的表达差异越大。排名靠前的基因则被认为与不同的条件有关。通过这种基因差异分析方法，研究者可以更深入地了解单细胞转录组数据中的基因表达模式，并从中挖掘出有意义的生物信息。
实施例2：深度压缩感知模型设计
压缩感知的基本思想是，对于一个具有稀疏性的实际信号，只需要进行较少的随机采样就可以在不失真的情况下恢复出完整的信号。这是因为信号在某些变换域内是稀疏的。因此，在压缩感知中，采样过程可以看作是在信号的稀疏表示下获得一组线性采样，也被称为“欠采样”。采用压缩感知还可以避免信号受到采样噪声的影响。
基于深度学习的压缩感知算法将深度学习技术应用到压缩感知领域中，进一步提高了数据的压缩效率和信号重建质量。该算法通常包括压缩和重构两个部分。压缩部分将原始信号转换为低维表示，而重构部分使用该低维表示来重建原始信号。相比传统的压缩感知方法，基于深度学习的压缩感知算法在信号重建方面表现更好。
本发明使用基于PyTorch开发的Pytorch-Lightning(Falcon W A.Pytorchlightning[J].GitHub,2019,3.)高级框架来实现该压缩感知算法。相比于直接使用PyTorch，Pytorch-Lightning提供了更高级的功能和框架。它增加了代码复用性，简化了代码结构，提高了训练效率，并支持自动化日志记录。这使得编写模型和训练代码更容易。
(一)模型框架
本发明开发了一个新的用于重构空间转录组的算法工具CSNet，该算法采用了自编码器(Autoencoder)(Hinton G E,Salakhutdinov R R.Reducing the dimensionalityof data with neural networks[J].science,,2006,313(5786):504-507.；Kingma D P,Welling M.Auto-encoding variational bayes[J].arXiv preprint arXiv:1312.6114,2013.)作为基本的网络框架。自编码器是一种常用的生成模型，它由编码器和解码器两个主要部分组成(如附图5所示)。在自编码器中，编码器将高维数据转换为低维表示，实现了数据的压缩。而解码器则将压缩后的数据重构成原始输入数据，实现数据的重构。这种方法不同于传统方法，它没有将压缩和重构的步骤分开进行，自编码器通过最小化输入x和输出之间的误差同时训练两个步骤。自编码器可以从数据中自动学习到其特征，具有更强的适应性和泛化性能。最小化输入x和输出/>之间的误差可表示为公式(1)：
编码器是由一个线性压缩层组成的神经网络。线性压缩层可以将基因表达矩阵X通过与采样矩阵A进行矩阵乘法进行线性变换得到复合测量矩阵Y。其中A是m×g的二维矩阵，X是g×n的二维矩阵，Y是m×n的二维矩阵。为了在实验中能测量基因的线性组合，采样矩阵A需服从伯努利分布，其中Aij＝1表示第i个测量中测量第j个基因，Aij＝0表示第i个测量中不测量第j个基因。这样就能够将高维的输入数据转换为低维的线性表示，在数据采集时实现压缩，将通过Gumbel-softmax激活函数训练模型得到采样矩阵A。编码器的输出可以表示为公式(2)：
Y＝AX (2)
解码器由多个全连接层组成，编码器采用一个函数g，，将编码器的输出结果Y作为解码器的输入，连接多个网络层，并在最后添加非线性激活函数，能够拟合复杂函数g，重构基因表达矩阵解码器的输出可以表示为公式(3)：
目前神经网络优化方法都是基于反向传播算法实现的，即损失函数计算的误差通过梯度反向传播的方式进行传递，更新网络权重。其中，将误差从最后一层向前传递的过程通过链式法则(ChainRule)来实现，，链式法则是一个连乘的形式，梯度将以指数形式传播。随着网络层数的增加会变得越来越明显，得到的梯度值接近0或特别大更新网络权重时，就会出现梯度消失或爆炸，从而导致难以训练深层神经网络。
残差网络(ResNet)(He K,Zhang X,Ren S,et al.Deep residual learning forimage recognition[C].Proceedings of the IEEE conference on computer visionand pattern recognition,2016:770-778.)的结构(如附图6所示)引入了跨层连接，一个残差模块直接将该模块的输入x与该模块的输出F(x)相加，得到下一个后续残差模块的输入H(x)，这种跨层连接能够保留原始信号信息，缓解了梯度消失和梯度爆炸的问题，并使得网络更易于训练和优化。每个残差块的输出可以表示为公式(4)：
H(x)＝F(x)+x (4)
为了使得模型更容易训练，需要对数据进行一些处理，包括：复合测量矩阵归一化和对数变换，批次标准化(Batch Normalization)(Ioffe S,Szegedy C.Batchnormalization:Accelerating deep network training by reducing internalcovariate shift[C].International conference on machine learning,2015:448-456.)。由于实验误差导致mRNA转录本的捕获率和测序深度有差异，需要对复合测量矩阵Y进行归一化操作，将细胞n的基因集g表达量Yng除以文库大小ln得到基因表达占比Png。且不同基因的平均表达量差异较大，服从负二项分布(Di Y,Schafer D W,Cumbie J S,etal.The NBP negative binomial model for assessing differential gene expressionfrom RNA-Seq[J].Statistical applications in genetics and molecular biology,2011,10(1).)，需要对归一化后的Png进行对数变换得到Zng，使其服从高斯分布，更有利于模型训练。归一化和对数变换可以分别表示为公式(5)和(6)：
此外，还在每一个残差块中添加了批次标准化层，批次标准化对每个数据批次batch进行标准化处理，它将每个输入特征按其均值和方差进行归一化，然后应用比例系数和偏移参数进行放缩和平移操作。经过批次标准化处理的输入具有相似的统计学特性，使得模型训练更稳定。
由于在将复合测量值输入到解码器前需要对数据进行归一化，故需要知道该数据集的文库大小，即需要测量一个全基因集表达量。通过增加一个全基因集测量值ln，不仅能得到m个基因集的表达占比，还能通过将文库大小ln减去一个复合测量值Yng得到一个新的复合测量值Yng′。因此通过m+1次测量，可以得到2m+1个基因集的表达量。
(二)模型训练与推理
在深度学习领域中，模型训练与推理是两个非常重要的过程。模型训练是指在已知数据集上，通过调整模型参数来最小化训练数据和真实结果之间的误差的过程。模型推理是指在模型训练完成后，使用训练好的模型对新的和未知的数据进行预测或分类的过程。本研究中模型训练与推理分为三个过程：编码器和解码器联合训练过程，解码器单独训练过程和解码器推理过程。为了能够客观评估和改进模型的性能，确保模型能够在实际应用中表现良好，每个模型训练过程又包括训练步骤、验证步骤以及测试步骤。
在训练步骤中，模型将根据反向传播算法计算(Rojas R,Rojas R.Thebackpropagation algorithm[J].Neural networks:a systematic introduction,1996:149-182.)的梯度更新每个节点的参数，以最小化损失函数，更好地拟合训练集数据。为了避免模型过拟合，需要对训练集进行随机打乱，同时还需要向解码器的输入数据添加随机噪声以模拟真实的实验数据。但是仅在训练集表现良好并不能保证模型能够很好地泛化到未知数据当中。因此，需要在每一轮训练结束后使用未知的数据集对模型性能进行验证。
在验证步骤中，使用验证集数据评估模型在未知的数据上的性能，无需进行反向传播来更新节点参数，也不需要添加噪声和打乱数据集。程序根据每轮验证集的评估结果选择最佳的模型参数并自动保存，这个过程称为超参数调优。验证步骤的另一个重要作用是提供早期停止策略(Ying X.An overview of overfitting and its solutions[C].Journal of physics:Conference series,2019:022022.)，当模型在多轮训练后性能不再改善时，可以及时停止训练以避免过拟合。本发明选择了5轮作为早期停止的参数。
在模型训练到最大轮次或执行早期停止策略后，程序将进入测试步骤。测试步骤与验证步骤类似，但所用的数据集相互独立。由于模型是根据在验证集上的评估结果进行挑选的，为了更加客观地评估模型性能，需要使用独立于训练集和验证集的数据集对模型进行测试。该步骤的评估结果将作为模型的最终评估结果保存。
编码器和解码器联合训练是指将编码器和解码器结合起来进行训练。在联合训练中，编码器将高维数据转换成低维数据表示，解码器将低维数据重构成原始输入数据。采样矩阵是m×g矩阵，表示有m个混合采样，每个混合采样向量由g个0或1两种元素组成，表示对哪些基因进行采样。后续将导出编码器的网络参数作为采样矩阵，因此编码层参数需要符合伯努利分布。伯努利分布是离散分布，无法求导，为了能够对符合离散分布的参数进行梯度计算使用Gumbel-Softmax(Jang E,Gu S,Poole B.Categorical reparameterizationwith gumbel-softmax[J].arXiv preprint arXiv:1611.01144,2016.)作为激活函数，它可以对离散分布进行重参数采样使其变得可导。
该方法的基本思想是基于Gumbel分布(Nadarajah S,Kotz S.The beta Gumbeldistribution[J].Mathematical Problems in engineering,2004,2004(4):323-332.)进行重参数采样在概率空间中引入噪声使之变得可训练，并使用softmax对argmax进行平滑近似使其可导，可以得到近似于极限分布的结果，从而得到符合伯努利分布的采样参数。使用Gumbel-Softmax作为激活函数，可以通过数据驱动方式学习到采样矩阵，提高模型的性能。
由于Gumbel-Softmax中引入了噪声，其输出结果具有一定的随机性。为了避免这种随机性影响到解码器的训练，并保证训练结果的可靠性，需要固定采样矩阵和冻结编码器的网络参数，并对解码器进行单独训练。在解码器单独训练过程中，仅训练模型的重构能力通过最小化解码器的重构误差来提高模型的性能。完成两个训练过程后，可以使用已经训练好的模型对新得到的数据进行推理。
具体来说，将采样得到的复合测量矩阵输入解码器即可重构基因表达矩阵。在这个过程中，只需要使用已经训练好的模型，无需再进行任何的训练或参数调整。因此，模型的推理步骤能够高效且快速地重构大规模基因表达矩阵。
(三)损失函数
常见的损失函数有零膨胀的负二项分布(ZINB)(Minami M,Lennert-Cody C E,Gao W,et al.Modeling shark bycatch:the zero-inflated negative binomialregression model with smoothing[J].Fisheries Research,2007,84(2):210-221.)和均方误差(MSE)。在实际使用中，ZINB可能会导致一些假阴性，而MSE则容易受到极大值的影响，因此本发明选择了一种新的损失函数对数损失函数(logMSE)，该函数计算的是对数转换后的预测值与原始值直接的MSE。LogMSE的计算方法如公式(7)：
此外，为了增加模型的鲁棒性和泛化能力，在损失函数中增加余弦损失约束项。该约束项是指在解码器单独训练过程中，针对同一批次数据的预测结果，计算该结果与真实标签之间的余弦相似度，并将其作为损失函数中的一部分。这个损失函数的目的是使得同一批次数据的预测结果之间相互靠近，从而减少模型的随机性，提高模型的稳定性。同时，余弦损失约束还可以帮助模型更好地学习到数据之间的相似性，从而提高模型的泛化能力。本发明一共添加了三种余弦损失约束项，分别为细胞间的余弦损失项，基因间的余弦损失项，全局余弦损失项。用predict表示n×g的二维重构矩阵，target表示n×g的二维原始矩阵，其中细胞间的余弦损失项和基因间的余弦损失项可分别表示为公式(8)和(9)：
分别将二维矩阵predict和target展平为b*g个元素的一维向量predic’和target’，其全局余弦损失项可分别表示为公式(10)：
最终，本发明构造的损失函数如公式(11)：
Loss＝LogMSE+Icos1+Lcos2+Lcos3 (11)
实施例3：空间转录组重构算法性能比较分析
(一)空间转录组与单细胞转录组数据收集
本发明分别收集了已发表的27个单细胞转录组数据集作为训练和验证数据集以及10个基于成像技术的空间转录组作为金标准测试数据集(参见表1和表2)(Li B,ZhangW,Guo C,et al.Benchmarking spatial and single-cell transcriptomicsintegration methods for transcript distribution prediction and cell typedeconvolution[J].Nat Methods,2022,19(6):662-670.)。所有数据集的物种来源，细胞数目，基因数目等注释信息如表3和表4所示。
表1单细胞转录组数据集信息
表2基于成像的空间转录组数据集信息
表3本发明中27个单细胞转录组数据集信息
表4本发明中10个基于成像的空间转录组数据集信息
(二)空间转录组与单细胞转录组数据预处理
本发明按照以下步骤对每个数据集进行预处理：
1、去除低质量细胞和基因
本发明使用scanpy工具去除scRNA-seq数据中捕获RNA种类少于200的细胞，少于50个细胞中表达的基因，捕获mRNA转录本数量少于2500的细胞以及线粒体基因表达占比大于5％的细胞。
2、标准化单细胞基因表达矩阵
本发明构建的基于深度学习的压缩感知算法接收的输入数据是原始表达矩阵，因此需要保留单细胞基因表达矩阵与空间转录组数据集的原始表达量。但是为了挑选高度可变的基因，需要对单细胞基因表达矩阵进行标准化处理。本文使用sc.pp.normalize_total函数归一化数据，使得每个细胞的文库大小相同，然后使用sc.pp.log1p函数进行对数转换。
3、挑选高度可变的基因
本发明使用了sc.pp.highly_variable_genes函数计算了每个基因的变异系数，并保留前20000个高变异基因及它们高可变性的排序结果，其计算方法如公式(12)：
其中CVi表示基因i的变异系数；δi表示单细胞中基因i表达分布的标准差；ui表示基因i在所有细胞中的平均表达水平。
(三)数据集构建
本发明比较了5种工具在使用相同测量数与同样的数据集的条件下重构空间转录组的性能。其中三种为基于单基因测序的插值算法，分别是：gimVI，SpaGE，Tangram。另外两种为基于复合测量的压缩感知算法，分别为：SMAF，CSNet。CSNet是本发明提出的算法工具。由于基于单基因测量的插值算法和基于复合测量的压缩感知算法的采样方法不同，需根据算法特性设计采样方案。为了使得单基因测量得到的信息最多，将根据highly_variable_rank函数选择最具有高可变性的基因。复合测量的采样矩阵则由CSNet和SMAF内置的函数方法获得，CSNet由于预设了一个全基因集的复合测量值，训练得到的采样矩阵A的测量数将比SMAF少1个。
本发明一共收集了10个基于成像的空间转录组数据集和27个单细胞数据集，并构造了包括训练集、验证集和测试集在内的10组数据集。为了保证算法评估的客观性，每组数据集中使用空间转录组数据作为测试集，并从27个单细胞数据集中选取与空间转录组数据集重叠的基因数最多的5个数据集作为训练集和验证集。这样的构建方案使得训练集和验证集的数据分布相似，模型能够正常训练，使用分布外的空间转录组数据集作测试集则符合实际应用场景。
每组数据集的细胞数和基因数都不同，为了保证公平性，每组实验将根据该组数据的基因数设置与MERFISH汉明编码方案相同的测量数(如表5所示)。为了观察模型在不同测量数下的表现，设一组数据集的汉明编码方案的测量数为m，在测量数分别为m/3、m/2、m、2m和3m的条件下比较模型的重构性能，总计50个任务。为了观察细胞数对模型性能的影响，将通过下采样细胞数量的方法在细胞比例分别为1％、5％、20％和100％条件下比较模型的重构性能，总计40个任务。
此外基于复合测量的压缩感知算法的性能对采样矩阵的稠密度较为敏感，因此需要为不同数据集搜索合适的稠密度参数，搜索范围设置为：1％，2％，5％，10％，20％和50％。然后根据验证集的表现情况选择最佳的模型。所有模型的任务总计1350个。
表5汉明编码方案中基因数与测量数对应关系
(四)算法工具的性能评估指标
X是n×g的二维原始表达矩阵，是n×g的二维重构表达矩阵，其中n代表细胞数，g代表基因数。X′是由二维矩阵X展开的n×g个元素的一维原始表达向量，/>是由二维矩阵展开的n×g个元素的一维重构表达向量。
在本发明中，使用以下四个相关性指标来衡量每种重构算法的预测准确度和调整后的兰德系数评估重构数据的聚类效果。
1、细胞间的皮尔逊相关系数(cellPCC)
其计算方法如公式(13)：
cellPCC是所有细胞的原始和重构表达向量间的相关性的平均值。其中Xi,j和分别表示细胞i中基因j的原始和重构表达向量；μ,i和/>分别表示细胞i在原始和重构表达向量中所有基因的平均表达水平；σi和/>分别表示细胞i原始和重构表达向量中所有基因的标准差。最终求得所有细胞的原始和重构表达向量间的相关性，再对其求平均即可计算出cellPCC。
2、基因间的皮尔逊相关系数(gene PCC)
其计算方法如公式(14)：
genePCC是所有基因的原始和重构表达向量间的相关性的平均值。其中Xi,j和分别表示细胞i中基因j的原始和重构表达向量；μ,j和/>分别表示基因j在原始和重构表达向量中所有细胞的平均表达水平；/>分别表示基因j原始和重构表达向量中所有细胞的标准差。最终求得所有细胞的原始和重构表达向量间的相关性，再对其求平均即可计算出genePCC。
3、基因表达矩阵的全局皮尔逊相关系数(overallPCC)
其计算方法如公式(15)：
overallPCC是一维原始和重构表达量间的相关性。其中X'i和分别表示原始和重构表达向量中的单个元素i的表达量；μ和/>分别表示原始和重构表达向量中所有表达量的平均表达水平；/>和/>分别表示原始和重构表达向量中所有表达量的标准差。
4、细胞距离矩阵的全局皮尔逊相关系数(cellDistPCC)
计算cellDistPCC之前要先计算细胞间的欧式距离(cellDistance)。其cellDistance和cellDistPCC的计算方法分别如公式(16)和公式(17)：
cellDistance是两个细胞之间的距离。其中x和y分别表示细胞i和细胞j中的表达向量；通过对矩阵X中所有细胞进行两两计算，得到一个有n*(n-1)/2个元素的一维距离向量d，同理得到的一维距离向量/>其中di和/>分别表示原始和重构矩阵的距离向量中的单个元素i的距离；μ和/>分别表示原始和重构表达矩阵中所有细胞间距离的平均值；O和/>分别表示原始和重构表达矩阵中所有细胞间距离的标准差。最终对得到的两个距离向量di和/>求皮尔逊相关性系数，即cellDistPCC。
5、调整后的兰德系数(ARI)
用于聚类模型的性能评估，其值域为[0，1]。
计算调整后的兰德系数需先计算兰德系数(RI)。RI和ARI的计算方法如公式(18)和(19)：
给定n个对象集合S＝{O1，O2，On}，假设T＝{a1，a2，an}和P＝{b1，b2，bn}表示表示S的两个不同划分，其中T为真实标签，P为聚类结果。设有四个统计量：
a为在T中为同一类且在P中也为同一类的数据点对数；
b为在T中为同一类且在P中不为同一类的数据点对数；
c为在T中不为同一类且在P中为同一类的数据点对数；
d为在T中不为同一类且在P中也为不同一类的数据点对数。
(五)算法工具的参数设置
每种算法的参数设置如下：
gimVI：根据gimVI的官方分析指南来重构空间转录组数据。输出的结果应为原始矩阵，将gimVI的推断函数“model.get_imputed_values”的参数设置为normalized＝False，其余均为默认参数。
SpaGE：根据SpaGE在GitHub代码存储库中的分析指南来重构空间转录组数据。用于训练集的基因数目(nGenes)小于200时，参数n_pv＝nGenes/2，否则，n_pv＝100，其余均为默认参数。
Tangram：根据Tangram在GitHub代码存储库中的分析指南代码存储库中的分析指南来重构空间转录组数据。将函数“tg.map_cells_to_space”的参数设置为：modes＝clusters，其余均为默认参数。
SMAF：根据SMAF在GitHub代码存储库中的分析指南代码存储库中的分析指南来重构空间转录组数据。SMAF将随机搜索30次采样矩阵并选择最优的一个，并将程序中的参数设置为：dictionary_size＝0.9dict_sparsity＝5，其余均为默认参数。
CSNet：解码层一共设置了3个残差块。每个残差块中的隐藏层神经元数量分别设置为：256,128,64。每个残差块中的隐藏层数设置2。其余重要参数设置为：learning_rate:0.001,batch_size＝128,max_epochs＝100。
(六)多节点分布式并行计算
由于任务较多，本发明采用多节点分布式计算的方法。进行测试的计算机群有25个计算节点，每个计算节点有2个Intel Xeon Scale 8358CPU(2.6GHz，48MB L3 Cache，32核)，1TB内存(DDR4 3200MHz)和8个NVIDIA A100显卡(SXM4 80GB显存)。
本发明采用wandb(Biewald L.Experiment tracking with weights andbiases,2020[J].Software available from wandb.com,2020,2(5).)工具来运行和管理多线程任务，这是一款用于机器学习实验跟踪的工具，可以用于多节点并行运算。wandb的主要功能包括：
1、多节点任务分发：可以将任务分配给多个节点并行执行。
2、程序执行与追踪：wandb会自动记录你的程序执行情况，并在实验结束后生成详细的日志，方便你进行后续分析。
3、网页可视化管理：可以在网页上实时查看实验的进度和结果，包括数据和图表等。
4、自动化日志记录：wandb可以自动记录你的实验结果和模型性能等信息，方便你进行后续分析。
首先，需要在服务器上搭建一个自己的网页管理平台，并通过秘钥与进行计算任务的服务器进行关联。然后在程序中调用wandb.init函数即可连接wandb网页管理平台并记录超参数，结合wandb.log将模型训练过程中变化的损失值和评估指标记录到日志系统中。最终，可以直接在网页端查看程序运行进度，并对记录的各类指标和超参数进行可视化分析。
wandb还提供了一种超参数调优工具sweep，它可以根据用户设置的参数列表自动分发任务。首先，需要设置一个config.yaml参数文件，其中包括一些参数列表(数据集，采样矩阵稠密度，细胞数等)；运行wandb sweep config.yaml命令将返回一个sweep_id；最后执行wandb agent sweep_id命令即可根据参数组合跑完所有任务。wandb允许在多节点中运行同一个sweep_id任务，且允许在单个节点中运行多次。该功能通过wandb任务管理平台对不同节点不同程序发来的执行参数进行扫描，保证相同参数的任务不会被重复执行，均衡了各个节点的计算资源。在5个节点中各自运行10个程序，即可将运行时间缩减到原来的1/50。
(七)五种重构算法的性能表现
由上文可知，本发明开发了一个只需用少量测量值即可重构空间转录组的模型CSNet，该模型是一个自编码器，由编码器和解码器组成。编码器将n个基因的表达量X压缩成m个基因集的表达量Y(n＞＞)，即Y＝AX，其中A为采样矩阵。解码器是一个ResNet结构的神经网络层，将m个基因集的表达量Y重构成n个基因的表达量模型需经过两次训练，第一次训练中同时训练编码器和解码器，该步骤可获得训练后的采样矩阵A第二次训练中冻结编码层，即固定了编码层的输出结果，只训练解码器，优化模型的重构性能。Y是多个基因集的复合测量值，通过采样矩阵A，可以将基因表达量X压缩为复合测量值Y，只需将采样得到的Y输入解码器即可重构出完整的基因表达矩阵。
为了验证CSNet模型的性能优越性，本发明选择了3种基于单基因测量的插值算法(SpaGE、Tangram和gimVI)与一种基于复合测量的压缩感知算法SMAF进行性能比较。
结果表明，本发明开发的模型CSNet在各个指标中性能表现最好，且压缩感知方法(CSNet，SMAF)要普遍优于单基因插值方法(SpaGE、Tangram和gimVI)。Tangram在genePCC指标上要比其余两种单基因插值方法表现更好。cellDistPCC指标反应了细胞邻居关系的保守性，压缩感知方法在该指标上的表现远优于单基因插值方法。这说明单基因插值的泛化能力弱，在处理外推问题时更容易过拟合，而压缩感知方法能更好地保留细胞间的邻居关系。此外，CSNet在任一指标中的表现也要优于现有的压缩感知算法SMAF，这表明CSNet比SMAF具有更强的学习能力和泛化能力(参见附图7)。
本发明通过下采样数据来改变原有数据集的细胞量，测试了5种工具在细胞占比分别为1％、5％、20％和100％时的性能。结果表明，CSNet的各个指标随着数据量的增加有明显提升，而其他工具在有的指标中有所提升、在有的指标中基本不变或有所下降(如：Tangram在genePCC指标中上升、在cellDistPCC指标下降)。且CSNet只有在训练集细胞量为1％时，它在genePCC指标上的表现不如SMAF，CSNet在cellPCC指标上的中位数也要低于SMAF；随着细胞量的增加，CSNet的性能逐渐超越SMAF。这可能是由于CSNet模型容量更大，能够拟合更复杂的函数，但是在数据量较少时容易导致过拟合；其他工具的容量较小，当数据集增加时，由于学习能力不足而导致欠拟合，从而性能变化不大。这说明其他算法更适用于小数据集的训练任务，而深度学习模型能够应用到大数据集的训练任务(参见附图8)。
由于不同的数据集的基因数不同，本发明选择每组数据集的基因数在汉明编码方案(如表5所示)中对应的测量数作为基数m，测试了5个算法工具在测量数分别为m/3、m/2、m、2m和3m的条件下的性能。结果表明，所有工具的各个指标随着测量数的增加有明显提升。其中，当测量数提高到原来3倍时，Tangram在genePCC指标的表现与CSNet相当。此外，压缩感知方法(CSNet和SMAF)在测量数较少的情况下，在cellDistPCC指标上的表现仍然较好，这说明哪怕只测量很少的基因集，压缩感知方法即使无法较好地重构出基因表达矩阵，但仍然能捕获到细胞的低秩结构信息，从而很好地保留细胞间的邻居关系。当测量数增加时，CSNet在cellDistPCC指标上的表现略微变差，这可能是由于CSNet训练具有随机性，不确定性较大(参见附图9)。
综上，本发明开发了一个基于深度压缩感知的空间转录组重构算法，并证明了比现有的其他工具的重构性能更好，这为本发明开发基于深度压缩感知的空间转录组技术打下了基础。并且，基于复合测量的压缩感知算法即使只使用较少的测量数重构数据，在cellDistPearson指标中表现仍然较好，这说明该方法只需极少量的测量便可还原细胞间的邻居关系。
实施例4：基于深度压缩感知的空间转录组实验方法
基于深度压缩感知的空间转录组实验技术是基于算法来重构空间转录组的新技术，简称为csFISH技术。其主要步骤包括模型训练、组织成像和数据重构三个步骤(如附图10所示)。首先，用单细胞数据集训练模型，得到采样方案。然后根据采样矩阵的结果将探针引物混合，设计每一次测量的引物池，将样本切片并用设计好的引物池中的荧光标记探针进行染色，通过低分辨率荧光显微镜对样本进行快速成像，并对图像数据进行滤波和细胞分割等处理，得到复合测量矩阵。最后，利用深度学习模型重构空间转录组数据。在本发明中，获得了测量小鼠海马体区域的113个基因的空间转录组图谱。
(一)单细胞数据预训练模型
在成像实验之前，需要训练好模型。首先，从NeMO和GEO数据库中获取4个基于smart-seq技术(Picelli S,Faridani O R,et al.Full-length RNA-seq from single cells using Smart-seq2[J]..Nature protocols,2014,9(1):171-181.)的鼠脑单细胞转录组数据集(如表6所示)，其中3个数据集用于构建训练集和验证集，剩下1个数据集作为测试集。然后进行模型调优，根据其皮尔逊相关系数指标和聚类指标ARI选择最优的模型参数和采样矩阵。最后，保存模型并根据模型输出的采样方案进行荧光原位杂交成像实验。所有smart-seq单细胞转录组数据集的物种来源、细胞数目和基因数目等注释信息如表7所示。
表6鼠脑单细胞转录组数据集信息
表7本发明中4个鼠脑单细胞转录组数据集信息
(二)图像数据获取
根据模型导出的采样方案设计实验所需引物池，并完成样本和其他实验试剂的制备，然后进行荧光原位杂交成像实验。本发明使用了共聚焦显微镜进行成像实验，其中使用20倍镜采集的图像像素大小为0.32μm×0.32μm，使用40倍镜采集的图像像素大小为0.16μm×0.16μm。实验需要进行多轮成像且一轮成像中有多个荧光通道，后续需要将图像对齐拼接。而由于使用低倍镜进行成像，得到的成像数据会较为模糊且有背景和随机噪声，需要对图像进行滤波将噪声去除。
(三)深度压缩感知模型重构空间转录组
在使用深度压缩感知模型重构空间转录组和进行下游分析之前，需要从图像中提取测量信息。本发明运行了一系列图像处理步骤，以对不同的荧光通道，不同的成像轮次和不同的视野中的图像进行归一化，拼接，对齐和分割。具体来说，本发明首先在z轴上进行了最大强度投影，然后使用DAPI通道在每一轮成像中拼接视野(使用Ashlar软件(Muhlich JL,Chen Y-A,Yapp C,et al.Stitching and registering highly multiplexed whole-slide images of tissues and tumors using ASHLAR[J].Bioinformatics,2022,38(19):4613-4621.)。通过宽度为8的中值滤波器，对每个轮次和每个通道的图像进行平滑处理。接着，使用Starfish管道(Axelrod S,Carr A,Freeman J,et al.Starfish:Opensource image based transcriptomics and proteomics tools[J].J.Open SourceSoftw,2018,6:2440.)进行点识别和背景去除。详细来说，调用了python图像处理库skimage中的peak_local_max函数，并分别设置参数“min_distance＝2和min_distance＝0.2来进行两轮点识别，并保留两轮点识别的结果，从而在信号拥挤的区域也能识别出足够密度的点对象。之后通过只保留点周围2个像素的值，并将其他像素的值置为0的方法来进行背景去除。去除背景后，将相同的缩放参数(固定的最大和最小强度阈值)应用到不同通道和轮次的图像上来进行重新缩放，调整亮度和对比度。对于细胞分割，使用基于深度学习的DeepCell memser模型(Greenwald N F,Miller G,Moen E,et al.Whole-cellsegmentation of tissue images with human-level performance using large-scaledata annotation and deep learning[J].Nat Biotechnol,2022,40(4):555-565.)对第一轮的DAPI通道的图像进行细胞核分割和分割结果修正，并计算图像掩膜，具体的细胞核分割参数为image_mpp＝0.32”，分割结果修正参数为“maxima_threshold＝0.001，maxima_model_smooth＝0，interior_model_smooth＝10，interior_threshold＝0.02，small_objects_threshold＝20，fill_holes_threshold＝20,radius＝1，pixel_expansion＝0”，使用每个细胞核分割掩膜中的像素强度和进行统计，即可得到每个细胞中不同轮复合测量的值，这些复合测量值被用于最终的解压缩步骤。
对经过图像处理得到的复合测量矩阵进行预处理(归一化和对数转换)，预处理后的复合测量矩阵作为输入数据传入训练好的模型的解码器中，解码器将根据输入推断出基因表达矩阵，结合细胞分割和识别得到的空间注释信息，最终得到单细胞精度的空间转录组数据。
(四)实验结果分析
如上所述，本发明挑选了113个基因进行实验，并为这113个基因设计了探针序列。本发明在同一个海马组织中同时进行csFISH和smFISH实验，smFISH的结果将作为金标准。csFISH技术将测量16个基因集来重构113个基因的表达量；smFISH技术则每次单独检测一个基因，一共测量113次。开始实验之前，需要训练模型获得采样矩阵。本发明收集了基于4个smart-seq技术的鼠脑单细胞转录组数据集(如表6所示)。其中3个作为训练集和验证集，剩下一个作为测试集。该测试集数据带有细胞亚型、实验批次等注释信息，本发明通过聚类指标ARI筛选最优模型。
本发明使用构造的鼠脑单细胞转录组数据集上进行训练和测试。得到了5种工具的重构数据后，进行聚类比较。结果表明，使用原始数据的113个基因进行聚类，其聚类结果与细胞亚型注释的ARI为0.68，SpaGE、gimVI、Tangram、SMAF和CSNet与细胞亚型注释的ARI分别为0.09、0.17、0.18、0.35和0.59。基于单基因插值的重构算法的性能明显差于基于压缩感知的重构算法，SpaGE的效果最差，没有准确分出任何亚群且分出许多小亚群。gimVI和Tangram性能要好一点，能分出三个亚群。SMAF则能分出四个亚群。CSNet性能最好，基本分出所有大亚群，数量较少的亚群则与原始数据一样聚在一起，这可能是由于这些小亚群的细胞数较少或选取基因数较少导致小亚群无法分开。这表明了CSNet重构的数据能够准确恢复基因表达量并保留细胞间的邻居关系(参见附图11)。
在验证了CSNe在鼠脑单细胞数据集上的重构性能后，保存训练好的模型并导出采样矩阵，根据采样矩阵设计和合成每轮成像的引物池。在进行荧光原位杂交成像实验之前，需要准备样本切片和调整好光学成像仪器。样本需要切成10μm的薄片以便探针引物能与目标mRNA转录本顺利杂交。
在完成实验准备工作后，开始荧光原位杂交成像实验。成像显微镜有4个荧光通道，所以每轮成像能测量4次。csFISH实验使用20倍物镜进行显微成像，只需进行4轮成像；smFISH实验使用40倍物镜进行荧光原位杂交成像实验，要进行29轮成像。由于一个完整的组织最终成像了多个视野，且每个视野成像了多轮通道，需要对其进行对齐和拼接。由于不同的荧光通道的光强不同，还需要对光强进行归一化。可以看到，smFISH和csFISH的成像结果都存在光污染(如附图12A所示)，因此需要对每一个视野进行滤波，去除背景和高斯噪声，得到一张高信噪比的图像(如附图12B所示)。
对图像进行预处理后，需要对图像进行细胞分割，得到细胞掩膜，然后才能统计细胞内的转录本数量。在进行成像时，为确定一个细胞的掩膜和位置，需要增加一个对细胞核进行DAPI染色的通道(如附图13A所示)。使用图像分割工具对该通道图像的细胞核进行图像分割，得到细胞核掩膜，然后调整参数控制细胞核掩膜向外膨胀得到细胞掩膜，最后将所有通道对齐，即可将mRNA转录本分配给对应细胞当中(如附图13B所示)。
smFISH实验的结果较为清晰，能够准确识别出单个mRNA转录本的空间位置，通过点识别方法统计单个细胞掩膜内单次测量到的mRNA转录本数作为此次测量的基因表达量，将所有视野中的细胞在每一次测量的转录本数进行统计，即可得到基因矩阵。
csFISH实验获得图像则较为模糊，多个mRNA转录本的荧光可能在一个位置上发生重叠，无法通过准确识别每一个转录本的具体位置来进行计数。于是，本发明通过统计细胞内的总光强值来代表该细胞在单次测量中测得基因集总表达量，将所有视野中的细胞在每一次测量的光强进行统计，即可得到复合测量矩阵。然而发现部分视野可能由于基因集表达量过高从而超出所能测得的光强最大值，即光强过曝。这可能会对后续重构的结果造成不利影响。
得到复合测量矩阵后，首先对数据进行质量控制，去除一些由于实验原因而没有捕获到mRNA转录本的细胞，然后对数据进行归一化处理和对数转换，最后将结果输入之前保存的模型，重构出基因表达矩阵。
过滤掉质量较差的细胞，然后对数据进行归一化和对数转换，最后进行无监督聚类并查看其聚类效果与基因表达情况是否符合预期。
使用leiden算法对细胞进行聚类，得到20个小亚群，本发明根据一些标记基因的表达情况对细胞进行注释，并把20个小亚群合并成11个亚群。其中CA1 neuron和CA2neuron是海马体CA区神经元细胞，DG neuron是海马体DG区神经元细胞，Ex neuron是兴奋性神经元细胞，Ctgf+Ex neuron是Ctgf阳性兴奋性神经元细胞，In neuron是抑制性神经元细胞，Astrocyte是星形胶质细胞，Microglia是小胶质细胞，Oligo是少突胶质细胞，Vascular是血管细胞(参见附图14A)。
神经细胞都在UMAP分布的左上部分，“Gad1，Scl17a7”都是神经细胞的标记基因，主要在左上部分表达；“Slc17a6”是兴奋性神经元的标记基因，但其分布较为弥散，没有集中兴奋性神经元细胞，可能是由于该基因恢复效果较差；在还有一小群特异性表达“Ctgf”的Ctgf阳性兴奋性神经元；“Gja1，Cldn11”则是星形胶质细胞和少突胶质细胞的标记基因，主要在下半部分表达；左下部分有一群小胶质细胞表达“Cd14”；中间则有一小群血管细胞表达“Esam”；还有一群细胞不表达任何标记基因，可能是由于实验误差或部分视野图像处理步骤没有做好，从而导致重构的基因表达结果有误(参见附图14B)。
经过初步判断csFISH实验测得的大部分细胞的基因表达量能够鉴定出相应的亚群和标记基因。后面将对smFISH实验结果进行聚类，并比较两个技术的细胞亚群和基因表达的空间分布是否一致。
对基于smFISH技术测得的空间转录组数据进行同样的预处理和聚类操作后，得到了smFISH数据的细胞亚群注释信息。本发明通过比较两种技术得到的细胞亚群和基因表达的空间分布是否一致来验证csFISH的准确性。
在csFISH和smFISH实验数据中，三种胶质细胞在海马体中都呈弥散分布；并且CA神经元和DG神经元都集中分布在海马体CA区和DG区；兴奋性神经元主要在右上角；抑制性神经元则呈弥散分布；而Ctgf阳性兴奋性神经元细胞则是下方某一区域成线状分布。两种技术测得的数据的细胞亚型分布较为一致，说明csFISH重构的数据较好地保留了细胞的低秩结构信息(参见附图15A和15C)。
“Cdn11，Gja1，Cd14，Esam”是胶质细胞和血管细胞的标记基因，它们的表达量在空间中弥散分布，Slc17a7在海马体CA区和DG区以及最上方和最下方中表达较高；Ctgf在两组数据中的表达有特异性分布且位置基本一致；在smFISH中，Gad1的表达量在CA区的分布存在不一致的区域，可能是由于该区域部分视野光强过曝导致；Slc17a6是兴奋性神经元的标记基因，兴奋性神经元在两组数据中分布较为一致，但是Slc17a6在smFISH实验数据中的兴奋性神经元区域表达则明显比csFISH更高，可能是由于模型重构的表达量较为平滑(参见附图15B和15D)。
对两组数据的所有细胞的基因表达量分别求平均值，再计算两组基因平均表达量的皮尔逊相关性系数，其相关性系数达到0.875，说明两组数据具有强相关性，结果比较吻合。再分别对每个亚群求平均值并计算两组数据中亚群间的相关性，两组数据中相同的亚群相关性较高，亚群的分类结果较为一致(参见附图15E和15F)。
csFISH和smFISH的细胞亚型的分类与空间分布较为一致，有少部分区域和基因的表达量不太相似。这些结果说明了csFISH技术基本能测量出一个质量较高的基因表达矩阵。
综上，本发明成功开发了一项新的空间转录组技术，该技术通过深度压缩感知算法在20倍镜的条件下通过16次测量重构出了小鼠海马体的113个基因空间转录组，并且还在40倍镜的条件下使用smFISH技术测量了113个基因作为金标准数据。通过相关性分析可以看到本发明开发的技术与smFISH的测量结果基本一致，但是csFISH技术在更低倍镜和更少的测量数条件下进行实验，使用20倍镜在每一次测量中的成像速度要比使用40倍镜快4倍，16次测量又比113次测量节约了6/7的时间，因此csFISH技术的成像速度是smFISH技术的28倍。
smFISH技术使用点识别方法来统计mRNA转录本数量，但当mRNA转录本数量较多时，则会导致光学拥挤，进而导致点识别方法失效。而csFISH技术则是通过统计总光强来估算mRNA转录本数量，该方法能够解决由于光学拥挤而无法统计基因表达量的问题，适用于测量表达量较高的基因。
以上结合了优选的实施方式对本申请进行了说明，不过这些实施方式仅是范例性的，仅起到说明性的作用。在此基础上，可以对本申请进行多种替换和改进，这些均落入本申请的保护范围内。

Claims (20)

1.一种自编码器，其是将输出恢复成输入的神经网络结构，包含编码器和解码器，其中：
编码器是由一个线性压缩层组成的神经网络，将高维数据转换为低维表示，实现数据的压缩；
解码器由多个全连接层组成，其采用残差网络ResNet作为网络结构，将压缩后的数据重构成原始输入数据，实现数据的重构。
2.根据权利要求1所述的自编码器，其通过最小化输入x和输出之间的误差同时完成压缩和重构两个步骤，其中，
编码器的线性压缩层可以将基因表达矩阵X通过与采样矩阵A进行矩阵乘法进行线性变换得到复合测量矩阵Y，其中A是m×g的二维矩阵，X是g×n的二维矩阵，Y是m×n的二维矩阵，编码器的输出可以表示为如下公式，
Y＝AX
通过所述公式，编码器将高维的输入数据转换为低维的线性表示，在数据采集时实现压缩；
解码器将编码器的输出结果Y作为解码器的输入，通过如下公式，重构基因表达矩阵
。
3.一种基于深度学习的压缩感知模型CSNet，其包含权利要求1或2所述的自编码器。
4.一种基于深度学习的压缩感知方法，其将深度学习技术应用到压缩感知中，包括以下步骤：
(1)构建权利要求1或2所述的自编码器，或者权利要求3所述的模型CSNet；
(2)采用步骤(1)构建的自编码器或模型CSNet进行模型训练和推理，在进行模型训练中，通过反向传播算法更新网络权重以最小化损失函数。
5.根据权利要求4所述的方法，其中步骤(2)包含编码器和解码器联合训练、解码器单独训练和解码器推理三个过程。
6.根据权利要求4或5所述的方法，步骤(2)中的模型训练包括训练步骤、验证步骤以及测试步骤。
7.根据权利要求4-6中任一项所述的方法，其使用基于PyTorch开发的Pytorch-Lightning高级框架来实现该压缩感知算法。
8.根据权利要求4-7中任一项所述的方法，在进行模型训练之前，还包括对数据进行处理的步骤，优选对复合测量矩阵进行归一化和对数转换。
9.权利要求1或2所述的自编码器或者权利要求3所述的模型CSNet在基于深度学习的压缩感知方法中的应用。
10.一种基于深度压缩感知的空间转录组方法，包括以下步骤：
(1)模型构建，构建权利要求1或2所述的自编码器，或者权利要求3所述的模型CSNet；
(2)模型训练，获得训练好的模型和采样矩阵；
(3)组织成像，进行荧光原位杂交成像实验得到数据图像；
(4)数据重构，从图像中提取复合测量矩阵并输入压缩感知模型中重构基因表达矩阵。
11.根据权利要求10所述的方法，其中步骤(2)通过单细胞数据集训练模型，得到采样方案。
12.根据权利要求10或11所述的方法，其中步骤(2)还包括根据皮尔逊相关系数指标和聚类指标ARI选择最优的模型参数和采样矩阵。
13.根据权利要求10-13中任一项所述的方法，其中步骤(3)还包括对图像进行滤波去噪以及图像处理步骤，所述图像处理优选为对图像进行归一化、拼接、对齐和分割。
14.根据权利要求10-13中任一项所述的方法，其中步骤(4)还包括对对经过图像处理得到的复合测量矩阵进行预处理的步骤，所述预处理优选为进行归一化和对数转换。
15.权利要求1或2所述的自编码器或者权利要求3所述的模型在基于深度压缩感知的空间转录方法中的应用。
16.一种基于深度压缩感知的空间转录组装置，包括以下模块：
(1)模型构建模块，用于构建权利要求1或2所述的自编码器，或者权利要求3所述的模型CSNet；
(2)模型训练模块，用于获得训练好的模型和采样矩阵；
(3)组织成像模块，用于进行荧光原位杂交成像实验得到数据图像；
(4)数据重构模块，用于从图像中提取复合测量矩阵并输入压缩感知模型中重构基因表达矩阵。
17.根据权利要求16所述的装置，其中模型训练模块通过单细胞数据集训练模型，得到采样方案。
18.根据权利要求16或17所述的装置，其中模型训练模块根据皮尔逊相关系数指标和聚类指标ARI选择最优的模型参数和采样矩阵。
19.根据权利要求16-18中任一项所述的装置，其中组织成像模块还包括滤波去噪和图像处理模块，滤波去噪用于对图像进行滤波将噪声去除，图像处理用于对图像进行处理，所述图像处理优选为对图像进行归一化、拼接、对齐和分割。
20.根据权利要求16-19中任一项所述的装置，其中数据重构模块还包括复合测量矩阵预处理模块，用于对经过图像处理得到的复合测量矩阵进行预处理，所述预处理优选为进行归一化和对数转换。
CN202310514516.7A 2023-05-08 2023-05-08 一种基于深度压缩感知的空间转录组方法及装置 Pending CN116842996A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310514516.7A CN116842996A (zh) 2023-05-08 2023-05-08 一种基于深度压缩感知的空间转录组方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310514516.7A CN116842996A (zh) 2023-05-08 2023-05-08 一种基于深度压缩感知的空间转录组方法及装置

Publications (1)

Publication Number Publication Date
CN116842996A true CN116842996A (zh) 2023-10-03

Family

ID=88169460

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310514516.7A Pending CN116842996A (zh) 2023-05-08 2023-05-08 一种基于深度压缩感知的空间转录组方法及装置

Country Status (1)

Country Link
CN (1) CN116842996A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117476114A (zh) * 2023-12-25 2024-01-30 墨卓生物科技(浙江)有限公司 一种基于生物多组学数据的模型构建方法与系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117476114A (zh) * 2023-12-25 2024-01-30 墨卓生物科技(浙江)有限公司 一种基于生物多组学数据的模型构建方法与系统
CN117476114B (zh) * 2023-12-25 2024-04-05 墨卓生物科技(浙江)有限公司 一种基于生物多组学数据的模型构建方法与系统

Similar Documents

Publication Publication Date Title
Hou et al. Sparse autoencoder for unsupervised nucleus detection and representation in histopathology images
Dorkenwald et al. Automated synaptic connectivity inference for volume electron microscopy
Fuchs et al. Computational pathology: challenges and promises for tissue analysis
Baygin et al. Exemplar Darknet19 feature generation technique for automated kidney stone detection with coronal CT images
CN114730463A (zh) 用于组织图像分类的多实例学习器
CN113454733A (zh) 用于预后组织模式识别的多实例学习器
CN106096506B (zh) 基于子类类间判别双字典的sar目标识别方法
JP7427080B2 (ja) 細胞検出およびセグメンテーションのための弱教師ありマルチタスク学習
CN116580848A (zh) 一种基于多头注意力机制的分析癌症多组学数据方法
CN116113993A (zh) 使用少样本学习的组织学染色模式和伪影分类
CN114864003A (zh) 基于混合实验组和对照组单细胞样本的差异分析方法及系统
Kulikov et al. DoGNet: A deep architecture for synapse detection in multiplexed fluorescence images
He et al. SAR target recognition and unsupervised detection based on convolutional neural network
Razavi et al. MiNuGAN: Dual segmentation of mitoses and nuclei using conditional GANs on multi-center breast H&E images
CN116842996A (zh) 一种基于深度压缩感知的空间转录组方法及装置
CN115881232A (zh) 一种基于图神经网络和特征融合的scRNA-seq细胞类型注释方法
Barrera et al. Automatic normalized digital color staining in the recognition of abnormal blood cells using generative adversarial networks
Ternes et al. Me-vae: Multi-encoder variational autoencoder for controlling multiple transformational features in single cell image analysis
Rethik et al. Attention Based Mapping for Plants Leaf to Classify Diseases using Vision Transformer
Lee et al. MorphNet predicts cell morphology from single-cell gene expression
Ram et al. Combined detection and segmentation of cell nuclei in microscopy images using deep learning
Rozendo et al. Classification of non-Hodgkin lymphomas based on sample entropy signatures
CN115862746A (zh) 一种精准的单细胞多组学匹配数据生成方法
CN115131628A (zh) 一种基于分型辅助信息的乳腺图像分类方法及设备
Shakya et al. Bacterial Image Segmentation through Deep Learning Approach

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination