CN109979528B - 一种单细胞免疫组库测序数据的分析方法 - Google Patents

Abstract

本发明公开了一种单细胞免疫组库测序数据的分析方法，包括如下步骤：步骤S1，原始数据进行数据质量统计与基本质控；步骤S2,对步骤S1所得的高质量的reads与参考V(D)J fragments序列进行序列比对及注释，得到Consensus序列；步骤S3，对步骤S2拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型分型；步骤S4，对步骤S3所得样本中的Clonotype进行CDR3、V/J基因、V‑J paired特征分析；步骤S5，将不同样本中的clonotype进行合并，然后进行多样本分析。该方法分析流程完整，内容丰富，结果能够可视化展示，能够全面揭示免疫组库信息，从而使用户能够全面了解评估免疫组库状况，可以更精确地反应生物体内的免疫功能。

Description

一种单细胞免疫组库测序数据的分析方法

技术领域

本发明涉及免疫组库测序技术领域，尤其是一种单细胞免疫组库测序数据的分析方法。

背景技术

免疫组库是指在特定时间段内，机体所有有功能的T细胞和B细胞的总和。人体的适应性免疫系统主要依靠于T淋巴细胞和B淋巴细胞表面的T细胞受体(TCR)和B细胞受体(BCR)上的互补决定区(CDR)与MHC-抗原肽分子进行识别并特异性结合。BCR/TCR存在一块区域叫作互补决定区(CDR),由CDR1、CDR2、CDR3组成，CDR1和CDR2都是由V基因编码，而CDR3是由V、D、J三个基因编码形成，所以CDR3多样性程度最高，是主要参与抗原识别的区域。

免疫组库测序(Immune Repertoire sequencing(IR-SEQ))是以T/B淋巴细胞为研究目标，通过多重PCR或5’RACE技术特异性扩增决定B淋巴细胞受体(BCR)或T淋巴细胞受体(TCR)多样性的互补决定区(CDR3区)，再结合高通量测序技术，全面评估免疫系统的多样性，深入挖掘免疫组库与疾病的关系。免疫组库测序技术应用非常广泛，目前已经在肿瘤预后评价、白血病微小残留病变监测、自身免疫性疾病诊断和疾病生物标志物开发等方面取得了一定的进展。

免疫系统本身就是一个非常复杂的体系，传统的方法下科学家们大都是关心免疫系统或某类免疫细胞的整体功能。但生物体内广泛存在的细胞异质性，相关因子表达的时效性、动态变化过程等对研究结果影响巨大。而单细胞测序是可以在单个细胞水平上非常准确地鉴定两条TCR/BCR链配对情况，从而使得分析达到更高的复杂水平，同时也可以更加精确地反应生物体内的功能。

近年来兴起的高通量测序技术虽然能够大规模地测序多样化的免疫基因，但是该技术还是存在测序错误、测序长度不足、测序深度不足等问题，而且高通量测序产生大量的序列信息，也需要建立一套准确完整的数据分析流程。

发明内容

基于上述问题，本发明的目的在于克服上述现有技术的不足之处而提供一种单细胞免疫组库测序数据的分析方法。

为实现上述目的，本发明采取的技术方案为：

一种单细胞免疫组库测序数据的分析方法，包括如下步骤：

步骤S1，对测序仪产生的T/B淋巴细胞基因组序列数据进行数据质量统计与基本质控，以去除含有测序接头、含有未知碱基及低质量的Reads,筛选出高质量的UMI和高质量的Barcode，得到高质量的reads；

步骤S2,对步骤S1所得的高质量的reads与参考V(D)J fragments序列进行序列比对及注释，得到Consensus序列(一致序列)；

步骤S3，对步骤S2拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型(克隆亚型)分型；

步骤S4，对步骤S3所得样本中的Clonotype进行CDR3、V/J基因、V-J paired特征分析；

步骤S5，将不同样本中的clonotype进行合并，然后进行多样本分析。优选地，利用10X Genomics官方分析软件Cell Ranger对原始数据进行数据质量统计与基本控制。需要说明的是，测序仪产生的T/B淋巴细胞序列数据即原始数据，是测序仪产生的原始图像数据经base calling转化为序列数据，包含reads的序列和reads的测序数据。

优选地，所述步骤S1中，高质量UMI筛选包括步骤：

根据每个UMI的reads支持数(RPU)，计算Barcode的N50RPU，在包含50％的reads中，选择最多reads支持的Barcodes，计算第99百分位数的N50RPU，过滤掉小于第99百分位数的N50RPU的Barcodes；利用二维高斯混合模型，计算log(N50RPU per barcode)；将较高均值所在的维度的最小RPU作为阈值，若UMI的RPU在阈值之上，则该UMI为高质量UMI；

所述高质量Barcode筛选的标准为：一个Barcode需要有2个以上的高质量UMI支持。需要说明的是，“较高均值”是先确定log(N50RPU per barcode)均值最高的维度，将这个维度的最小RPU作为阈值。

优选地，所述步骤S2中，序列比对与注释，具体包括：

数据比对：Cell Ranger首先将高质量reads与参考V(D)J fragments序列进行比对，将至少有15bp比对到该片段区域的paired reads进行保留；

Contig组装：Cell ranger对于通过mapping及UMI的reads可用于进行组装；

Barcode筛选：Cell ranger对Barcode进行有效性筛选；

组装Contig注释，包括步骤：(1)组装的Contig序列通过Smith-Waterman与germline序列比对，比对到v d j c区，并进行标记；(2)序列找起始密码子，寻找CDR3motif(Cys-FGXG/WGXG)；(3)根据以下条件筛选contig标记为productive，条件：跨越v-j区，v区包含起始密码子，包含cdr3区，v-j区不包含终止密码子；

Contig筛选：保留通过Barcode筛选的细胞，并且在Contig注释有2个以上的UMI支持的productive Contig的Contig；

Consensus序列组装注释:对于Contig筛选后得到的所有Contig进一步拼接，得到样本CDR3支持的Consensus序列(一致序列)。

优选地，所述步骤S3中，Clonotype(克隆亚型)分型是Cell Ranger根据拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型，其有效Barcode为该Clonotype的丰度。

优选地，所述步骤S4中，CDR3特征分析结果用CDR3丰度分布图和CDR3碱基序列长度分布图来展示；

V/J基因特征分析，包括对V/J基因丰度分布情况进行统计分析、CDR3V/J基因序列长度分布情况统计分析及样本间的V/J基因表达水平进行聚类分析；

V-J paired特征分析，包括对V-J paired丰度分布和频率分布情况进行统计分析。

优选地，所述步骤S5中，多样本分析包括样本间Clonotypes比较、OverlappingClonotype聚类分析及Overlapping Clonotype差异分析。

优选地，所述样本间Clonotypes比较，首先是将不同样本间相同的clonotype的名字进行合并，然后用韦恩图来展示样本间共有(overlapping)及特有(nonoverlapping)的Clonotype数目，若样本数量大于5个，用葵花图进行展示。

优选地，所述Overlapping Clonotype聚类分析结果通过绘制样本MDS图和样本聚类图来进行聚类分析。

优选地，所述Overlapping Clonotype差异分析，包括对样本间top 10的overlapping Clonotype进行丰度比较分析以及通过绘制样本间相同V-J pairs频率Ciros图进行差异分析。由此，同一Clonotype在不同样本中的丰度差异，反映出不同状态的免疫信息，为研究免疫系统提供方向。

综上所述，本发明的有益效果为：

(1)本发明提供了一种10X单细胞免疫组库测序数据分析方法，对于免疫组库测序技术的发展起到了很好地促进作用，对于免疫治疗和癌症预后监控有着重大意义；

(2)本发明提供的分析方法，是通过特异性提取富集T/B细胞TCR/BCR，通过比对、拼接、筛选、注释等构建免疫组库，分析得到TCR/BCR多样性序列，从而获得机体的免疫特征；该方法分析流程完整，内容丰富，结果能够可视化展示，能够全面揭示免疫组库信息，从而使用户能够全面了解评估免疫组库状况；

(3)本发明是在单细胞水平特异性富集T/B细胞TCR/BCR，进行测序，测序精度更高，能够获得成对的重链和轻链(B细胞)或α和β链(T细胞)的V(D)J全长序列信息，可以更精确地反应生物体内的免疫功能。

附图说明

图1为10X单细胞免疫组库测序数据分析流程示意图；

图2为本发明实施例中clonotype频率分布甜甜圈图；

图3为本发明实施例中CDR3碱基序列长度分布图；

图4为本发明实施例中V-J paired丰度分布热图；

图5为本发明实施例中样本间相同clonotype丰度比较图。

具体实施方式

单细胞免疫组库测序技术是刚发展起来的新技术，为了解决现有技术存在的不足，本发明提供一种10X单细胞免疫组库测序数据分析方法，该分析方法是通过对下机数据进行比对、拼接、筛选、注释等构建免疫组库，分析得到TCR/BCR多样性序列，从而获得机体的免疫特征。本发明的分析方法能够实现成对的重链和轻链(B细胞)或α和β链(T细胞)的V(D)J全长测序，分析内容丰富，也能更加精确地反应机体的内部功能。

在一些实施例中，本发明公开了一种10X单细胞免疫组库测序数据分析流程，包括如下步骤：步骤S1,对原始数据(即测序仪产生的T/B淋巴细胞基因组序列数据)进行数据质量统计与基本质控；步骤S2,进行序列比对及注释；步骤S3，Clonotype(克隆亚型)分型；步骤S4，特征分析，包括CDR3、V/J基因、V-J paired特征分析；步骤S5，多样本分析，包括样本间Clonotypes比较、Overlapping Clonotype聚类分析及Overlapping Clonotype差异分析。通过上述方法，能够实现成对的重链和轻链(B细胞)或α和β链(T细胞)的V(D)J全长测序，从而获得机体的免疫特征，研究精度更细，也能更加精确地反应机体的内部功能。

在一些实施例中，本发明提出了一种10X单细胞免疫组库测序数据分析方法，包括如下步骤：

步骤S1，对原始数据(即测序仪产生的T/B淋巴细胞基因组序列数据)进行数据质量统计与基本质控：利用10X Genomics官方分析软件Cell Ranger对原始数据进行数据质量统计与基本控制；

步骤S2，进行序列比对及注释：Cell Ranger首先将高质量reads与参考V(D)Jfragments序列进行比对，将至少有15bp比对到该片段区域的paired reads进行保留，然后对于每个细胞进行Contig组装，将组装得到的Contig通过Barcode筛选、Contig注释和Contig筛选，最终得到有效的Contig，对有效的Contig进一步组装，得到Consensus序列(一致序列)，并进行注释；

步骤S3，Clonotype(克隆亚型)分型：Cell Ranger根据拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型，其有效Barcode为该Clonotype的丰度；

步骤S4,特征分析，包括CDR3、V/J基因、V-J paired特征分析；CDR3区为TCR的高突变区，且直接识别抗原呈递细胞呈递的抗原氨基酸片段。CDR3特征分析主要是对CDR3丰度分布情况和碱基序列长度分布情况进行统计分析；V/J基因可以反映出clonotype的特征，V/J基因特征分析，主要是对V/J基因丰度分布情况进行统计分析、CDR3V/J基因序列长度分布情况统计分析及样本间的V/J基因表达水平进行聚类分析等；V-J paired直接反应出CDR3或免疫组库的变化，V-J paired特征分析，主要是对V-J paired丰度分布和频率分布情况进行统计分析；

步骤S5，多样本分析，包括：样本间Clonotypes比较、Overlapping Clonotype聚类分析及Overlapping Clonotype差异分析。

为更好的说明本发明的目的、技术方案和优点，下面将结合附图和具体实施例对本发明作进一步说明。本发明亦可通过其它不同的具体实例加以施行或应用，本说明书中的各项细节亦可基于不同观点与应用，在不背离本发明的精神下进行各种修饰与变更。

实施例1

本发明的单细胞免疫组库测序数据的分析方法的一种实施例，是一种10X单细胞免疫组库测序数据分析流程，图1为10X单细胞免疫组库测序数据分析流程图，具体包括步骤如下：

步骤S1，对原始数据(测序仪产生的原始图像数据经base calling转化为序列数据，包含reads的序列和reads的测序质量)进行数据质量统计与基本质控：利用10XGenomics官方分析软件Cell Ranger对原始数据进行数据质量统计与基本控制，以去除含有测序接头、含有未知碱基及低质量的Reads，筛选出高质量的UMI和高质量的Barcode，得到高质量的reads；

步骤S2,数据比对与注释，具体包括步骤：

S2.1数据比对：Cell Ranger首先将高质量reads与参考V(D)J fragments序列进行比对，将至少有15bp比对到该片段区域的paired reads进行保留；

S2.2Contig组装：Cell ranger对于通过mapping及UMI的reads可以用于进行组装。

S2.3Barcode筛选：Cell ranger对Barcode进行有效性筛选。首先对Barcode与10xGenomics试验中涉及的所有Barcode进行注释，去除人工产物等导致的错误Barcode信息的序列。

高质量Barcode筛选的标准为：对于可标识一个细胞的Barcode需要有2个以上的高质量UMI支持。高质量UMI筛选包括步骤：根据每个UMI的reads支持数(RPU)，计算Barcode的N50RPU。在包含50％的reads中，选择最多reads支持的Barcodes，计算第99百分位数的N50RPU，过滤掉小于第99百分位数的N50RPU的Barcodes。利用二维高斯混合模型(2-component Gaussian mixturemodel，GMM)，计算log(N50RPU per barcode)。将较高均值所在的维度的最小RPU作为阈值，若UMI的RPU在阈值之上，则该UMI为高质量UMI。

S2.4组装Contig注释：(1)组装的Contig序列通过Smith-Waterman与germline序列比对，比对到v d j c区，并进行标记；(2)序列找起始密码子，寻找CDR3motif(Cys-FGXG/WGXG)；(3)根据以下条件筛选contig标记为productive，条件:跨越v-j区，v区包含起始密码子，包含cdr3区，v-j区不包含终止密码子。

S2.5Contig筛选：保留通过Barcode筛选的细胞，并且在Contig注释有2个以上的UMI支持的productive Contig的Contig。

S2.6Consensus序列组装注释:对于Contig筛选后得到的所有Contig进一步拼接，得到样本CDR3支持的Consensus序列(一致序列)。

步骤S3，Clonotype(克隆亚型)分型:根据拼接得到的样本一致序列的中的CDR3氨基酸序列进行Clonotype分型。细胞支持数(Barcode个数)即为该Clonotype的丰度。可以用频数分布图或频率分布图(如图2所示)来展示clonotype丰度分布情况，样本clonotypes稀释曲线来展示样品clonotype的丰富程度，或者利用圈图来展示免疫组库的top10clonotype频率分布。

步骤S4，特征分析：主要包括CDR3特征分析、V/J基因特征分析及V-J paired特征分析。

其中，关于CDR3特征分析，样本中的CDR3丰度分布情况常用CDR3丰度分布图来展示，可以用频数或频率分布图。CDR3碱基序列长度分布情况统计分析结果常通过绘制CDR3碱基序列长度分布分布图来展示(如图3所示)。

V/J基因特征分析，通常包括对V/J基因丰度分布情况和CDR3V/J基因序列长度分布情况进行统计分析。分析结果分别用样本中V基因丰度分布图、J基因丰度分布图和CDR3(V Gene)碱基序列长度分布图、CDR3(J Gene)碱基序列长度分布图来进行结果展示。

V-J paired特征分析，V-J paired直接反应出CDR3或免疫组库的变化。分析结果通常用V-J paired丰度分布热图和V-J paired频率Ciros图来表示，如图4所示。

步骤S5，多样本分析：主要包含样本间Clonotypes比较、Overlapping Clonotype聚类分析及Overlapping Clonotype差异分析3部分。

样本间Clonotypes比较，由于Cell Ranger对各个样本单独处理，每个clonotype在不同样本间可能会被分配到不同的名字，首先需要这些名字进行合并处理。接着用韦恩图来展示样本间共有(overlapping)及特有(nonoverlapping)的Clonotype数目，如果样本数量大于5个，则用葵花图来进行展示。

Overlapping Clonotype聚类分析，使用样本间Overlapping Clonotype丰度进行聚类，反映样本间共有Clonotype的相似性(如图5所示)。结果可用R包中MASS中的isoMDS函数对多个样品进行分群，绘制出样本MDS图进行结果展示，样本间的距离越近，相似度越高。也可以R包stat的hclust函数进行层次聚类分析绘图，分析结果以样本聚类图的形式进行展示，图中相邻分分支间的关系最接近。

Overlapping Clonotype差异分析，通常是对样本间top 10的overlappingClonotype进行丰度比较分析，也可以通过绘制样本间相同V-J pairs频率Ciros图进行差异分析。由此，同一Clonotype在不同样本中的丰度差异，反映出不同状态的免疫信息，为研究免疫系统提供方向。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (8)

1.一种单细胞免疫组库测序数据的分析方法，其特征在于，包括如下步骤：

步骤S1，对测序仪产生的T/B淋巴细胞基因组序列数据进行数据质量统计与基本质控，以去除含有测序接头、含有未知碱基及低质量的Reads, 筛选出高质量的UMI和高质量的Barcode，得到高质量的reads；

步骤S2, 对步骤S1所得的高质量的reads与参考V(D)J fragments 序列进行序列比对及注释，得到Consensus序列；

步骤S3，对步骤S2拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型；

步骤S4，对步骤S3所得样本中的Clonotype进行CDR3、V/J基因、V-J paired特征分析；

步骤S5，将不同样本中的clonotype进行合并，然后进行多样本分析；

所述步骤S2中，序列比对与注释，具体包括：

数据比对：Cell Ranger首先将高质量reads与参考V(D)J fragments 序列进行比对，将至少有15bp比对到该片段区域的paired reads进行保留；

Contig组装：Cell ranger对于通过mapping及UMI的reads可用于进行组装；

Barcode 筛选：Cell ranger对Barcode进行有效性筛选；

组装Contig注释，包括步骤：(1)组装的Contig序列通过Smith-Waterman与germline序列比对，比对到v d j c 区，并进行标记；(2) 序列找起始密码子，寻找CDR3 motif Cys-FGXG/WGXG；(3) 根据以下条件筛选contig标记为productive，条件：跨越v-j区，v区包含起始密码子，包含cdr3区，v-j区不包含终止密码子；

Contig筛选：保留通过Barcode筛选的细胞，并且在Contig注释有2个以上的UMI支持的productive Contig的Contig；

Consensus序列组装注释: 对于Contig筛选后得到的所有Contig进一步拼接，得到样本CDR3支持的Consensus序列。

2.根据权利要求1所述的分析方法，其特征在于，所述步骤S1中，高质量UMI筛选包括步骤：

根据每个UMI的reads支持数，计算Barcode的N50 RPU，在包含50%的reads中，选择最多reads支持的Barcodes，计算第99百分位数的N50 RPU，过滤掉小于第99百分位数的N50 RPU的Barcodes；利用二维高斯混合模型，计算log(N50 RPU per barcode)；将较高均值所在的维度的最小RPU作为阈值，若UMI的RPU在阈值之上，则该UMI为高质量UMI；

所述高质量Barcode筛选的标准为：一个Barcode需要有2个以上的高质量UMI支持。

3.根据权利要求1所述的分析方法，其特征在于，所述步骤S3中，Clonotype分型是CellRanger根据拼接得到的Consensus中的CDR3氨基酸序列进行Clonotype分型，其有效Barcode为该Clonotype的丰度。

4.根据权利要求1所述的分析方法，其特征在于，所述步骤S4中，CDR3特征分析结果用CDR3丰度分布图和CDR3碱基序列长度分布图来展示；

V/J基因特征分析，包括对V/J基因丰度分布情况进行统计分析、CDR3 V/J基因序列长度分布情况统计分析及样本间的V/J基因表达水平进行聚类分析；

V-J paired特征分析，包括对V-J paired 丰度分布和频率分布情况进行统计分析。

5.根据权利要求1所述的分析方法，其特征在于，所述步骤S5中，多样本分析包括样本间Clonotypes比较、 Overlapping Clonotype聚类分析及Overlapping Clonotype差异分析。

6.根据权利要求5所述的分析方法，其特征在于，所述样本间Clonotypes比较，首先是将不同样本间相同的clonotype 的名字进行合并，然后用韦恩图来展示样本间共有及特有的Clonotype数目，若样本数量大于5个，用葵花图进行展示。

7.根据权利要求5所述的分析方法，其特征在于，所述Overlapping Clonotype聚类分析结果通过绘制样本MDS图和样本聚类图来进行聚类分析。

8.根据权利要求5所述的分析方法，其特征在于，所述Overlapping Clonotype差异分析，包括对样本间top 10 的overlapping Clonotype进行丰度比较分析以及通过绘制样本间相同V-J pairs 频率Circos图进行差异分析。

Images