CN112509639B - 一种单细胞基因融合检测方法 - Google Patents

一种单细胞基因融合检测方法 Download PDF

Info

Publication number
CN112509639B
CN112509639B CN202011451710.8A CN202011451710A CN112509639B CN 112509639 B CN112509639 B CN 112509639B CN 202011451710 A CN202011451710 A CN 202011451710A CN 112509639 B CN112509639 B CN 112509639B
Authority
CN
China
Prior art keywords
fusion
gene
potential gene
potential
sequences
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
CN202011451710.8A
Other languages
English (en)
Other versions
CN112509639A (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.)
Peking University
Original Assignee
Peking University
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 Peking University filed Critical Peking University
Priority to CN202011451710.8A priority Critical patent/CN112509639B/zh
Publication of CN112509639A publication Critical patent/CN112509639A/zh
Application granted granted Critical
Publication of CN112509639B publication Critical patent/CN112509639B/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
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search
    • 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/30—Detection of binding sites or motifs
    • 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

Landscapes

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

Abstract

本发明公开了一种单细胞基因融合检测方法。该方法包括如下步骤：将单细胞RNA测序数据进行序列比对和映射，分析映射结果，获得潜在基因融合；建立零点膨胀负二项分布统计模型和训练神经网络，获得每个潜在基因融合显著性和可靠性的评价；根据可靠性以及融合的特征，报告最终的融合基因。本发明可以敏感且精确地检测单细胞基因融合，降低假阳性或漏报率。本发明具有重大的应用价值。

Description

一种单细胞基因融合检测方法
技术领域
本发明属于生物信息学领域，具体涉及一种单细胞基因融合检测方法。
背景技术
融合基因是指由DNA序列重排产生的，包含2个或2个以上基因序列的嵌合体。融合基因是肿瘤基因组中的一种常见变异，很多融合基因是癌症发展的驱动变异，对肿瘤的发生、发展有重要作用。融合基因是一些肿瘤诊断和分类的主要依据，也是癌症药物的重要靶点。由于肿瘤细胞有很强的异质性，病人的不同癌细胞可能会有不同的变异，准确检测肿瘤单细胞的基因变异对研究肿瘤细胞的起源、进化、抗药性乃至肿瘤的治疗都有重要的意义。
单细胞RNA测序是在单个细胞的分辨率上进行高通量测序分析的一项技术，该技术的出现使得高通量大规模单细胞基因融合的检测成为可能。然而，利用已有分析普通RNA测序数据的基因融合检测方法直接检测单细胞基因融合仍有很多问题，主要包括：一、单细胞测序数据有很高的噪声，导致给出的结果的精度很低，得到很多假阳性；二、单独对每个单细胞测序数据独立检测会有敏感度低，漏检率高的问题；三、花费的计算时间和计算资源巨大。由此可见，现有的检测工具不能满足单细胞RNA测序数据的基因融合检测需求，亟需发展利用单细胞RNA测序检测单细胞基因融合的方法。
发明内容
本发明的目的是检测单细胞的基因融合。
本发明首先保护一种单细胞基因融合检测方法，可包括如下步骤：
(1)将单细胞RNA测序数据比对到参考基因组，得到的支持融合的分离序列和不统一序列；
(2)完成步骤(1)后，从分离序列和不统一序列中整合所有的潜在基因融合；
(3)完成步骤(2)后，采用统计模型刻画背景噪声的支持分离序列数的分布，获得每个潜在基因融合的显著性指标；
(4)完成步骤(2)后，训练神经网络学习背景噪声的序列特征，获得每个潜在基因融合为背景噪声的概率；
(5)根据步骤(3)获得的显著性指标和步骤(4)获得的概率，筛选获得显著可靠的潜在基因融合；
(6)从步骤(5)获得的显著可靠的潜在基因融合中筛选符合标准的基因融合，即单细胞的基因融合。
上述方法中，所述单细胞RNA测序数据可为2个以上单细胞的RNA测序数据。
上述任一所述单细胞可为T细胞。
所述步骤(1)中，所述参考基因组可为参考基因组hg19。
所述步骤(1)中，采用STAR软件(版本2.7.4a)进行比对。STAR软件(版本2.7.4a)比对后自动获得分离序列和不统一序列。
所述步骤(1)中，分离序列和不统一序列比对的位置均为基因组上连续75bp以上碱基序列唯一的区域。这样可以防止因不同区域的序列相同导致比对位置错误。
所述步骤(2)中，从分离序列和不统一序列中整合所有的潜在基因融合的步骤可如下：
(2-1)获得分离序列具体的比对位置(根据步骤(1)的结果获得)，即对应的潜在基因融合位点；
(2-2)将距离不超过3bp的融合位点看作是同一个融合位点，归并得到一个潜在基因融合的列表，并且记录支持每一个基因融合的全部分离序列与不统一序列的数量；
(2-3)去掉只有一个分离序列支持的潜在基因融合；
(2-4)从测序结果中统计每个潜在基因融合相关的两个基因在每个细胞中的表达量，以及融合位点附近200bp的GC含量。
所述步骤(2-4)中，表达量可定义为这个基因的覆盖序列数除以100再加1，然后取自然对数的值。
所述步骤(3)中，所述统计模型可为零点膨胀的负二项分布模型。所述零点膨胀的负二项分布模型中重参数化为取零值的概率与零点截断的负二项分布。
所述步骤(3)中，把潜在基因融合均看作背景噪声来训练出背景噪声的分布。这是由于潜在基因融合的数量远大于真实的基因融合，因此可以把潜在基因融合均看作背景噪声来训练出背景噪声的分布。
所述步骤(3)中，利用分离序列的支持总数来评价潜在基因融合的显著性；利用重抽样从零点膨胀的负二项分布中抽取多个背景噪声的分离序列的支持总数，建立正态分布，得到p值；将潜在基因融合划分为两个集合，通过较差集合中的数据来估计出总体错误发现率，然后选择合适的p值的阈值来控制错误发现率。
所述步骤(3)中，零点膨胀负二项分布模型对于每一个潜在基因融合包含三种参数pi，μij，λ，分别为取零值的概率pi、不取零值时分离序列支持数的期望μij以及分散系数(overdispersion)λ。对于同一个潜在基因融合，每个细胞取零值的概率是相同的，可以用关于GC含量的样条函数f(GCi)和两个相关基因的平均表达量
Figure BDA0002827266100000021
和
Figure BDA0002827266100000022
的可加模型来估计。对于同一个潜在基因融合，每个细胞不取零值分离序列支持数的期望是不同的，可以用关于GC含量的自由度为5的三次样条函数g(GCi)和两个相关基因在该细胞中的表达量ELij和EHij的可加模型来估计。分散系数设为一个未知常数。该模型可以写为如下表达：
Yij～ZINB(pi，μij，λ)；
Figure BDA0002827266100000031
log(μij)＝β20+g(Gci)+β21ELij+β22EHij。
所述步骤(3)中，采用零点膨胀的负二项分布模型刻画背景噪声的支持分离序列数的分布，获得每个潜在基因融合的显著性指标的步骤可如下：
(3-1)对所有潜在基因融合的每个细胞的分离序列支持数套用零点膨胀负二项分布模型，计算出关于未知参数的似然函数l：
Figure BDA0002827266100000032
Figure BDA0002827266100000033
(3-2)利用R中的optim函数以及默认参数求解出最优的参数；
(3-3)计算每个潜在基因融合在不同细胞中的分离序列支持总数；
(3-4)将步骤(3-2)中估计出的参数代入零点膨胀负二项分布中，然后从模型中重抽样1000个以上，计算这1000个以上分离序列支持总数的均值
Figure BDA0002827266100000034
和标准差
Figure BDA0002827266100000035
建立分离序列支持总数的正态分布，获得p值：
Figure BDA0002827266100000036
Φ为标准正态分布的分布函数；
(3-5)将有分离序列支持细胞数不小于1％的总细胞数并且有支持的细胞中平均支持数超过1.25的潜在基因融合划分到集合1中，其余的划分到集合2中；选取合适p值的阈值，使得集合2中的p值小于阈值的潜在基因融合数量乘以总潜在基因融合数除以集合2中的潜在基因融合数再除以两个集合中p值小于阈值的潜在基因融合数不超过给定的界；p值小于该阈值的所有潜在基因融合，被认为是显著的潜在基因融合。
所述步骤(4)中，神经网络为双向长短时记忆网络(bi-LSTM)。所述双向长短时记忆网络包括编码、四个长短时记忆层和两个全连接层。
所述步骤(4)中，输入是60碱基长的序列及其断点，输出是该潜在基因融合为假融合的概率。网络的建立是将序列编码(A编码为向量(1，0，0，0，0)，C编码为向量(0，1，0，0，0)，T编码为向量(0，0，1，0，0)，G编码为向量(0，0，0，1，0)，断点编码为向量(0，0，0，0，1))，每个序列编码为一个61*5维的矩阵。然后将这个矩阵依次通过三个bi-LSTM层(每个层分别有32、64、128个bi-LSTM单元)，然后通过一个输入为序列，输出为实数的bi-LSTM层，最后通过两个全连接层(分别是256维和2维)得到估计的概率值。
所述步骤(4)中，训练神经网络学习背景噪声的序列特征，获得每个潜在基因融合为背景噪声的概率的步骤可如下：
(4-1)根据输入的数据集建立正训练样本和负训练样本；正训练样本从潜在基因融合中抽取，负训练样本为没有噪声序列特征的随机拼接序列；
(4-2)在预训练模型的基础上进行重训练直至收敛；
(4-3)对每个潜在基因融合进行预测，得到概率。
所述步骤(4-1)中，正训练样本可从潜在基因融合中抽取10％的样本，得到断点附近的60bp长的序列，并且在断点处加上标记。负训练样本可从转录本中随机抽取的序列拼接而成，得到60bp长的相同数量的序列，也在断点处加了标记。
所述步骤(4-3)中，概率为该融合为背景噪声的概率。
所述步骤(6)中，所述标准同时满足①-⑥：①构成融合基因的基因不能是伪基因；②构成融合基因的基因不能是非编码RNA；③构成融合基因的任一基因没有被接受的代号；④融合基因涉及到的断点需要在外显子区域；⑤参与基因融合的基因不能出现在超过5个不同的基因融合中；⑥支持基因融合的不一致序列的数量不能大于分离序列的数量的10倍。
在本发明的实施例2中，所述单细胞RNA测序数据具体可为2355个T细胞的单细胞RNA测序数据。取步骤(3-5)中的界为0.05，得到显著的潜在基因融合；取步骤(4-3)中的界为0.75、保留噪声概率低于75％的潜在基因融合；此时筛选获得的即为显著可靠的潜在基因融合。最终获得单细胞的基因融合见表1。
本发明建立的检测单细胞的基因融合的方法，率先考虑了单细胞测序数据的噪声特点，利用零点膨胀的负二项分布模型结合多个细胞的数据来检测出显著的基因融合；与bulk RNA测序数据相比，极大地降低了假阳性，提高了精度。同时，利用神经网络学习噪声的序列特征，通过检测断点附近的序列是否有噪声的特征去掉假阳性，提高精度。本发明建立的方法花费的计算资源相比传统的bulk方法更加经济、高效。本发明在检测单细胞的基因融合上具有重要的价值，尤其是在数据集中存在多个细胞类型或无法获取大量细胞的情况下具有较高的应用前景。
具体实施方式
下面结合具体实施方式对本发明进行进一步的详细描述，给出的实施例仅为了阐明本发明，而不是为了限制本发明的范围。以下提供的实施例可作为本技术领域普通技术人员进行进一步改进的指南，并不以任何方式构成对本发明的限制。
下述实施例中的实验方法，如无特殊说明，均为常规方法，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、单细胞基因融合检测方法的建立
本发明的发明人经过大量实验，建立了一种单细胞基因融合检测方法。具体步骤如下：
1、采用STAR软件(版本2.7.4a)将单细胞RNA测序数据(双端)比对到参考基因组(如参考基因组hg19)上，得到的支持融合的分离序列(split-mappedread)和不统一序列(discordantread)。
STAR软件(版本2.7.4a)比对后自动获得分离序列和不统一序列。
进行步骤1时，为了防止因不同区域的序列相同导致比对位置错误，仅考虑连续75bp的碱基序列在整个基因组上唯一存在的区域。
2、完成步骤1后，从分离序列和不统一序列中整合出所有的潜在基因融合。
具体包括如下步骤：
(2-1)根据步骤1的结果，获得分离序列具体的比对位置，即对应的潜在基因融合位点；
(2-2)将距离不超过3bp的融合位点看作是同一个融合位点，归并得到一个潜在基因融合的列表，并且记录支持每一个基因融合的全部分离序列与不统一序列的数量。
(2-3)去掉只有一个分离序列支持的潜在基因融合。
(2-4)从测序结果中统计每个潜在基因融合相关的两个基因在每个细胞中的表达量，以及融合位点附近200bp的GC含量。
表达量定义为这个基因的覆盖序列数除以100再加1，然后取自然对数的值。
由于潜在基因融合的数量远大于真实的基因融合，因此可以近似地把全部的潜在基因融合都认为是背景噪声。把潜在基因融合均看作背景噪声来训练出背景噪声的分布。
3、完成步骤2后，采用零点膨胀的负二项分布模型(ZINB)刻画背景噪声的支持分离序列数的分布，获得每个潜在基因融合的显著性指标。
零点膨胀负二项分布模型对于每一个潜在基因融合包含三种参数pi，μij，λ，分别为取零值的概率pi、不取零值时分离序列支持数的期望μij以及分散系数(overdispersion)λ。对于同一个潜在基因融合，每个细胞取零值的概率是相同的，可以用关于GC含量的样条函数f(GCi)和两个相关基因的平均表达量
Figure BDA0002827266100000051
和
Figure BDA0002827266100000052
的可加模型来估计。对于同一个潜在基因融合，每个细胞不取零值分离序列支持数的期望是不同的，可以用关于GC含量的自由度为5的三次样条函数g(Gci)和两个相关基因在该细胞中的表达量ELij和EHij的可加模型来估计。分散系数设为一个未知常数。该模型可以写为如下表达：
Yij～ZINB(pi，μij，λ)；
Figure BDA0002827266100000061
log(μij)＝β20+g(GCi)+β21ELij+β22EHij。
具体包括如下步骤：
(3-1)对所有潜在基因融合的每个细胞的分离序列支持数套用零点膨胀负二项分布模型，计算出关于未知参数的似然函数l：
Figure BDA0002827266100000062
Figure BDA0002827266100000063
(3-2)利用R中的optim函数以及默认参数求解出最优的参数。
(3-3)计算每个潜在基因融合在不同细胞中的分离序列支持总数。
(3-4)将步骤(3-2)中估计出的参数代入零点膨胀负二项分布中，然后从此模型中重抽样1000次，计算这1000个分离序列支持总数的均值
Figure BDA0002827266100000064
和标准差
Figure BDA0002827266100000065
建立分离序列支持总数的正态分布，给出p值：
Figure BDA0002827266100000066
Φ为标准正态分布的分布函数。
(3-5)为了控制错误发现率(FDR)，将有分离序列支持细胞数不小于1％的总细胞数并且有支持的细胞中平均支持数超过1.25的潜在基因融合划分到集合1中，其余的划分到集合2中；选取合适p值的阈值，使得集合2中的p值小于阈值的潜在基因融合数量乘以总潜在基因融合数除以集合2中的潜在基因融合数再除以两个集合中p值小于阈值的潜在基因融合数不超过给定的界；p值小于该阈值的所有潜在基因融合，被认为是显著的潜在基因融合。
4、完成步骤2后，训练神经网络学习背景噪声的序列特征，获得每个潜在基因融合为背景噪声的概率，即可靠性指标。
神经网络是双向长短时记忆网络(bi-LSTM)。输入是60碱基长的序列及其断点，输出是该潜在基因融合为假融合的概率。网络的建立是将序列编码(A编码为向量(1，0，0，0，0)，C编码为向量(0，1，0，0，0)，T编码为向量(0，0，1，0，0)，G编码为向量(0，0，0，1，0)，断点编码为向量(0，0，0，0，1))，每个序列编码为一个61*5维的矩阵。然后将这个矩阵依次通过三个bi-LSTM层(每个层分别有32、64、128个bi-LSTM单元)，然后通过一个输入为序列，输出为实数的bi-LSTM层，最后通过两个全连接层(分别是256维和2维)得到估计的概率值。
具体包括如下步骤：
(4-1)根据输入的数据集建立训练样本。正训练样本从潜在基因融合中抽取10％的样本，得到断点附近的60bp长的序列，并且在断点处加上标记；负训练样本从转录本中随机抽取的序列拼接而成，得到60bp长的相同数量的序列，也在断点处加了标记。
(4-2)在预训练模型的基础上进行重训练直至收敛。
(4-3)对每个潜在基因融合进行预测，得到该融合为背景噪声的概率。
5、根据步骤3获得的显著的潜在基因融合和步骤4获得的可靠性指标，筛选获得显著可靠的潜在基因融合。
取步骤(3-5)中的界为0.05，得到显著的潜在基因融合；取步骤(4-3)中的界为0.75、保留噪声概率低于75％的潜在基因融合；此时筛选获得的即为显著可靠的潜在基因融合。
6、从步骤5获得的显著可靠的潜在基因融合中筛选符合标准的基因融合。
所述标准同时满足①-⑥：①构成融合基因的基因不能是伪基因；②构成融合基因的基因不能是非编码RNA；③构成融合基因的任一基因没有被接受的代号；④融合基因涉及到的断点需要在外显子区域；⑤参与基因融合的基因不能出现在超过5个不同的基因融合中；⑥支持基因融合的不一致序列的数量不能大于分离序列的数量的10倍。
实施例2、采用实施例1建立的方法检测2355个T细胞的基因融合
按照实施例1建立的方法检测2355个T细胞的基因融合。具体步骤如下：
1、获得2355个T细胞的单细胞RNA测序数据(双端)。
2355个T细胞的单细胞RNA测序数据(记载于如下文献中：QimingZhang,etal.Landscape and Dynamics of Single Immune Cells in HepatocellularCarcinoma.Cell.Volume 179,Issue 4,31 October 2019,Pages 829-845.)存放于中国科学院北京基因研究所的组学原始数据归档库(GSA)，编号为HRA000069。
2、完成步骤1后，将2355个T细胞双端测序数据文件重命名，命名形式为1_1.fastq、1_2.fastq、2_1.fastq、2_2.fastq、……2355_1.fastq和2355_2.fastq；之后采用STAR软件(版本2.7.4a)将测序序列映射到参考基因组hg19上，得到的支持融合的分离序列和不统一序列。
进行该步骤时，删除比对到基因组上的有重复序列超过75bp的区域的记录。
3、完成步骤2后，根据分离序列和不统一序列的列表，得到分离序列的断点位置以及不统一序列所支持的潜在基因融合。
将断点位置相同或距离不超过3bp的分离序列分为同一组，认为它们都是支持同一处的基因融合。将支持基因相同的不统一序列分为同一组，认为它们都是支持所有涉及这两个基因的基因融合的序列。将这些支持的信息统计为一个列表。列表中的每一行指代了一个潜在的基因融合的信息，包含基因名、断点在染色体上的坐标、每个细胞上支持的分离序列数和不统一序列数、每个细胞的这两个基因的表达、断点附近200bp的GC含量。
4、完成步骤3后，根据潜在基因融合的支持情况以及基因的表达与断点附近的GC含量，建立零点膨胀的负二项分布(ZINB)模型。计算整体的似然函数，利用R中的optim函数求解出最优参数，得到该分布的参数估计。
5、根据步骤4得到的分布，计算每个潜在基因融合的分离序列支持总数，作为描述这个基因融合的显著性指标。从估计好的ZINB模型中重抽样1000次，计算这1000个分离序列支持总数的均值和标准差。分离序列支持总数的分布可以认为是一个正态分布，这个正态分布的均值和标准差就用这1000个样本来估计。获得p值。
6、完成步骤5后，将有分离序列支持细胞数不小于1％的总细胞数并且有支持的细胞中平均支持数超过1.25的潜在基因融合划分到集合1中，其余的划分到集合2中。选取合适的阈值c，使得集合2中的p值小于c的潜在基因融合数量乘以总潜在基因融合数除以集合2中的潜在基因融合数再除以两个集合中p值小于c的潜在基因融合数不超过给定的0.05。之后报告出p值小于c的所有潜在基因融合，被认为是显著的潜在基因融合。
7、完成步骤6后，根据输入的数据集来建立训练样本。正训练样本从潜在基因融合中随机选取10％，抽取出断点附近的60bp长的序列，并且在断点出***标记；负训练样本从数据集的随机细胞中随机抽取的正常比对的部分序列，然后拼接成长度为60bp的序列，并且在断点处***标记。正负样本的数量相同。
8、完成步骤7后，加载预训练的网络，将正负样本传递给神经网络进行重训练。利用重训练之后的网络估计融合为噪声的概率。去除步骤6得到的显著的潜在基因融合中概率值为0.75以上的部分，此时获得的即为显著可靠的潜在基因融合。
9、完成步骤8后，根据对应的基因注释文件，可以知道每个基因的外显子区域、以及基因是否为非编码RNA、是否为伪基因。从步骤8获得的显著可靠的潜在基因融合中筛选符合标准的基因融合；
所述标准同时满足①-⑥：①构成融合基因的基因不能是伪基因；②构成融合基因的基因不能是非编码RNA；③构成融合基因的任一基因没有被接受的代号；④融合基因涉及到的断点需要在外显子区域；⑤参与基因融合的基因不能出现在超过5个不同的基因融合中；⑥支持基因融合的不一致序列的数量不能大于分离序列的数量的10倍。
最终输出结果见表1。
表1
序号 基因融合 位置1 位置2 噪声概率 FDR
1 BAG3--TRBC2 chr10:121429414 chr7:142499834 0.199 0
2 C15orf57--CBX3 chr15:40854180 chr7:26241365 0.484 0
3 TRAJ33--TRAV1-2 chr14:22977587 chr14:22111802 0.057 0
4 IGHV7-81--TRAC chr14:107282989 chr14:23016447 0.162 0
5 TRAJ12--TRAV1-2 chr14:23000892 chr14:22111802 0.068 0
6 TRAJ39--TRAV21 chr14:22970583 chr14:22521314 0.230 0
7 IGHV7-81--TRAJ41 chr14:107282788 chr14:22966644 0.148 0
8 TRAC--TRAV8-1 chr14:23016447 chr14:22265878 0.154 0
上述结果表明，采用实施例1建立的方法可以检测单细胞的基因融合。
实施例1建立的方法具有如下有益效果：
(1)精度高，假阳性低
由于实施例2使用的样本为免疫T细胞，因此基因融合的数量非常少，并且大部分都应该是和V(D)J重组相关的基因融合。一方面，本发明提供的方法只获得了8个基因融合，数量很少，并且其中6个都涉及到了V(D)J重组。而本发明的发明人采用四种常规方法(STAR-Fusion(v1.8.1)、Arriba(v1.0.1)、EricScript(v0.5.5b)和FusionCatcher(v1.10))分别报告了87、574、716和4306个基因融合，总数量远大于本发明的方法。另一方面，着重于TRAJ33--TRAV1-2融合，本发明提供的方法找到了128个细胞有该融合，STAR-Fusion(v1.8.1)找到了62个细胞，另外三个方法都没有报告这个基因融合。因此本发明提供的方法的精度高，假阳性低。
(2)运行时间短，更加经济、高效
本发明提供的方法使用了847个小时核心，STAR-Fusion(v1.8.1)、Arriba(v1.0.1)、EricScript(v0.5.5b)和FusionCatcher(v1.10)在同一个服务器上分别使用了3534个、1830个、8634个和22764个小时核心。由此可见，本发明提供的方法运行时间短，更加经济、高效。
以上对本发明进行了详述。对于本领域技术人员来说，在不脱离本发明的宗旨和范围，以及无需进行不必要的实验情况下，可在等同参数、浓度和条件下，在较宽范围内实施本发明。虽然本发明给出了特殊的实施例，应该理解为，可以对本发明作进一步的改进。总之，按本发明的原理，本申请欲包括任何变更、用途或对本发明的改进，包括脱离了本申请中已公开范围，而用本领域已知的常规技术进行的改变。按以下附带的权利要求的范围，可以进行一些基本特征的应用。

Claims (5)

1.一种单细胞基因融合检测方法，包括如下步骤：
(1)将单细胞RNA测序数据比对到参考基因组，得到的支持融合的分离序列和不统一序列；
(2)完成步骤(1)后，从分离序列和不统一序列中整合所有的潜在基因融合；
(3)完成步骤(2)后，采用统计模型刻画背景噪声的支持分离序列数的分布，获得每个潜在基因融合的显著性指标；
(4)完成步骤(2)后，训练神经网络学习背景噪声的序列特征，获得每个潜在基因融合为背景噪声的概率；
(5)根据步骤(3)获得的显著性指标和步骤(4)获得的概率，筛选获得显著可靠的潜在基因融合；
(6)从步骤(5)获得的显著可靠的潜在基因融合中筛选符合标准的基因融合；
所述步骤(3)中，所述统计模型为零点膨胀的负二项分布模型；所述零点膨胀的负二项分布模型中重参数化为取零值的概率与零点截断的负二项分布；把潜在基因融合均看作背景噪声来训练出背景噪声的分布；
所述步骤(3)中，利用分离序列的支持总数来评价潜在基因融合的显著性；利用重抽样从零点膨胀的负二项分布中抽取多个背景噪声的分离序列的支持总数，建立正态分布，得到p值；将潜在基因融合划分为两个集合，通过较差集合中的数据来估计出总体错误发现率，然后选择合适的p值的阈值来控制错误发现率；
所述步骤(3)中，采用零点膨胀的负二项分布模型刻画背景噪声的支持分离序列数的分布，获得每个潜在基因融合的显著性指标的步骤如下：
(3-1)对所有潜在基因融合的每个细胞的分离序列支持数套用零点膨胀负二项分布模型，计算出关于未知参数的似然函数l:
Figure FDA0003598567050000011
Figure FDA0003598567050000012
(3-2)利用R中的optim函数以及默认参数求解出最优的参数；
(3-3)计算每个潜在基因融合在不同细胞中的分离序列支持总数；
(3-4)将步骤(3-2)中估计出的参数代入零点膨胀负二项分布中，然后从模型中重抽样1000个以上，计算这1000个以上分离序列支持总数的均值
Figure FDA0003598567050000013
和标准差
Figure FDA0003598567050000014
建立分离序列支持总数的正态分布，获得p值：
Figure FDA0003598567050000021
Φ为标准正态分布的分布函数；
(3-5)将有分离序列支持细胞数不小于1％的总细胞数并且有支持的细胞中平均支持数超过1.25的潜在基因融合划分到集合1中，其余的划分到集合2中；选取合适p值的阈值，使得集合2中的p值小于阈值的潜在基因融合数量乘以总潜在基因融合数除以集合2中的潜在基因融合数再除以两个集合中p值小于阈值的潜在基因融合数不超过给定的界；p值小于该阈值的所有潜在基因融合，被认为是显著的潜在基因融合；
所述步骤(4)中，神经网络为双向长短时记忆网络；所述双向长短时记忆网络包括编码、四个长短时记忆层和两个全连接层；
所述步骤(4)中，训练神经网络学习背景噪声的序列特征，获得每个潜在基因融合为背景噪声的概率的步骤如下：
(4-1)根据输入的数据集建立正训练样本和负训练样本；正训练样本从潜在基因融合中抽取，负训练样本为没有噪声序列特征的随机拼接序列；
(4-2)在预训练模型的基础上进行重训练直至收敛；
(4-3)对每个潜在基因融合进行预测，得到概率。
2.如权利要求1所述的方法，其特征在于：所述单细胞RNA测序数据为2个以上单细胞的RNA测序数据。
3.如权利要求1所述的方法，其特征在于：所述步骤(1)中，分离序列和不统一序列比对的位置均为基因组上连续75bp以上碱基序列唯一的区域。
4.如权利要求1所述的方法，其特征在于：所述步骤(2)中，从分离序列和不统一序列中整合所有的潜在基因融合的步骤如下：
(2-1)获得分离序列具体的比对位置，即对应的潜在基因融合位点；
(2-2)将距离不超过3bp的融合位点看作是同一个融合位点，归并得到一个潜在基因融合的列表，并且记录支持每一个基因融合的全部分离序列与不统一序列的数量；
(2-3)去掉只有一个分离序列支持的潜在基因融合；
(2-4)从测序结果中统计每个潜在基因融合相关的两个基因在每个细胞中的表达量，以及融合位点附近200bp的GC含量。
5.如权利要求1所述的方法，其特征在于：所述步骤(6)中，所述标准同时满足①-⑥：①构成融合基因的基因不能是伪基因；②构成融合基因的基因不能是非编码RNA；③构成融合基因的任一基因没有被接受的代号；④融合基因涉及到的断点需要在外显子区域；⑤参与基因融合的基因不能出现在超过5个不同的基因融合中；⑥支持基因融合的不一致序列的数量不能大于分离序列的数量的10倍。
CN202011451710.8A 2020-12-10 2020-12-10 一种单细胞基因融合检测方法 Active CN112509639B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011451710.8A CN112509639B (zh) 2020-12-10 2020-12-10 一种单细胞基因融合检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011451710.8A CN112509639B (zh) 2020-12-10 2020-12-10 一种单细胞基因融合检测方法

Publications (2)

Publication Number Publication Date
CN112509639A CN112509639A (zh) 2021-03-16
CN112509639B true CN112509639B (zh) 2022-05-31

Family

ID=74971957

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011451710.8A Active CN112509639B (zh) 2020-12-10 2020-12-10 一种单细胞基因融合检测方法

Country Status (1)

Country Link
CN (1) CN112509639B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111564183A (zh) * 2020-04-24 2020-08-21 西北工业大学 融合基因本体和神经网络的单细胞测序数据降维方法
WO2020208181A1 (en) * 2019-04-12 2020-10-15 European Molecular Biology Laboratory Comprehensive detection of single cell genetic structural variations

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11127484B2 (en) * 2016-03-17 2021-09-21 Regents Of The University Of Minnesota Analysis of single cell transcriptomics
WO2017164936A1 (en) * 2016-03-21 2017-09-28 The Broad Institute, Inc. Methods for determining spatial and temporal gene expression dynamics in single cells

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020208181A1 (en) * 2019-04-12 2020-10-15 European Molecular Biology Laboratory Comprehensive detection of single cell genetic structural variations
CN111564183A (zh) * 2020-04-24 2020-08-21 西北工业大学 融合基因本体和神经网络的单细胞测序数据降维方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Meng Duan et al..Diverse modes of clonal evolution in HBV-related hepatocellular carcinoma revealed by single-cell genome sequencing.《Cell Research》.2018, *
李羽翡 等.单细胞分析技术的发展及应用.《生命的化学》.2016,(第05期), *
李贱成 等.单细胞转录组测序技术及其应用.《生命的化学》.2020,(第08期), *

Also Published As

Publication number Publication date
CN112509639A (zh) 2021-03-16

Similar Documents

Publication Publication Date Title
Sun et al. SLAF-seq: an efficient method of large-scale de novo SNP discovery and genotyping using high-throughput sequencing
NZ759804A (en) Deep learning-based techniques for training deep convolutional neural networks
CN110400602B (zh) 一种基于测序数据的abo血型***分型方法及其应用
CN110997936B (zh) 基于低深度基因组测序进行基因分型的方法、装置及其用途
Keel et al. Genome‐wide copy number variation in the bovine genome detected using low coverage sequence of popular beef breeds
CN115428088A (zh) 用于基因表达和dna染色质可及性的联合交互式可视化的***和方法
CN108256289A (zh) 一种基于目标区域捕获测序基因组拷贝数变异的方法
CN115083521B (zh) 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及***
CN108268752B (zh) 一种染色体异常检测装置
CN108710784A (zh) 一种基因转录变异几率及变异方向的算法
CN109929934B (zh) 免疫相关基因在结直肠癌预后的试剂盒和***中的应用
CN116230081A (zh) 一种用于肺腺癌预后预测的生物标志物、应用及模型构建方法
CN114913919A (zh) 一种单基因病遗传变异智能解读及报告的方法、***及服务器
Luebker et al. Comparing the genomes of cutaneous melanoma tumors to commercially available cell lines
CN112509639B (zh) 一种单细胞基因融合检测方法
CN108197428A (zh) 一种并行动态规划的下一代测序技术拷贝数变异检测方法
CN117037905A (zh) 基于祖先信息标记的鸡品种鉴定方法、***、设备及介质
CN116486913A (zh) 基于单细胞测序从头预测调控突变的***、设备和介质
CN116564410A (zh) 一种预测突变位点顺式调控基因的方法、设备和介质
CN116200490A (zh) 一种检测实体瘤微小残留病灶的方法
He et al. Effects of Mycobacterium tuberculosis lineages and regions of difference (RD) virulence gene variation on tuberculosis recurrence
US20240068041A1 (en) Free dna-based disease prediction model and construction method therefor and application thereof
Liu et al. CRSCNV: A cross-model-based statistical approach to detect copy number variations in sequence data
CN110675917B (zh) 一种个体癌症样本的生物标记物识别方法
CN112562787B (zh) 一种基于ngs平台的基因大片段重排检测方法

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