阅读说明：本技术 数据处理方法 (Data processing method ) 是由 杨圳 王文山 于 2019-08-27 设计创作，主要内容包括：本发明揭示了一种数据处理方法,具体来说涉及ICB-scSeq技术中的下机数据的处理方法。使用该数据处理方法,能够低成本地、简单地、快速地进行细胞测序,具有相当大的经济效益和安全效益。(Present invention discloses a kind of data processing method, the processing method of the lower machine data in ICB-scSeq technology is related in particular to.Using the data processing method, can at low cost, simply, rapidly carry out cell sequencing, there is sizable economic benefit and safety benefit.)

数据处理方法

技术领域

本发明涉及一种数据处理方法。具体来说，涉及ICB-scSeq(IntelligentCombinatorial Barcoding-single cell Sequencing，智能组合条形码法单细胞测序)技术中的下机数据的处理方法。

背景技术

在过去十年里，随着第二代测序(next generation sequencing,NGS)技术和第三代测序(third generation sequencing,TGS)技术的飞速发展，引起生命科学领域的巨大变革。以往的研究，需要从大量细胞中获得足够多的核酸进行测序，因此测序结果往往表示的是细胞群体的表征，而单个细胞独有的细胞特性往往被忽略。为了解决上述局限性问题，单细胞测序技术应运而生。

单细胞测序已经在肿瘤、发育生物学、神经科学等领域取得了丰富的成果。而单细胞的研究本可以更为迅速的扩大科学成果，但单细胞测序技术还存在着许多问题，比如需采用新鲜的细胞，样品利用率不高，昂贵的设备及相关试剂等，这对单细胞测序技术的研究和推广带来了很多不便，而对于单细胞生命科学研究的广泛开展也有诸多不利。因此优化开发新的单细胞测序技术，就显得刻不容缓。

ICB-scSeq(Intelligent Combinatorial Barcoding-single cell Sequencing，智能组合条形码法单细胞测序)是由本发明的发明人们研发的单细胞测序技术，它是一种基于SPLIT(split-pool ligation-based transcriptome，基于分割池连接的转录组测序)技术通过组合条形码(combinatorial barcoding)标记RNA细胞起源的单细胞测序的方法。

因此，由于存在上述技术缺陷，在ICB-scSeq的测序方法中，也需要找到一种更好的从原始下机数据到下游分析的数据处理方法，能够低成本地、简单地、快速地进行细胞测序。

发明内容

本发明的一个目的在于克服现有技术的缺陷，提供一种能够低成本地、简单地、快速地进行细胞测序的数据处理方法。

为实现上述目的，本发明提出如下技术方案：一种数据处理方法，其特征在于，包括：

原始数据获取步骤，进行双端测序而对用于智能组合条形码法单细胞测序的原始数据进行获取，第一端为cDNA部分，第二端为特异性分子标签以及细胞条形码部分；

质量控制以及过滤步骤，对所获取的原始数据进行过滤而得到过滤后数据；

比对步骤，将过滤后数据与参考基因组序列进行比对而得到比对后数据；

特异性分子标签去重步骤，对比对后数据中的特异性分子标签重复的部分进行去除而得到去重后数据；

基因定量步骤，对去重后数据进行基因定量而得到定量后数据；

表达矩阵构建步骤，根据定量后数据构建表达矩阵，该表达矩阵包括了每个细胞中的每个基因的原始计数值；

细胞筛选步骤，对表达矩阵的线粒体含量和表达基因数目进行筛选而得到筛选后矩阵；

标准化处理步骤，对筛选后矩阵的原始计数值进行标准化而得到标准化矩阵；

分析步骤，对标准化矩阵进行分析。

根据本发明所提供的数据处理方法，能够低成本地、简单地、快速地进行细胞测序，具有相当大的经济效益和安全效益。

附图说明

图1是本发明的第一实施方式的数据处理方法的示意图。

图2是图1的数据处理方法中所用的基因序列的示意图。

图3是图1的数据处理方法中的去重过程的示意图。

图4是图1的数据处理方法中的聚类分析的成果的展示图。

图5是图1的数据处理方法中的通路富集分析的成果的展示图。

图6是图1的数据处理方法中的通路富集分析的成果的另一个展示图。

具体实施方式

下面将结合本发明的附图，对本发明实施例的技术方案进行清楚、完整的描述。

本发明的第一实施方式是一种数据处理方法。

图1是本发明的第一实施方式的数据处理方法的示意图。如图1所示，该数据处理方法包括：原始数据获取步骤、质量控制以及过滤步骤、比对步骤、特异性分子标签去重步骤、基因定量步骤、表达矩阵构建步骤、细胞筛选步骤、标准化处理步骤、分析步骤。

在原始数据获取步骤中，进行双端测序(paired-end sequencing)而对用于智能组合条形码法单细胞测序(ICB-scSeq)的原始数据进行获取，第一端即read1端为cDNA部分，第二端即read2端为特异性分子标签以及细胞条形码(UMI+cell barcode)部分。cDNA是指具有与某RNA链呈互补碱基序列的DNA。UMI(Unique Molecular indentifier)是特异性分子标签。

在质量控制以及过滤步骤中，对所获取的原始数据进行过滤而得到过滤后数据。在本实施方式中，示例性地示出了质量控制以及过滤步骤包括如下子步骤：对所获取的原始数据的第二端的细胞条形码部分进行校正；构建细胞条形码的白名单；根据白名单提取第一端的序列；对所提取的第一端的序列进行筛选以进行过滤而得到过滤后数据。但是本发明并不限于此，质量控制以及过滤步骤也可以包括其他子步骤。

具体来说，在每一条read2中有三段cell barcode，分别是BC1,BC2,BC3,每段都是8bp的长度(如图2所示)。并且这些barcode的序列每次都是固定的。比如，如果barcode1使用了96种组合，那么说明barcode1的序列总共只有96种，每个都是8bp。因此根据hammingdistance(汉明距离)等于1的校正原则来对每条read进行校正。

在每一条read中，提取出BC1,BC2,BC3所在位置的三段8bp序列作为候选barcode序列(标记为barcode1-new,barcode2-new,barcode3-new)。然后依次对barcode1-new与已经确定好的barcode1的列表里的所有序列进行比对，计算hamming distance,记为hd。如果hd等于0，则不进行变动，如果hd等于1，将该barcode1-new的序列更改为相对应barcode1的序列。以完成barcode序列的校正过程。

在完成cell barcode的校正后，根据预估的细胞数目对cell barcode序列进行合并作为一个细胞的唯一标识(cell UID)，构建一个cell barcode的白名单。在这个白名单的列表里面，是所有能够鉴定到的细胞的UID。

根据上一步建好的cell barcode的白名单对read1端里面的cDNA序列进行提取。对任意一条read1的序列，如果它所对应的read2里面的cell UID在cell barcode白名单里面，那么这条read1就会被提取出来。构建白名单和提取序列可以使用开源工具umi-tools进行处理。

在提取完read1序列后，还需要对序列进一步筛选，主要是去除末尾的polyA结构(如下式所示)，以及序列两端的低质量值。下式中，上面一行为原始序列，下面一行为去除了末尾的polyA结构以及序列两端的低质量值的序列。

在比对步骤中，根据上述筛选得到的read1序列，与参考基因组的序列进行比对，该比对可以使用比对软件STAR来进行。

根据比对的结果，得到了一个已经经过排序的bam文件。对每一条比对上的序列根据参考基因组的GTF文件进行注释，也就是进行基因的指定。目的是明确每一条比对上的序列经过GTF文件的注释后属于哪个基因。该指定可以使用开源工具featureCounts来完成。

在特异性分子标签去重步骤中，根据上一步的结果，已经得知了每一条比对上的read属于哪一个基因。因为ICB-scSeq为了消除样本后续建库时的PCR-bias而在每条序列中引入了一段10bp长的UMI序列。这样，如果在同一个基因范围内，出现了两条相同的序列并且如果序列的UMI的10bp也是相同的话，那么我们认为这两条read来自于同一条cDNA分子，是需要去重的。如图3所示，图3的左侧示出了五条read，但是在这五条read中，上面三条read是重复的，下面两条read也是重复的，因此在去重之后，右侧的read只有两条。

在基因定量步骤中，对去重后数据进行基因定量而得到定量后数据。

在表达矩阵构建步骤中，根据定量后数据构建表达矩阵，该表达矩阵包括了每个细胞中的每个基因的原始计数值(raw counts)。在这个矩阵中，每一列代表一个细胞的UID,每一行代表一个基因的ID，如下表所示。

在细胞筛选步骤中，对表达矩阵的线粒体含量和表达基因数目进行筛选而得到筛选后矩阵。具体来说，对表达矩阵里面的每个细胞的数据进行计算，计算线粒体基因的所有表达值的比例，如果这个比例超过设定的阈值，那么就筛选掉这个细胞。阈值例如为5％，但是并不限于此，也可以设为其他阈值。另外还需要对表达矩阵中每个细胞中表达基因的数量进行筛选，一般的筛选标准例如是最低表达的数目是200，最高表达的数目是2500，但是并不限于此，也可以设为其他范围。筛选步骤可以使用Seurat来进行。经过两次筛选，得到了一个筛选后的表达矩阵，可以进行下一步的处理。

在标准化处理步骤中，对筛选后的表达矩阵的原始计数值进行标准化而得到标准化矩阵。由于单细胞测序过程中，每个细胞测得到reads的数目不均一，为了能够消除因为测序深度导致的定量误差，需要对raw counts进行标准化处理。标准化处理步骤可以使用Seurat来进行，其标准化计算式如下：

其中CountOfGene代表每个细胞中每个基因的raw counts，AllCount代表每个细胞中所有基因的raw counts的和。

在分析步骤中，对标准化矩阵进行分析。在本实施方式中，分析步骤在细胞层面上为聚类分析步骤，分析步骤在基因层面上为差异分析步骤以及通路富集分析步骤，但是并不限于此，也可以是其他合适的分析步骤。

聚类分析的分析方法如下。

首先进行特征提取，对所有测到的单细胞进行聚类分簇分析。先对标准化后的表达矩阵计算出高变异的特征，把这些特征提取出来进行后续的分析。

然后对矩阵数据进行尺度化分析，为了尽可能去除掉一些数据来源的误差(包括技术误差，批次误差以及一些生物学来源的误差)，对矩阵数据进行回归处理，排除掉这些误差，以提高后续降维和聚类的效果。

然后进行线性降维分析，利用PCA(principal component analysis,主成分分析)的方法对已经过尺度化分析的数据进行降维分析。

然后进行聚类分群分析，根据上一步鉴定出的显著性的PC(主成分),使用基于图形的聚类方法。该方法根据构建的KNN(K-nearest neighbor,K最近邻)图形以及Louvain算法来迭代地进行聚类，最终把所有的细胞聚到不同的类群里面。以上的分析过程可以使用Seurat来进行。

最后进行UMAP二维展示，如图4所示，根据上一步聚类的结果，使用UMAP(uniformmanifold approximation and projection,均匀流形近似和投影)的方法进行二维展示。该展示方法可以使用Seurat来进行分析。

聚类分析的分析方法如下。

根据聚类的结果，对所有的cluster使用Wilcoxon rank sum test(威尔科克森秩和检验)的方法进行差异基因筛选分析，得到一个关于所有cluster的差异表达基因的列表。结果如下表所示。

通路富集分析的分析方法如下。

第一种通路富集分析为GO(Gene Ontology,基因本体论)富集分析。如图5所示，根据上一步差异基因的结果，对每一个cluster的差异基因进行GO富集分析。

除了GO富集分析，还能进行KEGG(Kyoto Encyclopedia of Genes and Genomes,京都基因与基因组百科全书)通路富集分析，如图6所示，鉴定出在每个cluster里面差异基因显著性地富集到哪些通路里面，使用了气泡图进行展示。

如上所述，利用第一实施方式的数据处理方法，能够低成本地、简单地、快速地进行细胞测序，具有相当大的经济效益和安全效益。

需要说明的是，本发明各设备实施方式中提到的各单元都是逻辑单元，在物理上，一个逻辑单元可以是一个物理单元，也可以是一个物理单元的一部分，还可以以多个物理单元的组合实现，这些逻辑单元本身的物理实现方式并不是最重要的，这些逻辑单元所实现的功能的组合才是解决本发明所提出的技术问题的关键。此外，为了突出本发明的创新部分，本发明上述各设备实施方式并没有将与解决本发明所提出的技术问题关系不太密切的单元引入，这并不表明上述设备实施方式并不存在其它的单元。

需要说明的是，在本专利的权利要求和说明书中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

虽然通过参照本发明的某些优选实施方式，已经对本发明进行了图示和描述，但本领域的普通技术人员应该明白，可以在形式上和细节上对其作各种改变，而不偏离本发明的精神和范围。