CN116153404B - 一种单细胞ATAC-seq数据分析方法 - Google Patents

一种单细胞ATAC-seq数据分析方法 Download PDF

Info

Publication number
CN116153404B
CN116153404B CN202310182496.8A CN202310182496A CN116153404B CN 116153404 B CN116153404 B CN 116153404B CN 202310182496 A CN202310182496 A CN 202310182496A CN 116153404 B CN116153404 B CN 116153404B
Authority
CN
China
Prior art keywords
cell
dna
matrix
atac
characteristic peak
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.)
Active
Application number
CN202310182496.8A
Other languages
English (en)
Other versions
CN116153404A (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.)
Chengdu University of Information Technology
Original Assignee
Chengdu University of Information 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 Chengdu University of Information Technology filed Critical Chengdu University of Information Technology
Priority to CN202310182496.8A priority Critical patent/CN116153404B/zh
Publication of CN116153404A publication Critical patent/CN116153404A/zh
Application granted granted Critical
Publication of CN116153404B publication Critical patent/CN116153404B/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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • 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
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Medical Informatics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Molecular Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Biotechnology (AREA)
  • Databases & Information Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Epidemiology (AREA)
  • Bioethics (AREA)
  • Public Health (AREA)
  • Genetics & Genomics (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种单细胞ATAC‑seq数据分析方法，通过提取单细胞分辨率的染色质可达性特征峰序列中转录因子‑DNA结合基元的所属种类、相对位置、长距离依赖关系等众多转录调控语法规则，从而更全面地表示单个细胞的功能状态和高阶特征。此外，本发明方法利用获取的转录调控语法规则、细胞功能状态和高阶特征，一站式地实现染色质可达性预测、细胞类型注释、染色质可达性图谱降噪、转录因子活性推断等一系列下游分析任务。

Description

一种单细胞ATAC-seq数据分析方法
技术领域
本发明属于生物信息分析技术领域，具体涉及一种单细胞ATAC-seq数据分析方法。
背景技术
单细胞染色质可达性靶向测序(single-cell Assay for TargetingAccessible-Chromatin sequencing，scATAC-seq)方法的兴起为刻画单细胞分辨率的染色质可达性图谱奠定了坚实的基础，其已经成为揭示基因转录过程中细胞特异性调控机制的重要手段。近年来，随着单细胞ATAC-seq数据规模的快速增长，即scATAC-seq方法检测到的染色质可达性特征峰的规模的快速增长，生物信息学家们开始着手致力于利用单细胞ATAC-seq数据实现染色质可达性预测、细胞类型注释、染色质可达性图谱降噪、转录因子活性推断等下游分析任务。然而，单细胞ATAC-seq数据固有的高维性、二值性、稀疏性为计算分析带来了巨大挑战。因此，设计精准、鲁棒、高效的单细胞ATAC-seq数据分析方法已然成为生物信息学领域中亟待解决的关键问题。
截至目前为止，生物信息学家们已经提出了一系列单细胞ATAC-seq数据分析方法，它们的基本思想是利用数理统计、机器学习等手段从海量单细胞ATAC-seq数据中提取每个细胞的高阶特征，并将所提取的特征应用于下游分析任务中。例如，威库鲁文脑与疾病研究中心的斯坦因·艾茨教授领衔设计的cisTopic算法、清华大学张强峰教授领衔设计的SCALE算法分别采用隐式迪利克雷分配模型、高斯变分自编码器模型等生成式模型将单个细胞映射到低维的、连续的、广义的特征空间中，有效地克服了单细胞ATAC-seq数据的高维性、二值性、稀疏性，成功地实现了细胞类型注释、染色质可达性图谱降噪等下游分析任务。然而，cisTopic和SCALE算法仅通过染色体坐标表示染色质可达性特征峰，忽略了DNA序列中潜藏的转录调控语法规则，导致了单细胞特征表示的不准性。斯坦福大学医学院的威廉·格林利夫教授领衔设计的chromVAR算法、美国麻省理工学院-哈佛大学博德研究所的阿维夫·雷格夫教授领衔设计的BROCKMAN算法分别使用转录因子-DNA结合基元的数量、染色质可达性特征峰中k-mer碱基字符串等生物学先验知识刻画不同细胞的特征，相比于cisTopic和SCALE而言进一步提高了转录因子活性推断等下游分析任务的可用性、精准性、可解释性。“谷歌”公司旗下的计算生物学家大卫·凯利教授领衔设计的scBasset算法运用卷积神经网络从已知染色质可达性特征峰的一级序列中提取转录因子-DNA结合基元的基元信息，开创性地实现了单细胞染色质可达性的预测任务。
但上述chromVAR、BROCKMAN、scBasset等广泛使用的单细胞ATAC-seq数据分析方法依然存在一些不足之处。第一点，大量研究表明启动子和增强子等具有生物功能的DNA短序列能够通过发生相互作用控制各类细胞的状态，然而现有的单细胞ATAC-seq数据分析方法难以有效地捕捉启动子和增强子等转录因子-DNA结合基元的相对位置和长距离依赖关系，进而无法精准地刻画启动子和增强子之间的相互作用。第二点，现有的单细胞ATAC-seq数据分析方法呈现各司其职的状态，难以将染色质可达性预测、细胞类型注释、染色质可达性图谱降噪、转录因子活性推断等重要的下游分析任务集成到一个统一的框架中，导致了方法的通用性有一定局限。
发明内容
针对现有技术中的上述不足，本发明提供的单细胞ATAC-seq数据分析方法解决了现有分析方法中，无法精准地刻画启动子和增强子之间的相互作用，以及方法通用性存在局限的问题，本发明提供一种精准、统一、高效的基于DNA语言模型的单细胞ATAC-seq数据计算分析方法。
为了达到上述发明目的，本发明采用的技术方案为：一种单细胞ATAC-seq数据分析方法，包括以下步骤：
S1、采集ATAC-seq特征峰的DNA一级序列，作为ATAC-seq数据集；
S2、通过基于ProbDep Transformer的DNA语言模型对DNA一级序列进行分析，预测DNA一级序列中各细胞的染色质可达性，并学习每个细胞的功能状态和高阶特征；
S3、根据学习的功能状态和高阶特征，进行细胞类型注释；
S4、基于预测的染色质可达性，进行ATAC-seq特征峰的染色质可达性图谱降噪；
S5、通过DNA语言模型，分析ATAC-seq特征峰中每个转录因子在各细胞中的活跃性；
S6、将染色质可达性预测结果、细胞类型注释、降噪的染色质可达性图谱以及转录因子的活跃性分析结果作为单细胞ATAC-seq数据分析结果。
进一步地，所述步骤S2具体为：
S21、将长度为L的DNA一级序列采用独热编码映射至维数为4×L的隐式特征空间中，并将其转换为基元编码矩阵；
对DNA一级序列采用绝对位置编码生成维数为pos×2i的位置编码矩阵，将基元编码矩阵和位置编码矩阵相加作为DNA语言模型的输入数据；其中，pos为当前转录因子-DNA结合基元在DNA一级序列中的位置下标，2i为当前转录因子-DNA结合基元的位置编码向量的长度；
S22、在DNA模型中，采用长距离依赖性测量评估查询与键之间的依赖性的方法对输入数据进行分析，获得每个键向量聚焦于排名最高的u个查询向量，进而得到DNA模型的输出数据；
S23、将DNA模型的输出数据作为DNA一级序列的高维语义编码，通过序列高阶编码器将其映射到低维空间中，获得DNA一级序列的高阶特征；
S24、根据获得的高阶特征，通过染色质可达性预测器预测当前DNA一级序列在各细胞的染色质可达性大小，并学习得到每个细胞的功能状态和高阶特征。
进一步地，所述步骤S22中，长距离依赖性测量的表达式为：
式中，为长距离依赖性测量操作，/>表示qi向量与全部键向量之间经过Log-Sum-Exp操作之后的结果，/>表示Log-Sum-Exp结果的算术平均值，qi为查询矩阵Q的第i行，K为键矩阵，In为求对数操作，l为键矩阵K中当前行的下标，Lk为键矩阵K中行的个数，/>为键矩阵K中的第k行的转置，d为键(key)矩阵K中列的个数；
DNA模型的输出数据表示为：
式中，为自注意力机制操作的输出，/>为自注意力机制操作，Q,K,V分别为查询矩阵，键矩阵，值矩阵，Softmax(·)为激活函数，/>为与矩阵Q大小相同的稀疏矩阵，且只包含了了长距离依赖性测量中排名最高的u个查询向量。
进一步地，所述步骤S23中，DNA一级序列的高阶特征为：
式中，ELU(·)为ELU激活函数，Wf和bf分别为序列高阶编码器的权重矩阵和截距向量，Conv1(·)为1维卷积操作。
进一步地，当前DNA一级序列在各细胞的染色质可达性大小y为：
式中，σ(·)为Sigmoid激活函数，Wp和bp分为染色质可达性预测器的权重矩阵和截距向量，其中，权重矩阵作为全部细胞的功能状态和高阶特征矩阵。
进一步地，所述步骤S3具体为：
S31、根据学习的功能状态和高阶特征构建k-近邻图；
其中，在构建的k-近邻图中节点为细胞，边为细胞之间的相关性；
S32、在构建的k-近邻图中，划分每个细胞的类型，并进一步在二维空间中可视化细胞类型注释结果。
进一步地，所述步骤S4具体为：
S41、获取原始ATAC-seq数据，即染色质可达性特征峰矩阵；
S42、通过混合概率模型搜索染色质可达性特征峰矩阵中由实验误差导致的零计数，进而计算特征峰i在细胞类型k的细胞m中的丢失率；
S43、对于给定的细胞m，根据其每个特征峰的丢失率，将全部特征峰划分为待修正的特征峰集合Am和无需修正的特征峰集合Bm；
S44、将通过DNA语言模型预测的特征峰i在每个细胞中的染色质可达性大小作为候选修正计数，对特征峰集合Am中的特征峰进行修正，进而获得降噪后的ATAC-seq特征峰的染色质可达性图谱。
进一步地，所述步骤S42中的混合概率模型包括伽马分布模型和正态分布模型，其中伽马分布模型表示实验误差，正态分布模型表示scATAC-seq特征峰的真实计数。
进一步地，所述步骤S42中，特征峰i在细胞类型k的细胞m中的丢失率dim为：
式中，表示参数的估计值，/>为特征峰i在细胞类型k中的整体丢失率，和/>为伽马分布的形状参数和尺度参数，/>和/>为正态分布的均值和标准差。
进一步地，所述步骤S5具体为：
S51、随机选择s条scATAC-seq特征峰序列进行二核苷酸打乱生成s条的背景序列；
S52、设定一个转录因子，将其与DNA结合基元插入到s条背景序列中心位置，生成s条合成序列；
S53、通过DNA模型分别预测s条背景序列和s条合成序列的染色质可达性；
S54、将单个细胞的背景序列和合成序列的平均染色质可达性预测差值作为当前转录因子在该单个细胞中的活跃性分析结果。
本发明的有益效果为：
1)本发明利用Prob Transformer有效捕捉了scATAC-seq特征峰序列中各类转录因子-DNA结合基元的所属种类、相对位置、长距离依赖关系，从而更精准、更全面地刻画了单个细胞的功能状态和高阶特征；
2)本发明方法能够将染色质可达性预测、细胞类型注释、染色质可达性图谱降噪、转录因子活性推断任务集成到了一个的框架中，显著地提高了单细胞ATAC-seq数据计算分析框架的通用性、统一性；
3)本发明方法支持在多GPU并行运算，可用于超大规模单细胞ATAC-seq数据的分析。
附图说明
图1为本发明提供的单细胞ATAC-seq数据分析方法流程图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
本发明实施例提供了一种单细胞ATAC-seq数据分析方法，如图1所示，包括以下步骤：
S1、采集ATAC-seq特征峰的DNA一级序列，作为ATAC-seq数据集；
S2、通过基于ProbDep Transformer的DNA语言模型对DNA一级序列进行分析，预测DNA一级序列中各细胞的染色质可达性，并学习每个细胞的功能状态和高阶特征；
S3、根据学习的功能状态和高阶特征，进行细胞类型注释；
S4、基于预测的染色质可达性，进行ATAC-seq特征峰的染色质可达性图谱降噪；
S5、通过DNA语言模型，分析ATAC-seq特征峰中每个转录因子在各细胞中的活跃性；
S6、将染色质可达性预测结果、细胞类型注释、降噪的染色质可达性图谱以及转录因子的活跃性分析结果作为单细胞ATAC-seq数据分析结果。
本发明实施例的步骤S2中，对于给定的ATAC-seq特征峰的DNA一级序列，利用基于ProbDep Transformer的DNA语言模型从中捕捉转录因子-DNA结合基元的所属种类、相对位置、长距离依赖关系，从而预测给定DNA序列在各细胞中的染色质可达性大小，并学习每个细胞的功能状态和高阶特征，人类基因组中全部细胞的功能状态和高阶特征可视为一个数值型矩阵，其中每一列表示一种细胞，每一行表示一种特征。
基于此，本实施例步骤S2中采用DNA语言模型预测染色质可达性，确定细胞功能状态和高阶特征的过程具体为：
S21、将长度为L的DNA一级序列采用独热编码映射至维数为4×L的隐式特征空间中，并将其转换为基元编码矩阵；
对DNA一级序列采用绝对位置编码生成维数为pos×2i的位置编码矩阵，将基元编码矩阵和位置编码矩阵相加作为DNA语言模型的输入数据；其中，pos为当前转录因子-DNA结合基元在DNA一级序列中的位置下标，2i为当前转录因子-DNA结合基元的位置编码向量的长度；
S22、在DNA模型中，采用长距离依赖性测量评估查询与键之间的依赖性的方法对输入数据进行分析，获得每个键向量聚焦于排名最高的u个查询向量，进而得到DNA模型的输出数据；
S23、将DNA模型的输出数据作为DNA一级序列的高维语义编码，通过序列高阶编码器将其映射到低维空间中，获得DNA一级序列的高阶特征；
S24、根据获得的高阶特征，通过染色质可达性预测器预测当前DNA一级序列在各细胞的染色质可达性大小，并学习得到每个细胞的功能状态和高阶特征。
本实施例的步骤S21中，对于长度为L的DNA一级序列，如果L小于人工设定的阈值，则通过一次卷积操作将独热编码矩阵转换为基元编码矩阵，其刻画了DNA序列中各类转录因子-DNA结合基元的存在性，即所含结合基元的种类；如果L大于人工设定的阈值，为降低基元编码矩阵的维数，则通过交替的卷积操作和池化操作生成基元编码矩阵。
本实施例的步骤S21中，位置编码矩阵刻画了各类转录因子-DNA结合基元的位置信息，位置编码矩阵表示为：
本实施例的步骤S21中，生成基元编码矩阵和位置编码矩阵后，直接将它们相加作为模型的输入
式中，u为基元编码矩阵，α为权重平衡因子。
传统Transformer的基本思想是采用自注意力机制计算各类转录因子-DNA结合基元之间的依赖性大小(以注意力权重的形式表示)，并根据依赖性大小通过加权求和的方式整合各转录因子-DNA结合基元的编码特征，从而进一步刻画结合基元之间的长距离依赖关系。定义qi、ki、vi分别为矩阵Q、K、V的第i行，其中Q、K、V分别表示查询(query)、键(key)、值(value)。给定查询向量qi，传统Transformer中其注意力函数定义为概率形式的核平滑器：
其中，和κ(qi,kj)选择非对称指数核/>然而事实上，长距离依赖关系只存在于诸如启动子和增强子等少部分转录因子-DNA结合基元之间，导致了自注意力机制的分布稀疏性，即少量查询(query)-键(key)的注意力权重较大，而其它大量查询-健的注意力权重接近于零。
因此，本发明实施例设计了一种ProbDepTransformer捕捉转录因子-DNA结合基元之间真实存在的长距离依赖关系，从而降低计算消耗、提高模型精度。在本实施例中，为区分查询(query)-键(key)的注意力权重是否影响注意力机制的计算，采用长距离依赖性测量定量评估查询(query)与键(key)之间的依赖性。从上述核平滑器公式中可看出，查询(query)i-全部键(key)的注意力为p(kj|qi)，输出是它与向量v的组合。如果p(kj|qi)接近均匀分布自注意力将变为计算向量v的平均值，其对于模型的推理是冗余的。自然的，分布p(kj|qi)和q(κj|qi)之间的相似性可以用来区分重要的查询(query)，继而找到有生物学意义的κ(qi,kj)。
因此，在本实施例的步骤S22中，给定查询(query)i，其长距离依赖性测量的表达式为：
式中，为长距离依赖性测量操作，/>表示qi向量与全部键向量之间经过Log-Sum-Exp操作之后的结果，/>表示Log-Sum-Exp结果的算术平均值，qi为查询(query)矩阵Q的第i行，K为键(key)矩阵，In为求对数操作，l为键(key)矩阵K中当前行的下标，Lk为键(key)矩阵K中行的个数，/>为键(key)矩阵K中的第k行的转置，d为键(key)矩阵K中列的个数；如果查询向量i的/>越大，它的注意力权重分布p(kj|qi)越多样化，且有极大概率产生交大的注意力权重，反之亦然，根据/>的结果，选择排名最高的u个查询向量作为矩阵/>其余元素均用零填充。
根据长距离依赖性测量的结果，自注意力机制的计算过程中每个键向量只聚焦于排名最高的u个查询向量，得到DNA模型的输出数据表示为：
式中，为自注意力机制操作的输出，/>为自注意力机制操作，Q,K,V分别为查询(query)矩阵，键(key)矩阵，值(value)矩阵，Softmax(·)为激活函数，/>为与矩阵Q大小相同的稀疏矩阵，且只包含了了长距离依赖性测量中排名最高的u个查询向量；与Transformer类似，ProbDepTransformer中同样添加了残差连接和前馈神经网络。
本实施例的步骤S23中，输出可看作DNA序列的高维语义编码，其中有一些冗余的转录调控语法规则。为生成非冗余的序列编码/>通过以卷积层构建的“序列高阶编码器”将/>映射到低维空间中，得到DNA一级序列的高阶特征/>为：
式中，ELU(·)为ELU激活函数，Wf和bf分别为序列高阶编码器的权重矩阵和截距向量，Conv1(·)为1维卷积操作。
本实施例的步骤S24中，基于DNA序列的高阶特征通过以全连接层构建的“染色质可达性预测器”预测当前DNA序列在各个细胞中的染色质可达性大小为：
式中，σ(·)为Sigmoid激活函数，Wp和bp分为染色质可达性预测器的权重矩阵和截距向量，其中，权重矩阵作为全部细胞的功能状态和高阶特征矩阵。
本发明实施例的步骤S3具体为：
S31、根据学习的功能状态和高阶特征构建k-近邻图；
其中，在构建的k-近邻图中节点为细胞，边为细胞之间的相关性(由欧几里德距离定义)；
S32、在构建的k-近邻图中，划分每个细胞的类型，并进一步在二维空间中可视化细胞类型注释结果。
具体地，在步骤S32中，在k-近邻图的基础上运用鲁汶算法划分每个细胞的类型，从而实现细胞类型的注释，最终，使用UMAP算法在二维空间中可视化细胞类型注释的结果。
在本发明实施例的步骤S4中，为了修正原始scATAC-seq数据由实验误差(特征峰检测丢失)导致的“假零”计数，基于构建的DNA语言模型，本发明实施例设计了一种染色质可达性图谱降噪算法。具体而言，染色质可达性图谱降噪算法首先采用统计模型搜索由“特征峰检测丢失”现象导致的“假零”计数，然后通过DNA语言模型预测的染色质可达性修正搜索到的“假零”计数。因此，本发明实施例的步骤S4具体为：
S41、获取原始ATAC-seq数据，即染色质可达性特征峰矩阵；
S42、通过混合概率模型搜索染色质可达性特征峰矩阵中由实验误差导致的零计数，进而计算特征峰i在细胞类型k的细胞m中的丢失率；
S43、对于给定的细胞m，根据其每个特征峰的丢失率，将全部特征峰划分为待修正的特征峰集合Am和无需修正的特征峰集合Bm；
S44、将通过DNA语言模型预测的特征峰i在每个细胞中的染色质可达性大小作为候选修正计数，对特征峰集合Am中的特征峰进行修正，进而获得降噪后的ATAC-seq特征峰的染色质可达性图谱。
本实施例的步骤S42中的混合概率模型包括伽马分布模型和正态分布模型，其中伽马分布模型表示实验误差，正态分布模型表示scATAC-seq特征峰的真实计数。由于不同细胞类型中每个特征峰伽马分布(实验误差)和正态分布(真实计数)的占比有所差异，对不同的细胞类型分别构建混合概率模型。对于细胞类型k中的特征峰i，它的计数可通过随机变量的密度函数表示：
式中，表示特征峰i在细胞类型k中的整体丢失率，/>和/>表示伽马分布的形状参数和尺度参数，/>和/>表示正态分布的均值和标准差；通过最大期望值算法估计混合概率模型中的参数。
由此得到本实施例步骤S42中，特征峰i在细胞类型k的细胞m中的丢失率dim为：
式中，表示参数的估计值，/>为特征峰i在细胞类型k中的整体丢失率，和/>为伽马分布的形状参数和尺度参数，/>和/>为正态分布的均值和标准差。
本实施例的步骤S43中，待修正的特征峰集合Am＝{i:dim≥T}和无需修正的特征峰集合Bm＝{i:dim<T}，其中T是人工设置的阈值。
本实施例的步骤S44中，假定特征峰i在细胞m中的候选修正计数为ri,m，只对集合Am中的特征峰进行修正，表达式为：
本发明实施例的步骤S5中基于构建的DNA语言模型，采用转录因子-DNA结合基元插入法，分析每个转录因子在单个细胞中的活跃性；具体过程为：
S51、随机选择s条scATAC-seq特征峰序列进行二核苷酸打乱生成s条的背景序列；
S52、设定一个转录因子，将其与DNA结合基元插入到s条背景序列中心位置，生成s条合成序列；
S53、通过DNA模型分别预测s条背景序列和s条合成序列的染色质可达性；
S54、将单个细胞的背景序列和合成序列的平均染色质可达性预测差值作为当前转录因子在该单个细胞中的活跃性分析结果。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (9)

1.一种单细胞ATAC-seq数据分析方法，其特征在于，包括以下步骤：
S1、采集ATAC-seq特征峰的DNA一级序列，作为ATAC-seq数据集；
S2、通过基于ProbDep Transformer的DNA语言模型对DNA一级序列进行分析，预测DNA一级序列中各细胞的染色质可达性，并学习每个细胞的功能状态和高阶特征；
S3、根据学习的功能状态和高阶特征，进行细胞类型注释；
S4、基于预测的染色质可达性，进行ATAC-seq特征峰的染色质可达性图谱降噪；
S5、通过DNA语言模型，分析ATAC-seq特征峰中每个转录因子在各细胞中的活跃性；
S6、将染色质可达性预测结果、细胞类型注释、降噪的染色质可达性图谱以及转录因子的活跃性分析结果作为单细胞ATAC-seq数据分析结果；
所述步骤S2具体为：
S21、将长度为L的DNA一级序列采用独热编码映射至维数为4×L的隐式特征空间中，并将其转换为基元编码矩阵；
对DNA一级序列采用绝对位置编码生成维数为pos×2i的位置编码矩阵，将基元编码矩阵和位置编码矩阵相加作为DNA语言模型的输入数据；其中，pos为当前转录因子-DNA结合基元在DNA一级序列中的位置下标，2i为当前转录因子-DNA结合基元的位置编码向量的长度；
S22、在DNA模型中，采用长距离依赖性测量评估查询与键之间的依赖性的方法对输入数据进行分析，获得每个键向量聚焦于排名最高的u个查询向量，进而得到DNA模型的输出数据；
S23、将DNA模型的输出数据作为DNA一级序列的高维语义编码，通过序列高阶编码器将其映射到低维空间中，获得DNA一级序列的高阶特征；
S24、根据获得的高阶特征，通过染色质可达性预测器预测当前DNA一级序列在各细胞的染色质可达性大小，并学习得到每个细胞的功能状态和高阶特征。
2.根据权利要求1所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S22中，长距离依赖性测量的表达式为：
式中，为长距离依赖性测量操作，/>表示qi向量与全部键向量之间经过Log-Sum-Exp操作之后的结果，/>表示Log-Sum-Exp结果的算术平均值，qi为查询矩阵Q的第i行，K为键矩阵，In为求对数操作，l为键矩阵K中当前行的下标，Lk为键矩阵K中行的个数，/>为键矩阵K中的第k行的转置，d为键矩阵K中列的个数；
DNA模型的输出数据表示为：
式中，为自注意力机制操作的输出，/>为自注意力机制操作，Q,K,V分别为查询矩阵，键矩阵，值矩阵，Softmax(·)为激活函数，/>为与矩阵Q大小相同的稀疏矩阵，且只包含了了长距离依赖性测量中排名最高的u个查询向量。
3.根据权利要求2所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S23中，DNA一级序列的高阶特征为：
式中，ELU(·)为ELU激活函数，Wf和bf分别为序列高阶编码器的权重矩阵和截距向量，Conv1(·)为1维卷积操作。
4.根据权利要求3所述的单细胞ATAC-seq数据分析方法，其特征在于，当前DNA一级序列在各细胞的染色质可达性大小y为：
式中，σ(·)为Sigmoid激活函数，Wp和bp分为染色质可达性预测器的权重矩阵和截距向量，其中，权重矩阵作为全部细胞的功能状态和高阶特征矩阵。
5.根据权利要求1所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S3具体为：
S31、根据学习的功能状态和高阶特征构建k-近邻图；
其中，在构建的k-近邻图中节点为细胞，边为细胞之间的相关性；
S32、在构建的k-近邻图中，划分每个细胞的类型，并进一步在二维空间中可视化细胞类型注释结果。
6.根据权利要求1所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S4具体为：
S41、获取原始ATAC-seq数据，即染色质可达性特征峰矩阵；
S42、通过混合概率模型搜索染色质可达性特征峰矩阵中由实验误差导致的零计数，进而计算特征峰i在细胞类型k的细胞m中的丢失率；
S43、对于给定的细胞m，根据其每个特征峰的丢失率，将全部特征峰划分为待修正的特征峰集合Am和无需修正的特征峰集合Bm；
S44、将通过DNA语言模型预测的特征峰i在每个细胞中的染色质可达性大小作为候选修正计数，对特征峰集合Am中的特征峰进行修正，进而获得降噪后的ATAC-seq特征峰的染色质可达性图谱。
7.根据权利要求6所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S42中的混合概率模型包括伽马分布模型和正态分布模型，其中伽马分布模型表示实验误差，正态分布模型表示scATAC-seq特征峰的真实计数。
8.根据权利要求6所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S42中，特征峰i在细胞类型k的细胞m中的丢失率dim为：
式中，表示参数/>的估计值，/>为特征峰i在细胞类型k中的整体丢失率，/>和/>为伽马分布的形状参数和尺度参数，/>和/>为正态分布的均值和标准差。
9.根据权利要求2所述的单细胞ATAC-seq数据分析方法，其特征在于，所述步骤S5具体为：
S51、随机选择s条scATAC-seq特征峰序列进行二核苷酸打乱生成s条的背景序列；
S52、设定一个转录因子，将其与DNA结合基元插入到s条背景序列中心位置，生成s条合成序列；
S53、通过DNA模型分别预测s条背景序列和s条合成序列的染色质可达性；
S54、将单个细胞的背景序列和合成序列的平均染色质可达性预测差值作为当前转录因子在该单个细胞中的活跃性分析结果。
CN202310182496.8A 2023-02-28 2023-02-28 一种单细胞ATAC-seq数据分析方法 Active CN116153404B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310182496.8A CN116153404B (zh) 2023-02-28 2023-02-28 一种单细胞ATAC-seq数据分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310182496.8A CN116153404B (zh) 2023-02-28 2023-02-28 一种单细胞ATAC-seq数据分析方法

Publications (2)

Publication Number Publication Date
CN116153404A CN116153404A (zh) 2023-05-23
CN116153404B true CN116153404B (zh) 2023-08-15

Family

ID=86356181

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310182496.8A Active CN116153404B (zh) 2023-02-28 2023-02-28 一种单细胞ATAC-seq数据分析方法

Country Status (1)

Country Link
CN (1) CN116153404B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116884495B (zh) * 2023-08-07 2024-03-08 成都信息工程大学 一种基于扩散模型的长尾染色质状态预测方法

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110910950A (zh) * 2019-11-18 2020-03-24 广州竞远生物科技有限公司 一种联合分析单细胞scRNA-seq和scATAC-seq的流程方法
CN111312329A (zh) * 2020-02-25 2020-06-19 成都信息工程大学 基于深度卷积自动编码器的转录因子结合位点预测的方法
CN111353995A (zh) * 2020-03-31 2020-06-30 成都信息工程大学 一种基于生成对抗网络的宫颈单细胞图像数据生成方法
WO2020198942A1 (zh) * 2019-03-29 2020-10-08 中国科学技术大学 基于峰聚类的单细胞染色质可及性测序数据分析方法和系统
CN112735514A (zh) * 2021-01-18 2021-04-30 清华大学 神经网络提取调控dna组合模式的训练和可视化方法及系统
CN112992267A (zh) * 2021-04-13 2021-06-18 中国人民解放军军事科学院军事医学研究院 一种单细胞的转录因子调控网络预测方法及装置
CN113490986A (zh) * 2018-12-31 2021-10-08 辉达公司 使用深度学习对atac-seq数据进行去噪
CN114836536A (zh) * 2022-07-04 2022-08-02 北京大学第三医院（北京大学第三临床医学院） 一种基于malbac的单细胞高扩增区域筛选方法及系统
WO2022221593A1 (en) * 2021-04-15 2022-10-20 Illumina, Inc. Efficient voxelization for deep learning
CN115273966A (zh) * 2022-08-29 2022-11-01 西安交通大学 谱系树中可变剪接模式和染色质状态动态变化的分析方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200370112A1 (en) * 2019-05-23 2020-11-26 The Board Of Trustees Of The Leland Stanford Junior University Methods utilizing single cell genetic data for cell population analysis and applications thereof
US20210118524A1 (en) * 2019-10-16 2021-04-22 Camp4 Therapeutics Corporation Identifying Therapeutic Targets by Discovering Interactions between Non-Coding Variants and Candidate Genes
US20210332354A1 (en) * 2020-04-15 2021-10-28 10X Genomics, Inc. Systems and methods for identifying differential accessibility of gene regulatory elements at single cell resolution
US20220068438A1 (en) * 2020-08-27 2022-03-03 The Broad Institute, Inc. Deep learning and alignment of spatially-resolved whole transcriptomes of single cells
US20220399129A1 (en) * 2021-06-15 2022-12-15 Flagship Pioneering Innovations Vi, Llc Systems and methods for terraforming

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113490986A (zh) * 2018-12-31 2021-10-08 辉达公司 使用深度学习对atac-seq数据进行去噪
WO2020198942A1 (zh) * 2019-03-29 2020-10-08 中国科学技术大学 基于峰聚类的单细胞染色质可及性测序数据分析方法和系统
CN110910950A (zh) * 2019-11-18 2020-03-24 广州竞远生物科技有限公司 一种联合分析单细胞scRNA-seq和scATAC-seq的流程方法
CN111312329A (zh) * 2020-02-25 2020-06-19 成都信息工程大学 基于深度卷积自动编码器的转录因子结合位点预测的方法
CN111353995A (zh) * 2020-03-31 2020-06-30 成都信息工程大学 一种基于生成对抗网络的宫颈单细胞图像数据生成方法
CN112735514A (zh) * 2021-01-18 2021-04-30 清华大学 神经网络提取调控dna组合模式的训练和可视化方法及系统
CN112992267A (zh) * 2021-04-13 2021-06-18 中国人民解放军军事科学院军事医学研究院 一种单细胞的转录因子调控网络预测方法及装置
WO2022221593A1 (en) * 2021-04-15 2022-10-20 Illumina, Inc. Efficient voxelization for deep learning
CN114836536A (zh) * 2022-07-04 2022-08-02 北京大学第三医院（北京大学第三临床医学院） 一种基于malbac的单细胞高扩增区域筛选方法及系统
CN115273966A (zh) * 2022-08-29 2022-11-01 西安交通大学 谱系树中可变剪接模式和染色质状态动态变化的分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
不平衡分类的数据采样方法综述;刘定祥1a，乔少杰，张永清，韩楠等;重庆理工大学学报( 自然科学);第33卷(第7期);全文 *

Also Published As

Publication number Publication date
CN116153404A (zh) 2023-05-23

Similar Documents

Publication Publication Date Title
Chen et al. Efficient ant colony optimization for image feature selection
Guo et al. A feature fusion based forecasting model for financial time series
Zou et al. Regularized simultaneous model selection in multiple quantiles regression
CN108446741B (zh) 机器学习超参数重要性评估方法、系统及存储介质
CN116153404B (zh) 一种单细胞ATAC-seq数据分析方法
Ding et al. Intelligent optimization methods for high-dimensional data classification for support vector machines
US20230207064A1 (en) Inter-model prediction score recalibration during training
Pulikottil et al. Onet–a temporal meta embedding network for mooc dropout prediction
Tadayon et al. A clustering approach to time series forecasting using neural networks: A comparative study on distance-based vs. feature-based clustering methods
US20220336056A1 (en) Multi-channel protein voxelization to predict variant pathogenicity using deep convolutional neural networks
US11515010B2 (en) Deep convolutional neural networks to predict variant pathogenicity using three-dimensional (3D) protein structures
Lu et al. Entssr: a weighted ensemble learning method to impute single-cell RNA sequencing data
KR20230170680A (ko) 심층 콘볼루션 신경망들을 사용하여 변이체 병원성을 예측하기 위한 다중 채널 단백질 복셀화
KR20230171930A (ko) 3차원(3d) 단백질 구조들을 사용하여 변이체 병원성을 예측하기 위한 심층 콘볼루션 신경망들
Ye et al. Surrogate-based Global Optimization Methods for Expensive Black-Box Problems: Recent Advances and Future Challenges
Sulieman et al. A supervised feature selection approach based on global sensitivity
CN112735596A (zh) 一种相似患者的确定方法、装置、电子设备及存储介质
Jo et al. How Much a Model be Trained by Passive Learning Before Active Learning?
Wang et al. AutoST: Training-free Neural Architecture Search for Spiking Transformers
Islam et al. Leveraging cell-cell similarity for high-performance spatial and temporal cellular mappings from gene expression data
Ao et al. Gene expression time series modeling with principal component and neural network
US20230223100A1 (en) Inter-model prediction score recalibration
US20230047347A1 (en) Deep neural network-based variant pathogenicity prediction
Argiento et al. A Bayesian nonparametric mixture model for cluster analysis
Marion et al. VC-PCR: A Prediction Method based on Supervised Variable Selection and Clustering

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