CN117437973A - 一种单细胞转录组测序数据插补方法 - Google Patents

一种单细胞转录组测序数据插补方法 Download PDF

Info

Publication number
CN117437973A
CN117437973A CN202311763792.3A CN202311763792A CN117437973A CN 117437973 A CN117437973 A CN 117437973A CN 202311763792 A CN202311763792 A CN 202311763792A CN 117437973 A CN117437973 A CN 117437973A
Authority
CN
China
Prior art keywords
matrix
cell
gene
original
feature matrix
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.)
Granted
Application number
CN202311763792.3A
Other languages
English (en)
Other versions
CN117437973B (zh
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.)
Qilu University of Technology
Original Assignee
Qilu University of Technology
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 Qilu University of Technology filed Critical Qilu University of Technology
Priority to CN202311763792.3A priority Critical patent/CN117437973B/zh
Publication of CN117437973A publication Critical patent/CN117437973A/zh
Application granted granted Critical
Publication of CN117437973B publication Critical patent/CN117437973B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G16B40/20—Supervised data analysis
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G16B40/30—Unsupervised data analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Public Health (AREA)
  • Databases & Information Systems (AREA)
  • Artificial Intelligence (AREA)
  • Bioethics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Chemical & Material Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种单细胞转录组测序数据插补方法，涉及单细胞RNA测序数据处理领域。本发明包括以下步骤：S1：构建原始基因表达矩阵V，并对其进行归一化处理；S2：识别原始基因表达矩阵V中潜在的缺失值，构建系数矩阵P；S3：利用非负矩阵分解算法分解原始基因归一化矩阵，得到细胞潜在特征矩阵Q和基因潜在特征矩阵H，并分别对其进行L1正则化约束，而后对细胞潜在特征矩阵Q应用图正则化约束，而后，迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，直至收敛，得预测矩阵N；S4：获取插补后的矩阵。本发明所述方法插补后的测序数据，在下游的细胞聚类和细胞轨迹重建领域均有较佳表现。

Description

一种单细胞转录组测序数据插补方法
技术领域
本发明涉及单细胞RNA测序数据处理领域，尤其涉及一种单细胞转录组测序数据插补方法。
背景技术
单细胞RNA测序技术在细胞水平上提供了高分辨率的基因表达信息，不过，由于实验条件或技术限制，单细胞RNA测序数据中可能存在一些细胞的基因表达值缺失。而且，由于单细胞RNA测序数据还具有噪声、高度稀疏和海量的特点，使得现有的单细胞RNA测序插补方法很难仅针对基因表达值缺失的部分进行有效地插补，进而影响单细胞RNA测序数据在聚类以及重建细胞轨迹等下游任务中的效果。而聚类以及重建细胞轨迹等下游任务对于人类疾病的发病机制及其治疗又起着至关重要的作用。因此，如何提高单细胞RNA测序插补方法的插补效果一直是科研人员的关注重点。
现有技术中，虽然也有针对基因表达值缺失的部分进行插补的方法，不过，现有的针对基因表达值缺失的部分进行插补的方法，没有对测序技术导致的基因表达值缺失和生物变异导致的缺失进行区分，这就使得针对基因表达值缺失的部分进行插补时会导致原有的因生物变异导致的缺失值也被插补，造成过度插补，影响插补的效果。为解决上述问题，本申请提出了一种仅对测序技术导致的基因表达值缺失进行插补的单细胞转录组测序数据插补方法。
发明内容
为了弥补现有技术的不足，本发明提供了一种单细胞转录组测序数据插补方法。
本发明的技术方案为：
一种单细胞转录组测序数据插补方法，包括以下步骤：
S1：筛选变异特征基因，根据筛选出的变异特征基因及与其相对应的细胞类型构建原始基因表达矩阵V；而后，对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵；
S2：识别原始基因表达矩阵V中潜在的缺失值，构建用于仅对由测序技术原因所产生的技术零值的位置进行插补的系数矩阵P；
S3：利用非负矩阵分解算法分解原始基因归一化矩阵，得到细胞潜在特征矩阵Q和基因潜在特征矩阵H，而后，对细胞潜在特征矩阵Q和基因潜在特征矩阵H进行L1正则化约束，减小过拟合的风险，而后，对细胞潜在特征矩阵Q应用图正则化约束捕捉细胞在潜在空间中的几何关系，而后，迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，直至收敛，结束迭代更新过程，得到最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H，并利用最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H得到与原始矩阵相近的预测矩阵N；
S4：根据系数矩阵P、原始基因归一化矩阵以及与原始矩阵相近的预测矩阵N计算得到插补后的矩阵/>。
优选地，步骤S1具体为：
S1-1、从目标数据集中获取细胞类型和与其相对应的基因数据，而后，使用Seurat程序包对目标数据集中的基因数据进行裁剪，并保留基因数据裁剪处理后得到的一定数量的变异特征基因；
S1-2、利用步骤S1-1保留的变异特征基因，并结合其相适应的细胞类型的名称构成原始基因表达矩阵V；
S1-3、对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵。
优选地，步骤S1-1中，当目标数据集由多个不同数据集组成时，不同数据集基因数据裁剪处理后，针对不同数据集均保留一致数量的变异特征基因；
优选地，步骤S2具体包括如下步骤：
S2-1、使用谱聚类将原始基因表达矩阵V中的细胞进行大致聚类，分别得到k个细胞亚群；
S2-2、计算原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，基于每一个细胞亚群中的变异特征基因的零表达率/>、平均表达值/>和方差/>计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；当置信度β小于阈值t时，则将该变异特征基因的元素值为零的值，视为因为生物变异所产生的零值，在与原始基因表达矩阵V相对应的系数矩阵P中的对应位置的元素值设置为1；
S2-3、遍历原始基因表达矩阵V中每条变异特征基因的元素值，如果元素值大于零就将与该元素值在原始基因表达矩阵V中位置相对应的系数矩阵P中的位置处的元素值设置为1，表示为生物表达值；而后，再将原始基因表达矩阵V中剩余的元素值为零值的位置，在系数矩阵P中的对应位置处的元素值设置为0，表示因为技术原因所产生的技术零值，需要进行插补，如此便得到了完整的系数矩阵P。
优选地，步骤S2-2中，利用原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>计算用于区分变异特征基因元素为生物零和技术零的置信度β，如式(2)所示：
(2)
式(2)中，表示在第k个细胞亚群中第i个基因的置信度，/>表示在第k个细胞亚群中第i个基因的零表达率，/>表示在第k个细胞亚群中第i个基因的表达的平均值，表示在第k个细胞亚群中第i个基因的表达的方差。
优选地，步骤S3具体为：
S3-1、将原始基因归一化矩阵利用非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q和基因潜在特征矩阵H；
S3-2、对细胞潜在特征矩阵Q和基因潜在特征矩阵H施加L1正则化约束，得到L1正则化Y；
S3-3、计算原始基因归一化矩阵中相邻两列细胞的余弦距离，构建相似性矩阵S，并根据相似性矩阵S计算图拉普拉斯矩阵/>，根据图拉普拉斯矩阵/>计算相似性矩阵S的拉普拉斯矩阵L，而后对细胞潜在特征矩阵Q应用图正则化约束计算拉普拉斯矩阵L中矩阵特征值总和M；为了了解潜在的细胞亚群的结构，对细胞潜在特征矩阵Q应用图正则化约束，获得拉普拉斯矩阵L中矩阵特征值总和M，来学习潜在的细胞亚群的结构捕捉细胞在潜在空间中的几何关系；
S3-4、以步骤S3-1得到的细胞潜在特征矩阵Q和基因潜在特征矩阵H以及步骤S3-3得到的相似性矩阵S为基础，使用最小二乘法迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，在迭代交替更新的过程中，同时计算更新后的细胞潜在特征矩阵Q的相似性矩阵S来更新对更新后的细胞潜在特征矩阵Q的图正则化约束，利用跟新后的相似性矩阵S来更新图拉普拉斯矩阵以及更新细胞潜在特征矩阵Q，而后利用跟新后的图拉普拉斯矩阵/>来更新拉普拉斯矩阵L，而后利用跟新后的拉普拉斯矩阵L来更新拉普拉斯矩阵L中矩阵特征值总和M；而后利用更新后的矩阵特征值总和M更新近似矩阵N与原始基因归一化矩阵/>之间的F范数误差E，直至F范数误差E收敛，停止迭代，得到最终的细胞潜在特征矩阵Q、基因潜在特征矩阵H以及与原始矩阵相近的预测矩阵N；所述近似矩阵N是通过细胞潜在特征矩阵Q和基因潜在特征矩阵H相乘得到的。
优选地，步骤S3-1中，将原始基因归一化矩阵利用非负矩阵分解算法随机分解为细胞潜在特征矩阵Q和基因潜在特征矩阵H的计算方式，如式(4)和式(5)所示：
(4)
(5)
式(4)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>表示细胞潜在特征矩阵H的转置；/>是哈达玛积；/>是系数超参数。
式(5)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>是细胞潜在特征矩阵Q的转置；/>是哈达玛积；/>是系数超参数。
优选地，步骤S3-2中，L1正则化Y是通过对细胞潜在特征矩阵Q和基因潜在特征矩阵H的F范数之和乘以系数获得的。
优选地，步骤S3-4中，F范数误差E收敛的条件是F范数误差E小于10-5。
优选地，步骤S4中，插补后的矩阵的计算方式，如式(16)所示；该步骤S4能够只针对由于技术原因所产生的技术零值进行插补；
(16)
式(16)中，是插补后的矩阵，P表示系数矩阵，/>是哈达玛积，N表示预测矩阵；/>表示原始基因归一化矩阵。
与现有技术相比，本发明的有益效果如下：
本申请所述单细胞转录组测序数据插补方法，通过原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值和方差计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；当置信度β小于阈值t时，则将该变异特征基因的元素值为零值的情况视为因为生物变异所产生的零值，并将原始基因表达矩阵V中剩余的元素值为零值的情况，视为因为测序技术原因所产生的技术零值，确定因为测序技术原因所产生的技术零值需要进行插补，并基于区分的原始基因表达矩阵V中因为生物变异所产生的零值以及因为测序技术原因所产生的技术零值情况，构建用于仅对由测序技术原因所产生的技术零值位置进行插补的系数矩阵P，可有效避免过度插补；在构建系数矩阵P后，再使用图正则化，捕捉细胞在潜在空间中的几何关系，最后使用非负矩阵分解实现对单细胞RNA测序数据中因为测序技术原因所产生的技术零值进行插补。本申请所述的单细胞转录组测序数据插补方法，得到的插补后的单细胞RNA测序数据在下游分析任务比如细胞聚类和重建细胞轨迹等中均取得了最好的结果；比如，在li数据集中，本发明所述单细胞转录组测序数据插补方法插补后进行细胞聚类得到的ARI值和NMI值分别为0.9927305、0.9918196，相对于DrImpute插补方法得到的ARI值和NMI值分别提高了0.045%和0.122%；在liver数据集中，本发明所述单细胞转录组测序数据插补方法插补后进行细胞聚类得到的ARI值和NMI值分别为0.7835690、0.8123501，相对于ALRA插补方法得到的ARI值和NMI值分别提高了9.726%和3.363%；而本申请所述的单细胞转录组测序数据插补方法在Deng数据集上取得的KOR和POS两个评价指标分别达到了0.79和0.93，相对于MAGIC插补方法分别提高了8.219%和2.198%；这说明了本申请所述的单细胞RNA测序数据插补方法对单细胞RNA测序数据插补效果更好，更适用于下游任务的分析。
附图说明
图1为本发明的总体流程图；
图2为本发明所述的单细胞转录组测序数据插补方法以及ALRA插补方法、DrImpute插补方法、DeepImpute插补方法、MAGIC插补方法、SAVER插补方法、SCDD插补方法、scImpute插补方法、scVI插补方法和WEDGE插补方法这九种现有的插补方法对li数据集执行插补后的UMAP可视化散点图；
图3为本发明所述的单细胞转录组测序数据插补方法以及ALRA插补方法、DrImpute插补方法、DeepImpute插补方法、MAGIC插补方法、SAVER插补方法、SCDD插补方法、scImpute插补方法、scVI插补方法和WEDGE插补方法这九种现有的插补方法对liver数据集执行插补后的UMAP可视化散点图；
图4为本发明的单细胞转录组测序数据插补方法以及ALRA插补方法、DrImpute插补方法、DeepImpute插补方法、MAGIC插补方法、SAVER插补方法、SCDD插补方法、scImpute插补方法、scVI插补方法和WEDGE插补方法这九种现有的插补方法对Deng数据集执行插补后的可视化细胞轨迹重建图。
具体实施方式
为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
一种基于单细胞转录组测序数据插补方法，其总体流程图，如图1所示，具体包括以下步骤：步骤S1：筛选变异特征基因，根据筛选出的变异特征基因及与其相对应的细胞类型构建原始基因表达矩阵V；而后，对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵，具体包括如下步骤：
S1-1、从目标数据集中获取细胞类型和与其相对应的基因数据，而后，使用Seurat程序包对目标数据集中的基因数据进行裁剪，并保留基因数据裁剪处理后得到的一定数量的变异特征基因；其中，当目标数据集由多个不同数据集组成时，不同数据集基因数据裁剪处理后，针对不同数据集均保留一致数量的变异特征基因；
具体到本实施例来说，本实施例目标数据集由li数据集、liver数据集和Deng数据集构成，本实施例中分别从li数据集、liver数据集和Deng数据集中获取细胞类型和与其相对应的基因数据，而后，分别使用Seurat程序包对li数据集、liver数据集和Deng数据集中的基因数据进行裁剪，并保留li数据集中基因数据裁剪处理后得到的前2000条的变异特征基因、保留liver数据集中基因数据裁剪处理后得到的前2000条的变异特征基因、保留Deng数据集中基因数据裁剪处理后得到的前2000条的变异特征基因；
本申请选用的li数据集为现有的数据集，获取的网址为：https://zenodo.org/records/7504311#.Y7ei0naZMuK；
本申请选用的liver数据集为现有的数据集，获取的网址为：https://zenodo.org/records/7504311#.Y7ei0naZMuK；
本申请选用的Deng数据集为现有的数据集，获取的网址为：https://zenodo.org/records/6832477#.Yw9HYBxByHs；
S1-2、利用步骤S1-1保留的变异特征基因，并结合其相适应的细胞类型的名称构成原始基因表达矩阵V；
具体到本实施例中，具体为：利用li数据集中基因数据裁剪处理后得到的前2000条的变异特征基因以及与其相适应的细胞类型的名称构成原始基因表达矩阵V1；利用liver数据集中基因数据裁剪处理后得到的前2000条的变异特征基因以及与其相适应的细胞类型的名称构成原始基因表达矩阵V2；利用Deng数据集中基因数据裁剪处理后得到的前2000条的变异特征基因以及与其相适应的细胞类型的名称构成原始基因表达矩阵V3；
其中，本实施例中原始基因表达矩阵V1、V2和V3的第一列均代表变异特征基因的名称，第一行均代表与变异特征基因相对应的细胞类型的名称，原始基因表达矩阵V1、V2和V3的元素是基因表达值，原始基因表达矩阵V1行数为2001，原始基因表达矩阵V1列数为562，原始基因表达矩阵V2行数为2001，原始基因表达矩阵V2列数为778，原始基因表达矩阵V3行数为2001，原始基因表达矩阵V3列数为287；为了展示原始基因表达矩阵V1的式样，本申请特提供原始基因表达矩阵V1的前十行、前十列，如表1所示；
表1
A549 A549 A549 A549 A549 A549 A549 A549 A549
HBG2 0 1 0 0 0 0 0 0 0
HBG1 0 0 0 0 0 0 0 0 0
IGHM 0 0 0 0 0 0 0 0 0
ACSM3 74.9995 370 0 0 69 0 0 120 32
IGLV6-57 0 0 0 0 0 0 0 0 0
IGLC2 0 0 0 0 0 0 0 0 0
IGLC3 0 0 0 0 0 0 0 0 0
PRAME 204 5 0 345 498 39 225 4 1
AC079466.1 0 0 2 0 0 0 0 0 0
S1-3、对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵，如式(1)所示：
(1)
式(1)中，log10是计算以10为底的对数，V为原始基因表达矩阵，表示原始基因归一化矩阵。
本实施例中利用式(1)所述的方式对原始基因表达矩阵V1、V2和V3分别进行归一化处理后，分别得到原始基因归一化矩阵、/>和/>。
步骤S2：识别原始基因表达矩阵V中潜在的缺失值，构建用于仅对由测序技术原因所产生的技术零值位置进行插补的系数矩阵P；具体包括如下步骤：
步骤S2-1、使用谱聚类将原始基因表达矩阵V中的细胞进行大致聚类，得到k个细胞亚群；其中，不同原始基因表达矩阵V的细胞亚群数量k与原始基因表达矩阵V相对应的数据集中的细胞类型种类数量是相同的；
S2-2、计算原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，利用式(2)计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；置信度β越高，表明基因表达的零值为测序技术产生的缺失值的可能性越大，这主要是由于在细胞亚群中具有高表达和低变异性的基因表达的零值通常被认为是测序技术产生的零值；而在具有高变异性的细胞亚群中始终以低至中等水平表达的基因，如果其缺失值表达值为零，则被认为反映了真正的生物学变异性。
当置信度β小于阈值t时，则将该变异特征基因的元素值为零的值，视为因为生物变异所产生的零值，在与原始基因表达矩阵V相对应的系数矩阵P中的对应位置的元素值设置为1；
(2)
式(2)中，表示在第k个细胞亚群中第i个基因的置信度，/>表示在第k个细胞亚群中第i个基因的零表达率，/>表示在第k个细胞亚群中第i个基因的表达的平均值，表示在第k个细胞亚群中第i个基因的表达的方差。
S2-3、遍历原始基因表达矩阵V中每条变异特征基因的元素值，如果元素值大于零就将与该元素值在原始基因表达矩阵V中位置相对应的系数矩阵P中的位置处的元素值设置为1，表示为生物表达值；而后，再将原始基因表达矩阵V中剩余的元素值为零值的位置，在系数矩阵P中的对应位置处的元素值设置为0，表示因为技术原因所产生的技术零值，需要进行插补，如此便得到了完整的系数矩阵P；
其中，上述系数矩阵P中各元素的确定关系，也如式(3)示意：
(3)
式(3)中，表示系数矩阵P中各位置的元素值，i表示系数矩阵P中的第i行，j表示系数矩阵P中的第j列，/>表示在第k个细胞亚群中第i个基因的置信度，/>表示原始基因表达矩阵V1、V2或者V3中第i行第j列的变异特征基因的元素值，otherwise表示原始基因表达矩阵V1、V2或者V3中除了元素值大于零以及置信度/>小于阈值t时的变异特征基因为零值之外的剩余元素值为零值的情况。
具体到本实施例来说，步骤S2具体为：
识别原始基因表达矩阵V1、V2和V3中潜在的缺失值，构建系数矩阵P1、P2和P3，具体包括如下步骤：
S2-1、使用谱聚类将原始基因表达矩阵V1、V2和V3中的细胞分别进行大致聚类，分别得到k1、k2和k3个细胞亚群；其中，k1与li数据集中细胞类型种类数量相同，即k1为9；k2与liver数据集中细胞类型种量数据相同，即k2为7；k3与Deng数据集中细胞类型种类数量相同，即k3为10；
S2-2、计算原始基因表达矩阵V1中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，利用式(2)计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；此时，式(2)中的t取值为0.1；
计算原始基因表达矩阵V2中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，利用式(2)计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；此时，式(2)中的t取值为0.1；
计算原始基因表达矩阵V3中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，利用式(2)计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；此时，式(2)中的t取值为0.6。
S2-3、遍历原始基因表达矩阵V1中每条变异特征基因的元素，如果元素大于零就将与该元素在原始基因表达矩阵V中位置相对应的系数矩阵P中位置的元素值设置为1，表示生物表达值；而后，再将原始基因表达矩阵V1中剩余的零值，在系数矩阵P1中的对应位置的元素值设置为0，表示因为技术原因所产生的技术零，需要进行插补；如此，便得到完整的系数矩阵P1；
本申请中系数矩阵P2以及系数矩阵P3的获取方式，与系数矩阵P1的获取方式相同。
步骤S3：利用非负矩阵分解算法分解原始基因归一化矩阵，得到细胞潜在特征矩阵Q 和基因潜在特征矩阵H，而后，对细胞潜在特征矩阵Q和基因潜在特征矩阵H进行L1正则化约束，减小过拟合的风险，而后，对细胞潜在特征矩阵Q应用图正则化约束捕捉细胞在潜在空间中的几何关系，而后，迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，至更新后的细胞潜在特征矩阵Q和更新后的基因潜在特征矩阵H相乘得到的预测矩阵N与归一化的原始矩阵/>之间的F范数E小于10-5时，结束迭代更新过程，得到最终的细胞潜在特征矩阵Q、最终的基因潜在特征矩阵H以及最终的预测矩阵N，最终的预测矩阵N即为与原始矩阵相近的预测矩阵N；具体包括如下步骤：
步骤S3-1：将原始基因归一化矩阵利用非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q和基因潜在特征矩阵H，其中，细胞潜在特征矩阵Q中每一列均表示细胞的潜在特征，基因潜在特征矩阵H中每一列均表示变异特征基因的潜在特征。将原始基因归一化矩阵/>利用非负矩阵分解算法随机分解为细胞潜在特征矩阵Q和基因潜在特征矩阵H的计算方式，如式(4)和式(5)所示：
(4)
(5)
式(4)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>是细胞潜在特征矩阵H的转置；/>是哈达玛积；/>是系数超参数。
式(5)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>表示细胞潜在特征矩阵Q的转置；/>是哈达玛积；/>是系数超参数。
具体到本实施例中，将原始基因归一化矩阵利用式(4)和式(5)所述的非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q1和基因潜在特征矩阵H1，其中细胞潜在特征矩阵Q1中每一列均表示的是细胞的潜在特征，基因潜在特征矩阵H1中每一列均表示的是变异特征基因的潜在特征；将原始基因归一化矩阵/>利用式(4)和式(5)所述的非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q2和基因潜在特征矩阵H2，其中细胞潜在特征矩阵Q2中每一列均表示的是细胞的潜在特征，基因潜在特征矩阵H2中每一列均表示的是变异特征基因的潜在特征；将原始基因归一化矩阵/>利用式(4)和式(5)所述的非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q3和基因潜在特征矩阵H3，其中细胞潜在特征矩阵Q3中每一列均表示的是细胞的潜在特征，基因潜在特征矩阵H3中每一列均表示的是变异特征基因的潜在特征；
步骤S3-2：为了控制本申请所述单细胞转录组测序数据插补方法的复杂度以及减小过拟合的风险，本申请对细胞潜在特征矩阵Q和基因潜在特征矩阵H施加L1正则化约束，得到L1正则化Y，L1正则化Y是通过对细胞潜在特征矩阵Q和基因潜在特征矩阵H的F范数和乘以系数获得的。施加L1正则化约束的计算方式，如式(6)所示：
(6)
式(6)中，Q表示细胞潜在特征矩阵，H表示基因潜在特征矩阵，是F范数,/>是基因潜在特征矩阵H的转置，/>是系数超参数。
具体到本实施例中，利用式(6)计算细胞潜在特征矩阵Q1和基因潜在特征矩阵H1的F范数之和乘以系数时，λ1取值为0.2；利用式(6)计算细胞潜在特征矩阵Q2和基因潜在特征矩阵H2的F范数之和乘以系数/>时，λ1取值为2；利用式(6)计算细胞潜在特征矩阵Q3和基因潜在特征矩阵H3的F范数之和乘以系数/>时，λ1取值为0.2；
步骤S3-3：计算原始基因归一化矩阵中相邻两列细胞的余弦距离，构建相似性矩阵S，而后，根据相似性矩阵S计算相似性矩阵S的图拉普拉斯矩阵/>，而后，根据图拉普拉斯矩阵/>计算相似性矩阵S的拉普拉斯矩阵L，而后对细胞潜在特征矩阵Q应用图正则化约束计算拉普拉斯矩阵L中矩阵特征值总和M。本申请中对细胞潜在特征矩阵Q应用图正则化约束，获得拉普拉斯矩阵L中矩阵特征值总和M，能够有效学习潜在的细胞亚群的结构，有效捕捉细胞在潜在空间中的几何关系。其中，相似性矩阵S、图拉普拉斯矩阵/>、拉普拉斯矩阵L以及矩阵特征值总和M的计算方式，分别如式(7)至式(10)所示：
(7)
(8)
(9)
(10)
式(7)中，是原始基因归一化矩阵/>中相邻两列细胞之间的余弦距离，/>表示最接近细胞/>的细胞的集合，/>分别表示最接近细胞/>的细胞的集合；
式(8)中，是图拉普拉斯矩阵，D是相似性矩阵S的度矩阵，W是相似性矩阵S的邻接矩阵，本申请中邻接矩阵W=相似性矩阵S；
式(9)中，L是拉普拉斯矩阵，是图拉普拉斯矩阵，D是相似性矩阵S的度矩阵；
式(10)中，Tr是迹，表示细胞潜在特征矩阵Q的转置，L是拉普拉斯矩阵，Q表示细胞潜在特征矩阵，/>是系数超参数。
具体到本实施例中，计算原始基因归一化矩阵中相邻两列细胞之间的余弦距离，利用式(7)所述方式构建相似性矩阵S1，而后，利用式(8)计算相似性矩阵S1的图拉普拉斯矩阵/>，而后，利用式(9)计算相似性矩阵S1的拉普拉斯矩阵L1，而后，对细胞潜在特征矩阵Q1应用图正则化约束，利用式(10)获得拉普拉斯矩阵L1中矩阵特征值总和M1，来学习潜在的细胞亚群的结构，捕捉细胞在潜在空间中的几何关系；其中，在利用式(10)计算拉普拉斯矩阵L1中矩阵特征值总和M1时，λ2取值为0.4；
计算原始基因归一化矩阵中相邻两列细胞之间的余弦距离，利用式(7)所述方式构建相似性矩阵S2，而后，利用式(8)计算相似性矩阵S2的图拉普拉斯矩阵/>，而后，利用式(9)计算相似性矩阵S2的拉普拉斯矩阵L2，而后，对细胞潜在特征矩阵Q2应用图正则化约束，利用式(10)获得拉普拉斯矩阵L2中矩阵特征值总和M2，来学习潜在的细胞亚群的结构，捕捉细胞在潜在空间中的几何关系；其中，在利用式(10)计算拉普拉斯矩阵L2中矩阵特征值总和M2时，λ2取值为4；
计算原始基因归一化矩阵中相邻两列细胞之间的余弦距离，利用式(7)所述方式构建相似性矩阵S3，而后，利用式(8)计算相似性矩阵S3的图拉普拉斯矩阵/>，而后，利用式(9)计算相似性矩阵S3的拉普拉斯矩阵L3，而后，对细胞潜在特征矩阵Q3应用图正则化约束，利用式(10)获得拉普拉斯矩阵L3中矩阵特征值总和M3，来学习潜在的细胞亚群的结构，捕捉细胞在潜在空间中的几何关系；其中，在利用式(10)计算拉普拉斯矩阵L3中矩阵特征值总和M3时，λ2取值为0.4；
步骤S3-4：以步骤S3-1得到的细胞潜在特征矩阵Q和基因潜在特征矩阵H以及步骤S3-3得到的相似性矩阵S为基础，使用最小二乘法迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，在迭代交替更新的过程中，同时更新细胞潜在特征矩阵Q的相似性矩阵S来更新对更新后的细胞潜在特征矩阵Q的图正则化约束，并利用跟新后的相似性矩阵S来更新图拉普拉斯矩阵以及更新细胞潜在特征矩阵Q，而后利用跟新后的图拉普拉斯矩阵/>来更新拉普拉斯矩阵L，而后利用跟新后的拉普拉斯矩阵L来更新拉普拉斯矩阵L中矩阵特征值总和M；而后利用更新后的矩阵特征值总和M更新近似矩阵N与原始基因归一化矩阵/>之间的F范数误差E，当F范数误差E小于10-5时，停止迭代，得到最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H；并将最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H相乘，得到最终的预测矩阵N，该最终的预测矩阵N与原始矩阵相近的预测矩阵N；其中，利用相似性矩阵S更新后的细胞潜在特征矩阵Q作为下一个迭代过程的基础；近似矩阵N的计算方式，如式(11)所示；本申请中F范数误差E表示的是近似矩阵N与原始基因归一化矩阵/>之间的F范数误差，F范数误差E的计算方式，如式(12)所示；而更新细胞潜在特征矩阵Q、更新基因潜在特征矩阵H以及更新后的细胞潜在特征矩阵Q的相似性矩阵S的计算方式，分别如式(13)、式(14)和式(15)所示：
(11)
(12)
(13)
(14)
(15)
式(11)中，Q表示细胞潜在特征矩阵，HT是基因潜在特征矩阵H的转置，N是预测矩阵；
式(12)中，P表示构建的系数矩阵，是哈达玛积，/>表示原始基因归一化矩阵，Q表示细胞潜在特征矩阵，HT是基因潜在特征矩阵H的转置，/>是F范数，Y表示L1正则化；M表示拉普拉斯矩阵L中矩阵特征值的总和。
式(13)中，P表示构建的系数矩阵，是哈达玛积，/>表示原始基因归一化矩阵，H表示基因潜在特征矩阵，D是相似性矩阵S的度矩阵，S表示相似性矩阵，Q表示细胞潜在特征矩阵，HT是基因潜在特征矩阵H的转置，/>为系数超参数，式(13)中的/>与式(4)至式(6)中的均相同，/>是系数超参数，式(13)中的/>与式(10)中的/>相同。
式(14)中，H表示基因潜在特征矩阵，是哈达玛积，/>是基因潜在特征矩阵Q的转置，P表示构建的系数矩阵，/>表示原始基因归一化矩阵，/>是基因潜在特征矩阵H的转置，Q表示细胞潜在特征矩阵，/>为系数超参数，式(14)中的/>与式(4)至式(6)中的/>均相同。
式(15)中，代表更新后的细胞潜在特征矩阵Q中相邻两列细胞的潜在特征之间的余弦距离，/>表示与细胞的潜在特征/>最接近的细胞的潜在特征的集合，/>表示与细胞的潜在特征/>最接近的细胞的潜在特征的集合。
具体到本实施例而言，以步骤S3-1得到的细胞潜在特征矩阵Q1和基因潜在特征矩阵H1以及步骤S3-3得到的相似性矩阵S1为基础，使用如式(12)至式(15)所述的最小二乘法迭代更新细胞潜在特征矩阵Q1和基因潜在特征矩阵H1，停止迭代后，得到最终的细胞潜在特征矩阵Q1和基因潜在特征矩阵H1；并将最终的细胞潜在特征矩阵Q1和基因潜在特征矩阵H1相乘，得到最终的预测矩阵N1，该最终的预测矩阵N1与原始矩阵相近的预测矩阵N1；
预测矩阵N2和预测矩阵N3的获取方式与预测矩阵N1的获取方式相同。
步骤S4：根据系数矩阵P、原始基因归一化矩阵以及与原始矩阵相近的的预测矩阵N通过计算得到插补后的矩阵/>，插补后的矩阵/>的计算方式，如式(16)所示；该步骤S4能够只针对由于技术原因所产生的技术零值进行插补；
(16)
式(16)中，是插补后的矩阵，P表示系数矩阵，/>是哈达玛积，N表示预测矩阵，/>表示原始基因归一化矩阵。
具体到本实施例中，根据系数矩阵P1以及步骤S3-4得到的预测矩阵N1通过式(16)计算得到插补后的矩阵；插补后的矩阵/>和插补后的矩阵/>的获取方式与插补后的矩阵获取方式相同。
为了对比本申请所述的单细胞转录组测序数据插补方法的优异性，本申请特地利用本申请所收集的公开数据集对ALRA插补方法(出自论文《Zero-preservingimputationof single-cell RNA-seq data》)、DrImpute插补方法（出自论文《“DrImpute： imputingdropout events in single cell RNA sequencing data》）、DeepImpute插补方法（出自论文《“DeepImpute： an accurate, fast, and scalable deep neural network method toimpute single-cell RNA-seq data》）、MAGIC插补方法（出自论文《Recovering geneinteractions from single-cell data using datadiffusion》）、SAVER插补方法（出自论文《SAVER：gene expression recovery for single-cell RNA sequencing》）、SCDD插补方法（出自论文《SCDD： a novel single-cell RNA-seq imputation method withdiffusion and denoising》）、scImpute插补方法（出自论文《An accurate androbustimputation method scImpute for single-cell RNA-seq data》）、scVI插补方法（出自论文《Interpretable dimensionality reduction of single cell transcriptomedata with deep generative models》）、WEDGE插补方法（出自论文《WEDGE：imputation ofgene expression values from single-cell RNA-seqdatasets using biased matrixdecomposition》）这九种现有的插补方法以及本申请所述的单细胞转录组测序数据插补方法分别对所收集的公开数据集（即li数据集和liver数据集）的进行插补处理，而后，再对插补后的单细胞转录组测序数据使用Seurat进行聚类后测试以分析利用本申请所述的单细胞转录组测序数据插补方法插补得到的单细胞转录组测序数据相对于上述其他九种现有的插补方法插补得到的单细胞转录组测序数据对下游细胞聚类任务的提升效果，测试结果，如表2、表3以及图2和图3所示。
表2
ARI NMI
Raw 0.8850048 0.9023147
Ours 0.9927305 0.9918196
ALRA插补方法 0.9822363 0.9822784
DeepImpute插补方法 0.8958311 0.9068955
DrImpute插补方法 0.9922798 0.9906099
MAGIC插补方法 0.7189163 0.7662982
SAVER插补方法 0.8678724 0.8957816
SCDD插补方法 0.9100630 0.9257576
scImpute插补方法 0.7948199 0.8497699
scVI插补方法 0.8297630 0.9016815
WEDGE插补方法 0.7383484 0.7899693
表2表示本发明以及上述九种现有的插补方法对单细胞转录组测序公开数据集li数据集进行插补之后再进行细胞聚类的得分。表2中，Raw 表示没有执行插补的原始数据矩阵进行细胞聚类的得分，Ours表示本发明所述插补方法对单细胞转录组测序公开数据集li数据集进行插补之后再进行细胞聚类的得分。
ARI是调整兰德指数，是一种综合考虑聚类结果相似性和随机性的度量方式，ARI的值越大，表示聚类结果越相似，与随机分类或者完全不同的情况差距越大，ARI对簇的编号敏感，适用于有真实标签的情况下的评估。NMI是归一化互信息，是一种用于度量两个聚类结果之间相似性的指标，它衡量了两个聚类结果之间的相互信息量，NMI数值越高表示聚类的精确性越好，NMI对簇的编号不敏感，因此不依赖于簇的标签，它适用于评估在没有真实标签的情况下的聚类的性能。本申请基于ARI（调整兰德指数）和NMI（归一化互信息）这两个指标评价本发明所述的单细胞转录组测序数据插补方法以及上述九种现有的插补方法在li数据集上插补之后的细胞聚类效果。
从表2中可以看出，本申请所述的单细胞转录组测序数据插补方法对li数据集进行插补之后进行细胞聚类得到的ARI和NMI两个评价指标，相对于DrImpute插补方法（上述九种插补方法中ARI和NMI两个评价指标最高的）对li数据集进行插补之后进行细胞聚类得到的ARI和NMI两个评价指标分别提高了0.045%和0.122%。这说明本申请所述的单细胞转录组测序数据插补方法对li数据集进行插补处理后得到的插补后的单细胞转录组测序数据在下游聚类任务方面表现更为优异。
从图2中还可以看出：本申请所述的单细胞转录组测序数据插补方法对li数据集进行插补处理，而后，再对插补后的单细胞转录组测序数据使用Seurat进行聚类后得到的UMAP降维的可视化散点图来看优于上述九种现有的插补方法插补后的单细胞转录组测序数据使用Seurat进行聚类后得到的UMAP降维的可视化散点图，具体来说就是，本发明所述的插补方法插补后的单细胞转录组测序数据使用Seurat进行聚类后得到的UMAP降维的可视化散点图中细胞轮廓更为清晰，细胞亚群内部更为紧凑，细胞亚群之间的距离更为明显，在插补的同时更好地保留了原始数据，更符合细胞亚群的分布。
表3表示本发明以及上述九种现有的插补方法对单细胞转录组测序公开数据集liver数据集进行插补之后再进行细胞聚类的得分。表3中，Raw 表示没有执行插补的原始数据矩阵进行细胞聚类的得分，Ours表示本发明所述插补方法对单细胞转录组测序公开数据集li数据集进行插补之后再进行细胞聚类的得分；
表3
ARI NMI
Raw 0.6790377 0.7555714
Ours 0.7835690 0.8123501
ALRA插补方法 0.7141175 0.7838701
DeepImpute插补方法 0.6683665 0.7495774
DrImpute插补方法 0.6790224 0.7541622
MAGIC插补方法 0.6237190 0.7282580
SAVER插补方法 0.6792750 0.7545047
SCDD插补方法 0.6240431 0.7130164
scImpute插补方法 0.7007967 0.7640005
scVI插补方法 0.6811972 0.7604257
WEDGE插补方法 0.7029150 0.7778239
从表3中可以看出：本申请所述的单细胞转录组测序数据插补方法对liver数据集进行插补后细胞聚类得到的ARI和NMI两个评价指标，相对于ALRA插补方法（上述九种插补方法中ARI和NMI两个评价指标最高的）对liver数据集进行插补后细胞聚类得到的ARI和NMI两个评价指标分别提高了9.726%和3.363%。这说明本申请所述的单细胞转录组测序数据插补方法对liver数据集进行插补处理后得到的插补后的单细胞转录组测序数据在下游聚类任务方面表现更为优异。
从图3中还可以看出：本申请所述的单细胞转录组测序数据插补方法对liver数据集进行插补处理后，使用Seurat进行聚类后测试得到的UMAP降维的可视化散点图要优于上述九种现有的插补方法对liver数据集进行插补处理后，使用Seurat进行聚类后测试得到的UMAP降维的可视化散点图，具体来说，本发明所述的插补方法对liver数据集进行插补处理后，使用Seurat进行聚类后测试得到的UMAP降维的可视化散点图生成的UMAP降维的可视化散点图细胞轮廓更为清晰，细胞亚群之间的距离更为明显，在插补的同时更好地保留了原始数据，更符合细胞亚群的分布。
此外，为了对比本申请所述的单细胞转录组测序数据插补方法对下游重建细胞轨迹任务的提升，本申请特地利用上述九种插补方法对Deng数据集进行插补，并使用Monocle3进行轨迹推断测试，测试结果，如表4和图4所示：
表4
KOR POS
Raw 0.31 0.19
scCAN 0.79 0.93
ALRA插补方法 0.11 0.04
DeepImpute插补方法 0.67 0.87
DrImpute插补方法 0.46 0.66
MAGIC插补方法 0.73 0.91
SAVER插补方法 0.26 0.02
SCDD插补方法 0.12 0.17
scImpute插补方法 0.41 0.36
scVI插补方法 0.8 0.14
WEDGE插补方法 0.67 0.87
表4表示本发明以及上述九种现有插补方法在单细胞转录组测序公开数据集Deng数据集上进行插补之后再进行重建细胞轨迹的得分。表4中，Raw 表示没有执行插补的原始数据矩阵进行重建细胞轨迹的得分，scCAN表示本发明所述插补方法对单细胞转录组测序公开数据集Deng数据集进行插补之后再进行重建细胞轨迹的得分；
表4中，POS是伪时间排序，是用于评估细胞运动轨迹中伪时间排序准确性的指标，POS通过分析基因表达的动态变化来估计细胞的伪时间顺序。POS衡量的是估计的伪时间与实际细胞发育或状态之间的一致性或相关性。POS分数值越高表示估计的伪时间与实际发育或状态的一致性更好。KOR是肯德尔等级相关，在细胞运动轨迹的上下文中，用于比较两个估计的运动轨迹的相似性，还可用来评估不同的细胞的运动轨迹之间的一致性。KOR分数值越高，表明如果两个轨迹的相对顺序越相似。
本申请基于POS（伪时间排序）和KOR（肯德尔等级相关）这两个指标评价本发明所述的单细胞转录组测序数据插补方法以及上述九种现有插补方法对Deng数据集进行插补之后再进行重建细胞轨迹的效果。
从表4可知，本申请所述的单细胞转录组测序数据插补方法对Deng数据集进行插补之后再进行重建细胞轨迹得到的KOR和POS两个评价指标分别达到了0.79和0.93，相对于MAGIC插补方法（上述九种插补方法中KOR和POS两个评价指标最高的）分别提高了8.219%和2.198%。这说明本申请所述的单细胞转录组测序数据插补方法对Deng数据集进行插补处理后得到的插补后的单细胞转录组测序数据在下游轨迹推测任务方面表现更为优异。
此外，从图4中还可以看出，本发明所述的单细胞转录组测序数据插补方法对Deng数据集进行插补之后再进行重建细胞轨迹生成的细胞轨迹最为清晰可靠，最符合真实的细胞发育轨迹，插补后的总体的细胞发育轨迹连续且不间断。
综上可知，本发明所述的单细胞转录组测序数据插补方法相较于上述九种现有插补方法，本发明对li数据集和liver数据集进行插补处理后得到的插补后的单细胞转录组测序数据在下游的细胞聚类任务能够取得更优的效果；而且，本发明所述插补方法对Deng数据集进行插补处理后得到的插补后的单细胞转录组测序数据在下游的重建细胞轨迹的任务中也能取得更优的效果。这也说明了本申请所述的单细胞转录组测序数据插补方法在进行插补时，能够同时在两个下游任务中都表现出较佳的插补效果。

Claims (10)

1.一种单细胞转录组测序数据插补方法，其特征在于：包括以下步骤：
S1：筛选变异特征基因，根据筛选出的变异特征基因及与其相对应的细胞类型构建原始基因表达矩阵V；而后，对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵；
S2：识别原始基因表达矩阵V中潜在的缺失值，构建用于仅对由测序技术原因所产生的技术零值的位置进行插补的系数矩阵P；
S3：利用非负矩阵分解算法分解原始基因归一化矩阵，得到细胞潜在特征矩阵Q 和基因潜在特征矩阵H，而后，对细胞潜在特征矩阵Q和基因潜在特征矩阵H进行L1正则化约束，而后，对细胞潜在特征矩阵Q应用图正则化约束捕捉细胞在潜在空间中的几何关系，而后，迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，直至收敛，结束迭代更新过程，得到最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H，并利用最终的细胞潜在特征矩阵Q和基因潜在特征矩阵H得到与原始矩阵相近的预测矩阵N；
S4：根据系数矩阵P、原始基因归一化矩阵以及与原始矩阵相近的预测矩阵N计算得到插补后的矩阵/>。
2.根据权利要求1所述的单细胞转录组测序数据插补方法，其特征在于：所述步骤S1具体为：
S1-1、从目标数据集中获取细胞类型和与其相对应的基因数据，而后，使用Seurat程序包对目标数据集中的基因数据进行裁剪，并保留基因数据裁剪处理后得到的一定数量的变异特征基因；
S1-2、利用步骤S1-1保留的变异特征基因，并结合其相适应的细胞类型的名称构成原始基因表达矩阵V；
S1-3、对原始基因表达矩阵V进行归一化处理，得到原始基因归一化矩阵。
3.根据权利要求2所述的单细胞转录组测序数据插补方法，其特征在于：步骤S1-1中，当目标数据集由多个不同数据集组成时，不同数据集基因数据裁剪处理后，针对不同数据集均保留一致数量的变异特征基因。
4.根据权利要求1所述的单细胞转录组测序数据插补方法，其特征在于：步骤S2具体包括如下步骤：
S2-1、使用谱聚类将原始基因表达矩阵V中的细胞进行大致聚类，分别得到k个细胞亚群；
S2-2、计算原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值/>和方差/>，而后，基于每一个细胞亚群中的变异特征基因的零表达率/>、平均表达值/>和方差/>计算得到用于区分变异特征基因元素为生物零和技术零的置信度β；当置信度β小于阈值t时，则将该变异特征基因的元素值为零的值，视为因为生物变异所产生的零值，在与原始基因表达矩阵V相对应的系数矩阵P中的对应位置的元素值设置为1；
S2-3、遍历原始基因表达矩阵V中每条变异特征基因的元素值，如果元素值大于零就将与该元素值在原始基因表达矩阵V中位置相对应的系数矩阵P中的位置处的元素值设置为1；而后，再将原始基因表达矩阵V中剩余的元素值为零值的位置，在系数矩阵P中的对应位置处的元素值设置为0，得到完整的系数矩阵P。
5.根据权利要求4所述的单细胞转录组测序数据插补方法，其特征在于：步骤S2-2中，利用原始基因表达矩阵V中的每一个细胞亚群中的变异特征基因的零表达率、平均表达值和方差/>计算用于区分变异特征基因元素为生物零和技术零的置信度β，如式(2)所示：
(2)
式(2)中，表示在第k个细胞亚群中第i个基因的置信度，/>表示在第k个细胞亚群中第i个基因的零表达率，/>表示在第k个细胞亚群中第i个基因的表达的平均值，/>表示在第k个细胞亚群中第i个基因的表达的方差。
6.根据权利要求1所述的单细胞转录组测序数据插补方法，其特征在于：步骤S3具体为：
S3-1、将原始基因归一化矩阵利用非负矩阵分解算法随机分解为两个基因表达矩阵，分别为细胞潜在特征矩阵Q和基因潜在特征矩阵H；
S3-2、对细胞潜在特征矩阵Q和基因潜在特征矩阵H施加L1正则化约束，得到L1正则化Y；
S3-3、计算原始基因归一化矩阵中相邻两列细胞的余弦距离，构建相似性矩阵S，并根据相似性矩阵S计算相似性矩阵S的图拉普拉斯矩阵/>，根据图拉普拉斯矩阵/>计算相似性矩阵S的拉普拉斯矩阵L，而后对细胞潜在特征矩阵Q应用图正则化约束计算拉普拉斯矩阵L中矩阵特征值总和M；
S3-4、以步骤S3-1得到的细胞潜在特征矩阵Q和基因潜在特征矩阵H以及步骤S3-3得到的相似性矩阵S为基础，使用最小二乘法迭代更新细胞潜在特征矩阵Q和基因潜在特征矩阵H，在迭代交替更新的过程中，同时计算更新后的细胞潜在特征矩阵Q的相似性矩阵S来更新对更新后的细胞潜在特征矩阵Q图正则化约束，而后，利用跟新后的相似性矩阵S来更新图拉普拉斯矩阵以及更新细胞潜在特征矩阵Q，而后利用跟新后的图拉普拉斯矩阵/>来更新拉普拉斯矩阵L，而后利用跟新后的拉普拉斯矩阵L来更新拉普拉斯矩阵L中矩阵特征值总和M；而后利用更新后的矩阵特征值总和M更新近似矩阵N与原始基因归一化矩阵/>之间的F范数误差E，直至F范数误差E收敛，停止迭代，得到最终的细胞潜在特征矩阵Q、最终的基因潜在特征矩阵H以及与原始矩阵相近的预测矩阵N；所述近似矩阵N是通过细胞潜在特征矩阵Q和基因潜在特征矩阵H相乘得到的。
7.根据权利要求6所述的单细胞转录组测序数据插补方法，其特征在于：步骤S3-1中，将原始基因归一化矩阵利用非负矩阵分解算法随机分解为细胞潜在特征矩阵Q和基因潜在特征矩阵H的计算方式，如式(4)和式(5)所示：
(4)
(5)
式(4)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>表示细胞潜在特征矩阵H的转置；/>是哈达玛积；/>是系数超参数；
式(5)中，是表示原始基因归一化矩阵；P表示构建的系数矩阵；Q表示细胞潜在特征矩阵；H表示基因潜在特征矩阵；/>是细胞潜在特征矩阵Q的转置；/>是哈达玛积；/>是系数超参数。
8.根据权利要求6所述的单细胞转录组测序数据插补方法，其特征在于：步骤S3-2中，L1正则化Y是通过对细胞潜在特征矩阵Q和基因潜在特征矩阵H的F范数之和乘以系数获得的。
9.根据权利要求6所述的单细胞转录组测序数据插补方法，其特征在于：步骤S3-4中，F范数误差E收敛的条件是F范数误差E小于10-5。
10.根据权利要求1所述的单细胞转录组测序数据插补方法，其特征在于：步骤S4中，插补后的矩阵的计算方式，如式(16)所示：
(16)
式(16)中，是插补后的矩阵，P表示系数矩阵，/>是哈达玛积，N表示预测矩阵；/>表示原始基因归一化矩阵。
CN202311763792.3A 2023-12-21 2023-12-21 一种单细胞转录组测序数据插补方法 Active CN117437973B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311763792.3A CN117437973B (zh) 2023-12-21 2023-12-21 一种单细胞转录组测序数据插补方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311763792.3A CN117437973B (zh) 2023-12-21 2023-12-21 一种单细胞转录组测序数据插补方法

Publications (2)

Publication Number Publication Date
CN117437973A true CN117437973A (zh) 2024-01-23
CN117437973B CN117437973B (zh) 2024-03-08

Family

ID=89546491

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311763792.3A Active CN117437973B (zh) 2023-12-21 2023-12-21 一种单细胞转录组测序数据插补方法

Country Status (1)

Country Link
CN (1) CN117437973B (zh)

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109906276A (zh) * 2016-11-07 2019-06-18 格里尔公司 用于检测早期癌症中体细胞突变特征的识别方法
US20190316209A1 (en) * 2018-04-13 2019-10-17 Grail, Inc. Multi-Assay Prediction Model for Cancer Detection
CN110797089A (zh) * 2019-10-30 2020-02-14 华东交通大学 一种基于单细胞rna测序数据识别细胞类型的方法
CA3145700A1 (en) * 2019-07-25 2021-01-28 Jeanne F. Loring Methods of identifying dopaminergic neurons and progenitor cells
WO2021109419A1 (zh) * 2019-12-05 2021-06-10 东南大学 大规模mimo波束域鲁棒预编码传输方法与系统
CN113537358A (zh) * 2021-07-19 2021-10-22 华南理工大学 一种基于多组学数据集的癌症亚型识别方法及系统
US20220122250A1 (en) * 2020-10-19 2022-04-21 Northwestern University Brain feature prediction using geometric deep learning on graph representations of medical image data
CN116779034A (zh) * 2023-05-15 2023-09-19 湖南科技大学 miRNA与疾病关联预测方法、设备及存储介质
CN117133359A (zh) * 2022-05-19 2023-11-28 天津科技大学 一种基于拉普拉斯正则化的单细胞RNA-seq数据的插补方法
CN117253550A (zh) * 2023-09-08 2023-12-19 湖南工业大学 一种空间转录组数据聚类方法

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109906276A (zh) * 2016-11-07 2019-06-18 格里尔公司 用于检测早期癌症中体细胞突变特征的识别方法
US20190316209A1 (en) * 2018-04-13 2019-10-17 Grail, Inc. Multi-Assay Prediction Model for Cancer Detection
CA3145700A1 (en) * 2019-07-25 2021-01-28 Jeanne F. Loring Methods of identifying dopaminergic neurons and progenitor cells
CN110797089A (zh) * 2019-10-30 2020-02-14 华东交通大学 一种基于单细胞rna测序数据识别细胞类型的方法
WO2021109419A1 (zh) * 2019-12-05 2021-06-10 东南大学 大规模mimo波束域鲁棒预编码传输方法与系统
US20220122250A1 (en) * 2020-10-19 2022-04-21 Northwestern University Brain feature prediction using geometric deep learning on graph representations of medical image data
CN113537358A (zh) * 2021-07-19 2021-10-22 华南理工大学 一种基于多组学数据集的癌症亚型识别方法及系统
CN117133359A (zh) * 2022-05-19 2023-11-28 天津科技大学 一种基于拉普拉斯正则化的单细胞RNA-seq数据的插补方法
CN116779034A (zh) * 2023-05-15 2023-09-19 湖南科技大学 miRNA与疾病关联预测方法、设备及存储介质
CN117253550A (zh) * 2023-09-08 2023-12-19 湖南工业大学 一种空间转录组数据聚类方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
JUNCHEN YANG ET.AL.: ""Multi-model Differentiable Unsupervised Feature Selection"", 《ARXIV》, 16 March 2023 (2023-03-16) *
卢志强;莫晓晖;吴松松;: "基于离散指示矩阵优化的多视图子空间聚类算法", 金陵科技学院学报, vol. 35, no. 01, 30 March 2019 (2019-03-30), pages 11 - 15 *
王珊珊: ""低秩非负矩阵分解算法在聚类中的研究"", 《万方在线出版》, 1 September 2023 (2023-09-01) *
谢豪杰: ""基于低秩张量补全与张量分解的交通数据填补方法"", 《万方在线出版》, 17 October 2023 (2023-10-17) *
陆聪海: ""稀疏低秩表示模型的研究及在癌症测序数据中的应用"", 《中国优秀硕士学位论文全文数据库 医药卫生科技辑》, vol. 2021, no. 01, 15 January 2021 (2021-01-15), pages 072 - 80 *

Also Published As

Publication number Publication date
CN117437973B (zh) 2024-03-08

Similar Documents

Publication Publication Date Title
CN110827921B (zh) 一种单细胞聚类方法、装置、电子设备及存储介质
Zhou et al. Subspace segmentation-based robust multiple kernel clustering
Yu et al. Self-paced learning for k-means clustering algorithm
Maulik et al. Simulated annealing based automatic fuzzy clustering combined with ANN classification for analyzing microarray data
CN111785329A (zh) 基于对抗自动编码器的单细胞rna测序聚类方法
Data Machine learning
WO2013067461A2 (en) Identifying associations in data
CN106886793B (zh) 基于判别信息和流形信息的高光谱图像波段选择方法
CN115240772A (zh) 一种基于图神经网络的解析单细胞多组学中活性通路的方法
CN108764276A (zh) 一种鲁棒自动加权多特征聚类方法
CN110348287A (zh) 一种基于字典和样本相似图的无监督特征选择方法和装置
CN116580848A (zh) 一种基于多头注意力机制的分析癌症多组学数据方法
Vengatesan et al. The performance analysis of microarray data using occurrence clustering
Ding et al. Dance: A deep learning library and benchmark for single-cell analysis
CN117437973B (zh) 一种单细胞转录组测序数据插补方法
Vignes et al. Gene clustering via integrated Markov models combining individual and pairwise features
CN111863135A (zh) 一种假阳性结构变异过滤方法、存储介质及计算设备
CN116705151A (zh) 一种空间转录组数据的降维方法及系统
Gan et al. DSAE-Impute: Learning discriminative stacked autoencoders for imputing single-cell rna-seq data
CN115881232A (zh) 一种基于图神经网络和特征融合的scRNA-seq细胞类型注释方法
CN115661498A (zh) 一种自优化单细胞聚类方法
Yang et al. Robust landmark graph-based clustering for high-dimensional data
CN104573727A (zh) 一种手写体数字图像降维方法
CN112768001A (zh) 一种基于流形学习和主曲线的单细胞轨迹推断方法
CN109215741B (zh) 基于双超图正则化的肿瘤基因表达谱数据双聚类方法

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
GR01 Patent grant
GR01 Patent grant