CN109658981B

CN109658981B - 一种单细胞测序的数据分类方法

Info

Publication number: CN109658981B
Application number: CN201811501781.7A
Authority: CN
Inventors: 谢尚潜; 刘宇枭; 林加论; 邢剑锋
Original assignee: Hainan University
Current assignee: Hainan University
Priority date: 2018-12-10
Filing date: 2018-12-10
Publication date: 2022-10-04
Anticipated expiration: 2038-12-10
Also published as: CN109658981A

Abstract

本发明属于生物信息学分析技术领域，公开了一种单细胞测序的数据分类方法；包括有第一条序列(Read1.fastq)的信息识别模块；第二条序列(Read2.fastq)的信息识别模块；条形码列表信息加载(barcodeList)模块；引物信息加载(primerList)模块。本发明主要对单细胞测序技术SPLiT‑seq的数据进行分类，在分类过程中充分考虑条码barcode的信息，这是第一个对单细胞SPLiT‑seq技术的数据分类的方法。在barcode和特征序列中加入了容错比对机制，并使用碱基转换函数，将字符转换成数字进行运算，使单细胞测序数据的分类更加快速和高效。

Description

一种单细胞测序的数据分类方法

技术领域

本发明属于生物信息学技术领域，尤其涉及一种单细胞测序的数据分类方法。

背景技术

高通量测序技术(Next generation sequencing，NGS)是生命科学研究领域重要技术之一，近年来基于高通量测序技术的生命科学研究普遍应用于个体、组织等群体水平，如人的全基因组测序(Whole genome sequencing，WGS)、转录组测序(RNA sequencing，RNA-seq)等。由于多细胞的组织中广泛存在细胞异质性，即相同表型细胞的基因组、转录组等遗传信息可能存在显著性差异，因此有必要从单细胞水平对生物体进行解析和研究。虽然已有一些早期单细胞的研究方法，如免疫荧光，单细胞PCR和单细胞实时荧光定量PCR法，但这些方法普遍通量低，无法从单个细胞的水平上精确地解析细胞的分化、再生、衰老以及病变等复杂生物学现象。结合高通量测序技术的特点，单细胞测序技术便应运而生。单细胞测序(Single cell sequencing，SCS)是在单细胞水平对全基因组、转录组进行扩增和测序的一项新技术。现在称为年度6大领域最值得关注的技术，并且将其评为最重要的方法学研究成果。相对于常规的测序技术，单细胞测序技术具有明显优势，尤其对于组织水平研究中易丢失的低丰度的基因信息。单细胞测序技术能够揭示单个细胞的基因结构和基因表达状态，反映细胞间的异质性，在肿瘤、发育生物学、微生物学、神经科学等领域发挥重要作用。随着单细胞测序研究的深入，开发了多个用于单细胞测序的技术，比如DroNC-seq，CROP-seq，LIANTI，和scCOOL-seq等。DroNC-seq技术将低通量能有效观察单细胞类型和动态过程的sNuc-seq技术与利用微滴包括标记单细胞DNA的Drop-seq技术相结合，表现了明显的高效灵敏的细胞分类能力。CROP-seq又称为CRISPR液滴测序技术，结合CRISPR-Cas9在细胞筛选中的优势与单细胞液滴的方法，具有通量高、成本低的单细胞测序特点。LIANTI主要优势体现在单细胞的全基因组扩增方法上，通过转座子插入对DNA进行大量线性扩增，从而达到单细胞的DNA测序深度和分辨率的显著提高。scCOOL-seq是一种单细胞多重测序技术，可以完成对单细胞染色质状态、核小体定位、甲基化、拷贝数变异和染色体倍性等5个层面的特征解析。尽管上述单细胞测序技术各有优势和特点，但都需要微流体或者微孔做细胞分选，从而造成高昂的单细胞测序费用。现有技术中关于SPLiT-seq技术通过DNA的条码(barcode)的4轮扩增标记，可以获得潜在的21233664个条码组合，足以特异性的标记所有细胞，避免了单个细胞的分离。SPLiT-seq大大降低了单细胞测序的成本和实验设备要求，使单细胞测序技术的广泛推广和使用成为现实。由于SPLiT-seq不经过细胞分选进行测序，所有的细胞信息都存在了测序后的下机数据中。通过条码序列信息对所有的数据进行分类，以便确定相应的细胞和后续分析；针对SPLiT-seq单细胞测序的数据分类，由于该技术是2018年才提出的新技术，相应的数据分类的方法还未见报道。

综上所述，现有技术存在的问题是：针对SPLiT-seq单细胞测序的数据分类，目前还没有相应的方法。

解决上述技术问题的难度和意义：本发明的提出能有效的解决基于SPLiT-seq单细胞测序的数据分类和提取问题，为人类疾病尤其是肿瘤医学的基础科研工作提供重要的方法支持和帮助。

发明内容

针对现有技术存在的问题，本发明提供了一种单细胞测序的数据分类方法。

本发明是这样实现的，一种单细胞测序的数据分类系统，所述单细胞测序的数据分类系统包括：

第一条序列阅读模块和第二条序列阅读模块，作为实际测序数据的输入；

条形码列表信息模块，用于存放3轮条码序列信息，生成相应的表来加速查询过程；

引物信息加载模块，用于存放附件的引文序列信息，生成相应的表来加速查询过程。

进一步，所述第二条序列阅读模块2包括：UMI单元，3轮标签单元和cDNA单元；

UMI单元，作为标识，对不同的细胞来源进行分类处理；

3轮标签单元，标签作为标识，对不同的细胞来源进行分类处理；

cDNA单元，最终要提取的序列信息。

本发明的另一目的在于提供一种运行所述单细胞测序的数据分类系统的单细胞测序的数据分类方法，所述单细胞测序的数据分类方法包括：

步骤一，对实际测序数据分类提取；

步骤二，存放3轮条码序列信息；生成相应的表来加速查询过程；

步骤三，存放附件的引文序列信息；生成相应的表来加速查询过程。

进一步，所述步骤一具体包括：

(1)采用K-mer方法，提供容错机制，查找到特征序列的位置；

(2)从序列中提取出3轮的barcode：根据序列中特征序列的位置，向前偏移8位，提取出相对应的barcode条码；

(3)提取出3轮barcode后，通过Barcode Table将barcode转换为3组数字，共同作为唯一标识，确定一个细胞；

(4)将UMI追加到标识之后。

进一步，所述步骤一具体包括：

(1)对primerList进行建表，产生primer table；

(2)将primerList中的整条长序列，每次对其取长度为k的片段，从起始处开始取片段，每次向后偏移1位，记录子序列在整体中出现的位置；

(3)对每个fragment进行一次转换，将其视为一个4进制数，得到碱基到数字的具体映射关系；

(4)将映射关系通过十进制转换，获得的数值作为数组的下标，得到最终的primertable结构。

本发明的另一目的在于提供一种实现所述单细胞测序的数据分类方法的计算机程序。

本发明的另一目的在于提供一种实现所述单细胞测序的数据分类方法的信息数据处理终端。

本发明的另一目的在于提供一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行所述单细胞测序的数据分类方法。

综上所述，本发明的优点及积极效果为：本发明对单细胞测序技术SPLiT-seq的数据进行分类，在分类过程中充分考虑条码barcode的信息，是第一个对单细胞SPLiT-seq技术的数据分类的方法。本发明在barcode和特征序列中融入了容错比对机制，并使用碱基字母ATCG的转换函数，使单细胞测序数据的分类更加快速和高效。

附图说明

图1是本发明实施例提供的单细胞测序的数据分类系统结构示意图；

图中：1、第一条序列阅读模块；2、第二条序列阅读模块第；3、条形码列表信息模块；4、引物信息加载模块。

图2是本发明实施例提供的单细胞测序的数据分类方法流程图。

图3是本发明实施例提供的Read2数据内容示意图。

图4是本发明实施例提供的primer table的生成示意图。

图5是本发明实施例提供的primer table的结构示意图。

图6是本发明实施例提供的根据BarcodeList的三轮信息生成3张表示意图。

图7是本发明实施例提供的根据PrimerList生成PrimerTable示意图。

图8是本发明实施例提供的最终primertable的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对SPLiT-seq单细胞测序的数据分类，目前还没有相应的方法。本发明在barcode和特征序列中融入了容错比对机制，并使用碱基字母ATCG的转换函数，使单细胞测序数据的分类更加快速和高效。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例提供的单细胞测序的数据分类系统包括：

第一条序列阅读模块1和第二条序列阅读模块2，作为实际测序数据的输入；

条形码列表信息模块3，用于存放3轮条码序列信息，生成相应的表来加速查询过程；

引物信息加载模块4，用于存放附件的引文序列信息，生成相应的表来加速查询过程。

在本发明的优选实施例中，本发明实施例提供的第二条序列阅读模块2包括：UMI单元，3轮标签单元和cDNA单元；

UMI单元，作为标识，对不同的细胞来源进行分类处理；

cDNA单元，最终要提取的序列信息。

如图2所示，本发明实施例提供的单细胞测序的数据分类方法包括：

S201：对实际测序数据分类提取；

S202：存放3轮条码序列信息；生成相应的表来加速查询过程；

S203：存放附件的引文序列信息；生成相应的表来加速查询过程。

在本发明的优选实施例中，步骤S201具体包括：

(1)采用K-mer方法，提供一位容错机制，查找到特征序列的位置；

(3)提取出3轮barcode后，通过Barcode Table将barcode转换为3组数字，共同作为唯一标识，来确定一个细胞；

(4)将UMI追加到标识之后。

在本发明的优选实施例中，步骤S201具体包括：(1)对primerList进行建表，产生primer table；

下面结合具体实施例对本发明的应用原理作进一步的描述。

实施例1；

1、Read2.fastq数据的分类提取

如图3所示，为Read2数据内容，Read2一共分为5个部分，UMI，3轮标签和cDNA，其中UMI和3轮标签作为标识来对不同的细胞来源进行分类，cDNA则是最终要提取的序列信息。

(1)首先从序列中提取出3轮的barcode，具体方法为先找出序列中特征序列的位置，之后向前偏移8位便可提取出相对应的barcode条码。在查找特征序列的位置时，采用的是K-mer方法，并提供一位容错的机制。提取出3轮barcode之后，通过Barcode Table将barcode转换为3组数字，共同作为唯一标识，来确定一个细胞。然后将UMI追加到标识之后。

获取到特征序列之后，便可提取出barcode。以上述Read2为例，提取出的barcode则为：

·Barcode1:CCATCCTC；

·Barcode2:GCGAGTAA；

·Barcode3:AATGTTGC；

之后对barcode做一次修正(一位容错)。

假设Barcode1的值为ACATCCTC，第一位在测序的时候发生了错误，错把C测成了A。本发明会将Barcode1在barcodeList中进行比对，通过距离函数找出一个与Barcode1最相似的一个barcode修正Barcode1，如果得出的距离(差异性)大于等于阈值(阈值为2)时，视为该barcode无效，丢弃整个Read2。如果在阈值范围内，则经过修正后的Barcode1就变成了CCATCCTC。

(2)结果：

真实Read2数据情况：

@A00268:72:H52W3DSXX:2:1101:2085:10002:N:0:ACTTGAA

CCTTGTTGCGAATGTTGCGTGGCCGATGTTTCGCATCGGCGTACGACTGCGAGTAAATCCACGTGCTTGAGAGGCCAGAGCATTCGCCATCCTCGCCCCTTGATCCCAACTCATCTCTCATTTATTTCGGCTTCTTTTATTCCAGGATTA+

FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF

分类提取后的Read2为：

@85_51_73:CCTTGTTGCG A00268:72:H52W3DSXX:2:1101:2085:1000 2:N:0:ACTTGAA

GCCCCTTGATCCCAACTCATCTCTCATTTATTTCGGCTTCTTTTATTCCAGGATTA

+

FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF

第一行为分类信息(855173)+UMI(CCTTGTTGCG)+原有测序信息(A00268:72:H52W3DSXX:2:1101:2085:1000 2:N:0:ACTTGAA)

第二行为提取出来的cDNA信息：(GCCCCTTGATCCCAACTCATCTCTCATTTATTTCGGCTTCTTTTATTCCAGGATTA)

第三行保留，第四行根据cDNA做相应的碱基质量信息。

2、R1的分类与提取

在R2中的序列有效的情况下才会对R1中对应的序列进行分类和提取，如果R2无效。则直接丢弃与之相对应的R1序列。关于R1的分类与提取只需要保留中间的有效部分，剔除首尾的无效部分即可。而如何确定哪些是"有效部分"是根据primerList中序列来确定的，只有出现在primerList中的序列(正序逆序均可)才作为有效正确序列。

首先是对primerList进行建表产生primer table。

如图4所示，primer table的生成是根据输入的文本文件primerList来生成的一个链表数组。将primerList中的数据视为一整条长序列，每次对其取长度为k的片段，从起始处开始取片段，每次向后偏移1位，主要是记录子序列在整体中出现的位置。

其中对每个fragment进行一次转换，将其视为一个4进制数，碱基到数字的具体映射关系如下:

'A'->0；

'T'->1；

'G'->2；

'C'->3；

以上述1中的fragment AGTCGTAC为例，转换完成后的四进制数为02132103再将其转换为十进制为10131；

然后十进制转换后的数值作为数组的下标，其出现的位置为0，添加到链表10131的尾端，最终primer table的结构如图5所示；

如图5所示，最终primer table的结构示意图。

因此，通过下标直接查询，时间复杂度为O(1)，从而实现快速查询过程。得出有效部分的起始位点和结束位点之后，只需要直接切割便可完成R1的分类和提取。

实施例2：

步骤一，加载实际数据以及相关文件：

实际数据文件2个：

R1.fastq

R2.fastq

三轮barcode信息文件：

BarcodeList

特征信息：

PrimerList

步骤二，根据BarcodeList和PrimerList生成相应的表来加速查询过程：根据BarcodeList的三轮信息生成3张表，如图6：

根据PrimerList生成PrimerTable如图7：

PrimerTable的生成是根据输入的文本文件PrimerList来生成的一个链表数组。将PrimerList中的数据视为一整条长序列，每次对其取长度为k的片段，从起始处开始取片段，每次向后偏移1位，主要是记录子序列在整体中出现的位置。

其中对每个fragment进行一次转换,将其视为一个4进制数,碱基到数字的具体映射关系如下:

'A'->0

'T'->1

'G'->2

'C'->3

以上述1中的fragment AGTCGTAC为例,转换完成后的四进制数为02132103再将其转换为十进制为10131

然后intSeq作为数组的下标，其出现的位置为0，添加到链表10131的尾端，

seq＝AGTCGTAC

n＝len(seq)

最终primer table的结构如图8。

这样就可以通过下标来直接查询，以这样一种用空间换时间的方式来加速查询过程。

得出有效部分的起始位点和结束位点之后，只需要直接切割便可完成R1的清洗与变换。

步骤三，读取R2的一条数据，将其清洗与变换：

以下列Read2为例：

Read2一共分为5个部分，UMI，3轮标签和cDNA，其中UMI和3轮标签作为标识来区分不同的细胞，cDNA则是最终要取得的部分。首先要做的是从序列中提取出3轮的barcode，具体方法为先找出序列中特征序列1的位置，之后向前(或向后)偏移8位便可提取出相对应的barcode。在查找特征序列的位置时，采用的是K-mer2方法，并提供一位容错的机制。

Barcode1:CCATCCTC

Barcode2:GCGAGTAA

Barcode3:AATGTTGC

之后会对barcode做一次修正(一位容错)。这里假设Barcode1的值为ACATCCTC，第一位在测序的时候发生了错误，错把C测成了A。程序会拿着Barcode1在barcodeList中进行比对，通过距离函数3找出一个与Barcode1最相似的一个barcode来修正Barcode1，如果得出的距离(差异性)大于等于阈值(阈值为2)时，视为该barcode无效，丢弃整个Read2。如果在阈值范围内，则经过修正后的Barcode1就变成了CCATCCTC。

提取出3轮barcode之后，通过Barcode Table将barcode转换为3组数字，共同作为唯一标识，来确定一个细胞。然后将UMI追加到标识之后。

步骤四，读取R1的一条数据，将其清洗与变换：

只有在R2中的序列有效的情况下才会对R1中对应的序列进行清洗与变换，如果R2无效。则直接丢弃与之相对应的R1序列。关于R1的清洗与变换只需要保留中间的有效部分，剔除首尾的无效部分即可。而如何确定哪些是"有效部分"是根据步骤一中加载的PrimerList来确定的，只有出现在PrimerList中的序列(正序逆序均可)才视为“有效部分”。

成功读取到一条Read1之后，将其分割为n个片段，每段长度为k，并将其通过步骤一中所述的方法转换为二进制数，结合步骤一中生成的PrimerTable来进行查询判断有效部分的具体位置。

步骤五，将数据输出到指定文件。

步骤六，重复步骤三，四，五，直到处理完所有的数据为止。

在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用全部或部分地以计算机程序产品的形式实现，所述计算机程序产品包括一个或多个计算机指令。在计算机上加载或执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(DSL)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输)。所述计算机可读取存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，DVD)、或者半导体介质(例如固态硬盘SolidState Disk(SSD))等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims

1.一种单细胞测序数据的分类系统，其特征在于，所述单细胞测序的数据分类系统包括：

所述第二条序列阅读模块用于实现：(1)采用K-mer方法，提供容错机制，查找到特征序列的位置；

(3)提取出3轮barcode后，通过Barcode Table将barcode转换为3组数字，共同作为唯一标识，确定数据来源；

(4)将UMI追加到标识之后；

所述第一条序列阅读模块用于实现：只有在第二条序列有效的情况下对第一条序列中对应的序列进行阅读，具体包括：

(1)对primerList进行建表，产生primertable；

(4)将映射关系通过十进制转换，获得的数值作为数组的下标，得到最终的primertable结构；

2.如权利要求1所述的单细胞测序数据的分类系统，其特征在于，所述第二条序列阅读模块2包括：UMI单元，3轮标签单元和cDNA单元；

UMI单元，作为标识，对不同的细胞来源进行分类处理；

cDNA单元，最终要提取的序列信息。

3.一种运行权利要求1所述单细胞测序数据的分类系统的单细胞测序的数据分类方法，其特征在于，所述单细胞测序的数据分类方法包括：

步骤一，将实际测序数据输入，包括第一测序数据和第二条测序数据；

4.一种实现权利要求3所述单细胞测序的数据分类方法的信息数据处理终端。

5.一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行如权利要求3所述单细胞测序的数据分类方法。