CN107463801A - 一种Drop‑seq数据质量控制和分析方法 - Google Patents

Abstract

本发明公开了一种Drop‑seq数据质量控制和分析方法，属于生物学技术领域。其包括以下步骤：第一步、使用两个配对的测序文件作为输入文件；第二步、对测序文件进行处理；第三步、测序片段的层面、多细胞的层面、单个细胞的层面以及细胞聚类的层面对Drop‑seq数据进行质量控制流程的开发。第四步、给用户提供一份报告文档。本发明以Drop‑seq数据研究作为出发点，主要涵盖测序片段层面的质量控制、多细胞层面的质量控制以及单个细胞层面的质量控制，借助Drop‑seq数据对单个细胞的标识，针对细胞异质性的特点与来源进行生物信息方法开发，探究细胞异质性对细胞基因表达与表观特征的影响。

Description

一种Drop-seq数据质量控制和分析方法

技术领域

本发明属于生物学技术领域，特别涉及基于Drop-seq数据特征的生物信息分析技术领域，开展基于Drop-seq数据的质量控制与数据分析流程的探究和开发。

背景技术

随着第二代测序（Next-generation sequencing）技术的飞速发展，逐步增长的通量和逐年降低的测序成本使得人们利用测序技术对转录组学与表观遗传组学的研究逐渐深入，并开发了一系列定量研究手段如用于检测特定条件下各基因的表达水平的转录组测序技术（RNA-seq）以及定性研究手段如用于观测DNA与转录因子蛋白相互作用的染色质免疫共沉淀技术（ChIP-seq）。随着第二代测序的通量和分辨率不断提高，一些其它方面的问题也逐渐出现，为了解决细胞异质性以及稀有细胞无法检测的问题，近年来开发了基于单细胞的测序技术。scRNA-seq技术，可以获得单个细胞的全基因组表达谱，从而使scRNA-seq技术成为同时研究上千个细胞之间的表达差异的的工具。伴随着scRNA-seq技术飞速的发展，与此同时其数据分析的方法也越来越成熟。单细胞测序数据分析的生物信息学方法也与传统的多细胞测序数据的分析方法有着较大的差别。单细胞测序数据由于其细胞用量少的原因使得测序中会丢失大量的序列片段。与基于多细胞的第二代测序方法相比，单细胞测序数据具有零值多，精度低等特点，为第二代测序所开发的众多生物信息学算法以及工具都需要进一步的考量与改进，甚至是重新开发以适应单细胞测序数据。

获得单个细胞之间的转录本的差异，是了解复杂组织和功能响应的一个重要方法，但是由于单细胞测序对于同时进行多个单细胞的文库制备所花费的时间和成本的限制，以及单细胞测序数据分析算法的缺少，限制了人们的研究。而基于液滴微流控技术的转录组测序方法（Drop-seq）可以在一次实验中得到大批量的单个细胞的基因表达谱数据，有效的弥补了传统单细胞测序技术的不足。Drop-seq是通过把细胞封装在微小液滴中进行成千上万个单细胞的mRNA并行测序，Drop-seq技术的优势使其可以在广泛的领域中使用。例如，全基因组范围的遗传学研究主要目的是确定哪些基因的改变会对疾病的风险有影响，但是目前生物学上还缺乏类似把这些基因关联到特定疾病和某一功能响应上的高通量方法，Drop-seq技术可以用来探索这些基因是如何作用在各个组织的不同类型细胞的机制。另外，Drop-seq可以得到突变体、病原体或者其它刺激源的扰动信息，进而获得其在许多细胞类型上的包含扰动影响的一个信息丰富、多维的测序结果。基于液滴微流控分选细胞技术还可以与其它的表观遗传组学研究手段相结合从而实现了大批量的单细胞表观遗传组测序技术。

随着细胞分选技术的发展（例如液滴微流控技术），大批量单细胞测序技术（例如Drop-seq）将会越来越多，随之意外的、深远的发现将在短时期内大量涌现出来，包括对新的细胞亚型的鉴定，识别基因表达模式来预测细胞的状态以及研究随机转录对功能的影响。但是这些结论都是建立在稳定的计算方法的基础上。这些技术目前普遍存在目标片段的捕获率低、转换效率低以及测序深度低的问题，使得测序结果不足以覆盖每个细胞中的大部分信息，并且相比多细胞测序和传统单细胞测序增加了测序的噪音水平。因此，大批量单细胞测序数据的好坏与否对于下游分析和研究有极大的影响。质量控制应该是在分析流程中的第一步骤，以确保数据的质量，并为随后的分析提供一个坚实基础。目前尚无针对Drop-seq数据的质量控制算法，Drop-seq数据的低mRNA捕获率、超高维、低精度的特征，开发质量控制流程是极其必要的，针对Drop-seq这样一个全新的技术手段及独特的数据特点，开发生物信息学质量控制和分析方法也迫在眉睫。

发明内容

本发明的目的开发针对基于液滴微流控技术的转录组测序数据的质量控制体系确，保数据可信度。

为达到上述目的，本发明的具体技术方案为：

一种Drop-seq数据质量控制和分析方法，该方法包括以下步骤：

第一步、使用两个配对的测序文件作为输入文件，其中一个文件包含转录本的信息，另一个文件包含细胞条形码和UMI的信息。即Drop-seq数据的结构由细胞条形码、UMI和cDNA三部分组成。其中细胞条形码是用来标识不同的细胞，由12个碱基组成；UMI，由8个碱基组成，用来识别同一个细胞中不同的cDNA。

第二步，对测序文件进行处理：对所得的原始测序数据：去除细胞条形码碱基质量小于10的测序片段对，第二个测序片段（50bp）在5’端去除TSO的接头片段，3’端移除长度为6个碱基或者更长的多聚腺苷酸尾（polyA tails）。在测序序列比对中，STAR被定为默认的比对工具，使用STAR的默认设置参数，并且把Bowtie2作为可选择比对工具。此外，为了进一步获得更高质量的比对结果，使用Samtools移除测序片段的测序质量小于30（Q30）的测序片段。

第三步，把Drop-seq数据看作是多细胞RNA-seq数据进行分析。在Drop-seq一次测序产生的上千万个测序片段的数据，随机抽样获得5百万个测序片段进行质量控制分析，使用软件包RseQC进行计算，提供以下层面的质量控制的测量；

（1）测序片段层面的质量控制包括测序片段的质量、核苷酸的组成以及测序片段的GC含量。具体为对测序质量分布、每个测序位点的碱基组成以及每个测序片段的GC含量这几个测度进行计算。

（2）多细胞层面的质量控制包括基因区域的覆盖程度和序列回帖率；具体方法首先，列出测序片段的比对结果总结，包括序列的回帖率和全基因组范围的序列分布，并通过绘制基因区域的序列覆盖程度来估计5’端和3’端的测序偏差。然后生成表达矩阵：使用bedtools工具整合测序片段、细胞条形码以及基因组的注释信息；进行测序片段的合并，如果具有相同细胞条形码的测序片段位于基因组同一个位置并且具有相同的UMI，代表这些测序片段是由同一个片段经过PCR扩增而来的，因此这些冗余测序片段将被去除；把合并后去重的测序片段通过基因的注释信息和细胞条形码信息生成基因表达矩阵，并且利用不同细胞条形码信息把测序片段分配到不同的细胞中。

（3）单个细胞层面的质量控制包括，挑选出单个细胞中连接到微磁珠上的转录组，即STAMPs，在（2）中获得的表达矩阵是包含所有细胞的表达矩阵，即同时包括STAMPs和“empty”细胞（没有有效信息的细胞），从混有“empty”细胞的表达矩阵中区分出STAMPs。根据细胞条形码信息把测序片段分组到不同的细胞条形码中。

首先，以单个的细胞为研究对象，利用UMI和每个测序片段的基因组上位置的共同信息，计算每个细胞的独特测序片段数目、测序片段的重复率以及覆盖的基因数目。

其次，筛选STAMPs条形码：挑选的细胞条形码需要满足其覆盖的基因数目大于用户定义的阈值（1000为设定的默认参数值）。当细胞条形码的测序片段重复率太低时（0.1为设定的默认参数值）将在这一步筛选中被舍弃。

最后，评估在单个细胞的层面上mRNA的捕获效率：通过覆盖基因数目的分布和属于内含子测序片段的比例（这个比例定义为测序片段回帖到内含子和基因间区的比）。

第四步、给用户提供一份质量控制的报告文档，其中描述了第三步中质量控制的测量值。

优选地，

所述第三步中的质量控制细胞聚类层面，具体方法如下：

a.在区分出STAMPs的基础上，选择STAMPs条形码的方法产生相对应的表达矩阵，这个表达矩阵以STAMPs条形码为列基因为行，并且表达矩阵中的每个值是由原表达值进行log转换的值。

b.计算表达矩阵中每个基因在所有单细胞中表达的均值和离差值（方差除以均值），并且根据基因的平均表达值把所有的基因分配到20个bin中。

c.对每个bin中的所有基因的离差进行Z-标准化，定义z-score的阈值为1.7来确定高可变基因，这些高可变基因将作主成分分析（PCA）的输入，以确保准确的识别数据的主要结构。

d.在获得的结果中挑选最主要的主成分，这些主成分可以解释超过50%的基因方差。

e.使用t-SNE（t-分布随机近邻嵌入）降维方法进一步对这些主成分进行降维到2维层面。

f.对t-SNE的输出结果进行K均值（k-meams）聚类，并且在聚类中运用gap统计（gapstatistics）确定类的数目，并将STAMPs条形码分配到不同的类中。

g.利用gap statistics、silhouette score和t-SNE的可视化输出结果来评估细胞聚类的好坏以及样本的异质性。

其中，使用gap统计确定K-means聚类的K的具体方法如下:

进行gap statistic，当用户输入一个最大的候选k值（记为kmax）后，计算gap（k_max+1）, gap（k_max+2）, …, gap（k_max+20）。计算这20个gap得分的平均值和方差，提供一个稳定的gap得分的参考值。当首次ks满足下面的两个条件时，确定该ks为k均值聚类的k。

这里，代表在给定下k均值聚类的gap score。

或者，可选择基于密度的聚类算法DBSCAN（Density-Based Clustering ofApplications with Noise，DBSCAN），它将簇定义为密度相连点的最大集合，能够把具有足够高密度的区域划分为簇，并且可以在噪声的空间数据库中识别任意形状的聚类。

本发明以Drop-seq数据研究作为出发点，主要涵盖测序片段层面的质量控制、多细胞层面的质量控制以及单个细胞层面的质量控制。借助Drop-seq数据对单个细胞的标识，针对细胞异质性的特点与来源进行生物信息方法开发，探究细胞异质性对细胞基因表达与表观特征的影响。最终用户将会得到测序样本的异质性分析结果。

附图说明

图1为本发明的Drop-seq数据质量控制与数据分析流程控制方法的步骤流程图；

图2小鼠视神经细胞的Drop-seq 数据的数据和参数描述；

图3小鼠视神经细胞的Drop-seq 数据的基因测序片段覆盖程度累积图；

图4小鼠视神经细胞的Drop-seq 数据的覆盖基因的数目柱状图。

具体实施方式

为了能够更清楚地理解本发明的技术内容，特举以下实施例详细说明。

下面参照附图用本发明的示例性实施例对本发明进行更全面的描述及说明，但并不意味着本发明仅限于此。

如图1 所示，为本发明的Drop-seq数据质量控制与数据分析流程控制方法的步骤流程图。

在实施方式中，该质量控制和分析方法，使用小鼠视神经细胞的Drop-seq 数据（GSE1626793，总共包含510,210,716个reads）作为示例（如图1 所示），包括以下步骤：

第一步、使用两个配对的测序文件（FASTQ或者SAM格式）作为输入文件。这两个文件其中一个是包含转录本的信息（read1），另外一个文件则包含细胞条形码和UMI的信息。

第二步、处理测序文件：把转录本的测序文件比对到参考基因组上，并且只留取具有高测序质量的比对上的测序片段。

第三步、生成该数据的基本描述和软件参数表，如图2所示。

（1）测序片段层面的质量控制包括测序片段的质量、核苷酸的组成以及测序片段的GC含量；

a.测序片段层面的质量控制报告：针对测序片段的质量检测，衡量测序质量的好坏。从总的测序片段中随机抽样5百万个片段进行分析。

b.测序片段的质量：用测序片段上每个位置的Phred Quality Score绘制一个广泛的分布图来测量用户提供的数据的基本测序质量。可以观察到在近3’端出现质量的下降现象，这是因为持续长时间的测序会导致质量的降低。如果测序质量的降低影响了回帖率，常见的补救措施是对测序片段进行截断，可以根据测序片段的平均质量进行截取，或者可以直接在近3’端进行截断。如果测序数据处理后对回帖率还是没有改善，或许用户应该考虑是否是所采用的Drop-seq数据的质量太差。

c.测序片段的核苷酸组成：对样本的核苷酸组成偏好进行评估。理论上，四种核苷酸在测序片段的每个位置上的比例应该是相似的。但是，用户将会在测序片段的3’端观察到腺嘌呤A/胸腺嘧啶T具有更高的比例，这是因为在测序的cDNA文库中3’端的多聚腺苷（polyA）尾产生的，否则其他位置的腺嘌呤A/胸腺嘧啶T的数目应该和胞嘧啶C/鸟嘌呤G的数目相近。在任何情况下，用户至少在测序片段的3’端可以观察到稳定的模式。钉状的模式（不稳定的模式）出现在测序片段的中间或者末端的位置，表示测序质量低。如果用户同时还观察到测序片段的回帖率低（见“多细胞层面的质量控制报告”），可以从测序片段的3’端修剪一些不稳定的碱基。如果测序数据处理后对回帖率还是没有改善，或许应该考虑是否是Drop-seq数据的质量太差。用户需要注意的是腺嘌呤A/胸腺嘧啶T与胞嘧啶C/鸟嘌呤G的比例会因为物种的不同而产生差异。

d.测序片段的GC含量：测量Drop-seq测序中构建文库的整体质量。

（2）多细胞层面的质量控制报告：对Drop-seq的全部测序片段的质量进行测量。在这个步骤中，不区分细胞条形码是否是“empty”的，把Drop-seq样本作为多细胞RNA-seq样本进行分析。

a.测序片段比对总结：显示所有Drop-seq的测序片段的回帖率和分布。

b.基因区域的覆盖程度：所有基因的测序片段覆盖程度累积图，这个模块目的在于测量Drop-seq数据的整体质量（如图3所示）。

（3）单个细胞层面的质量控制报告：

a.测序片段的重复率分布：给用户展示了每个细胞条形码中的测序片段重复率，从而帮助去除具有低测序片段重复率的细胞条形码（具有低测序片段重复率的细胞条形码通常是由“empty”细胞条形码和环境RNA造成的）。

b.测序片段重复率与覆盖基因数的累积分布：检测用户数据中每个单细胞是否被测序到，以及是否可以从“empty”的细胞条形码中清楚的区分出来。

c.UMI与覆盖的基因数目：在单个细胞层面质量控制，测量Drop-seq实验质量，并从“empty”细胞条形码中区分STAMPs。

d.基因覆盖数目的分布：筛选得到的STAMPs的覆盖基因数目的柱状图，显示筛选得到的STAMPs是否具覆盖足够的基因数目。（如图4所示）

e.内含子比率分布：内含子比率是衡量一个转录组测序样本质量好坏的有效方法，绘制每个STAMPs条形码的内含子比率的柱状图，用来检测来自每个STAMPs的测序片段是否富集在外显子区域。

（4）细胞聚类层面的质量控制报告：

a.Gap statistics：基于t-SNE降维的结果进行k均值聚类来测量样本被区分成不同细胞亚型的能力。

b.聚类图：散点图所展示的是选择得到的STAMPs条形码的t-SNE降维结果的可视化结果图。

c.聚类的Silhouette score：用于解释和验证之前步骤中定义的聚类的类别数的一致性。

在此说明书中，本发明已参照其特定的实施例作了描述，是为便于该技术领域的普通技术人员能理解和使用发明。熟悉本领域技术的人员显然可以容易地对这些实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，不脱离本发明范畴所做出的改进和修改都应该在本发明的保护范围之内。

Claims (9)

1.一种Drop-seq数据质量控制和分析方法，其特征在于：包括以下步骤：

第一步、使用两个配对的测序文件作为输入文件；

第二步、对测序文件进行处理；

第三步、提供以下层面的质量控制的测量；

（1）测序片段层面的质量控制包括测序片段的质量、核苷酸的组成以及测序片段的GC含量；

（2）多细胞层面的质量控制包括基因区域的覆盖程度和序列回帖率；

（3）单个细胞层面的质量控制包括，通过对挑选出的STAMPs计算其测序片段的重复率分布、位于内含子的测序片段的比例和覆盖基因的数目，从而在单细胞层面上对mRNA的捕获效率进行评估；

（4）所述细胞聚类层面的质量控制包括，通过计算gap statistic和silhouettescore，对样本的异质性进行评估；

第四步、给用户提供一份质量控制的报告文档，其中描述了第三步中质量控制的测量值。

2.根据权利要求1所述的Drop-seq数据质量控制和分析方法，其特征在于：还包括以下步骤，

第五步，提供产生以下一项或多项的分析结果，

（1）表达指数；

（2）pair-wise相关表；

（3）主成分分析（PCA）和t-SNE降维的输出结果；

（4）筛选出的STAMPs的聚类簇的分配结果；

（5）t-SNE和聚类结果的可视化输出。

3.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于，所述第一步中的其中一个文件包含转录本的信息，另一个文件包含细胞条形码和UMI的信息。

4.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于，所述第二步中对测序数据的处理方法为把转录本的测序文件比对到参考基因组上，且只留取具有高测序质量的比对上的测序片段。

5.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于：第三步（1）中的测序片段层面的质量控制方法为，将Drop-seq数据看作是多细胞RNA-seq数据进行分析；对测序质量分布、每个测序位点的碱基组成以及每个测序片段的GC含量这几个测度进行计算。

6.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于：所述第三步（2）中的多细胞层面的质量控制部分方法为，列出测序片段的比对结果总结，包括序列的回帖率和全基因组范围的序列分布，并通过绘制基因区域的序列覆盖程度来估计5’端和3’端的测序偏差；使用bedtools工具整合测序片段、细胞条形码以及基因组的注释信息；进行测序片段的合并，把合并后去重的测序片段通过基因的注释信息和细胞条形码信息生成基因表达矩阵，并且利用不同细胞条形码信息把测序片段分配到不同的细胞中；获得的表达矩阵是包含所有细胞的表达矩阵，即同时包括STAMPs和“empty”细胞。

7.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于：所述第三步（3）从混有“empty”细胞的表达矩阵中区分出STAMPs；根据细胞条形码信息把测序片段分组到不同的细胞条形码中；以单个的细胞为研究对象，利用UMI和每个测序片段的基因组上位置的共同信息，计算每个细胞的独特测序片段数目、测序片段的重复率以及覆盖的基因数目；筛选STAMPs条形码：挑选的细胞条形码需要满足其覆盖的基因数目大于用户定义的阈值，评估在单个细胞的层面上mRNA的捕获效率：通过覆盖基因数目的分布和属于内含子测序片段的比例。

8.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于，所述第三步（4）细胞聚类层面质量控制方法为：选择STAMPs条形码的方法产生相对应的表达矩阵，这个表达矩阵以STAMPs条形码为列基因为行，并且表达矩阵中的每个值是由原表达值进行log转换的值；计算表达矩阵中每个基因在所有单细胞中表达的均值和离差值，并且根据基因的平均表达值把所有的基因分配到20个bin中；对每个bin中的所有基因的离差进行Z-标准化，定义z-score的阈值为1.7来确定高可变基因，这些高可变基因将作主成分分析的输入，以确保准确的识别数据的主要结构；在获得的结果中挑选最主要的主成分；使用t-SNE降维方法进一步对这些主成分进行降维到2维层面；对t-SNE的输出结果进行K均值聚类，并且在聚类中运用gap统计确定类的数目，并将STAMPs条形码分配到不同的类中；利用gapstatistics、silhouette score和t-SNE的可视化输出结果来评估细胞聚类的好坏以及样本的异质性；使用gap统计确定K-means聚类的K。

9.根据权利要求1所述Drop-seq数据质量控制和分析方法，其特征在于：所述第一步中的测序文件为FASTQ或SAM格式。