本发明公开了一种基于分子标记的免疫组库生物信息分析流程：选取序列相同的单链免疫球蛋白基因，在选取的单链免疫球蛋白基因的3’端或5’端添加引物和结构相同的UMI序列，添加完成后对免疫基因进行PCR扩增得到测试序列，对测试序列进行过滤，过滤之后通过建立有向无环树对免疫基因所带的UMI序列进行校正，之后对经过校正的UMI序列所标记的免疫基因进行校正，组装校正之后的免疫基因，最后对组装的免疫基因进行统计和报告。本发明通过对UMI自身的校正和免疫基因的校正，有效去除扩增错误和测序错误对免疫基因测序的影响，提高了测序数据的准确度。
1.基于分子标记的免疫组库生物信息分析方法，其特征在于，包括以下步骤：（1）构建双末端测序所需的测试序列：在单链免疫基因上引入扩增引物和UMI标记序列，进行PCR扩增，得到测试序列；
（2）剔除不完整的测试序列：根据是否带有UMI和引物分离步骤（1）所得的测试序列，保留含有UMI和引物的测试序列，去除不带引物或UMI的测试序列；
（3）测序质量控制：根据测序质量值对步骤（2）中保留下来的测试序列进行过滤；
（4）UMI自身校正：根据不同UMI序列之间的汉明距离和每个UMI所标记的免疫基因序列的种类数，将步骤（3）所得的测试序列分成不同的团簇，团簇的数目就是纠正后UMI的种类数；
（5）每个团簇中免疫基因的数目统计：纠正前各个UMI所标记的免疫基因序列的种类数之和就是该UMI所在团簇的免疫基因序列的种类数；
（6）每个团簇中免疫基因序列的校正：对同一个团簇中的免疫基因序列采用多序列对比软件muscle进行相互之间的序列对比，如果某个位点的一致性碱基的比例大于0.6，则该位点为该一致性碱基，反之用N代替，得出每个团簇中免疫基因的序列；
（7）免疫基因序列组装：将各个团簇上的免疫基因进行组装；
（8）统计免疫基因的真实表达量：对步骤（7）中组装好的免疫基因序列进行过滤，除去没有组装序列，并分析不同团簇中免疫基因的相似性，如果一致，则合并，并记录相关的数目信息，即为该基因的真实表达量；
（9）统计与报告：采用igblast进行数据库比对注释，分析注释结果分别统计有功能的和没有功能的基因类型。
2.根据权利要求1所述的分析方法，其特征在于，所述步骤（1）中的UMI是结构为UNNNNUNNNNUNNNNU的四个U碱基框架，两两U碱基之间包含四个随机的碱基，该UMI理论种类有412种。
3.根据权利要求1所述的分析方法，其特征在于，所述步骤（1）中的免疫基因是免疫球蛋白M基因、免疫球蛋白G基因、免疫球蛋白A基因、免疫球蛋白D基因和免疫球蛋白E基因中的一种或多种。
4.根据权利要求1所述的分析方法，其特征在于，所述步骤（1）中在且仅在单链免疫基因的一端添加引物和UMI序列，且单链免疫基因3’端和5’端添加相同结构的UMI序列，一条
3’端带有引物和UMI的免疫基因序列和一条5’端带引物和UMI的相同序列的免疫基因序列称为一对测试序列。
5.根据权利要求1所述的分析方法，其特征在于，所述步骤（3）对整个测试序列进行测序，通过过滤，保留测序质量值在20以上的测试序列。
6.根据权利要求1所述的分析方法，其特征在于，所述步骤（4）中UMI自身校正过程包括以下步骤：①建树：将不同的UMI作为不同的节点，连接汉明距离为1的UMI节点，建成多棵有向无环树；
②赋值：对步骤①中建成的有向无环树的节点进行赋值，所赋数值为该UMI所标记的免疫基因序列的种类数；
③砍树：当节点A所赋数值大于与之相连的节点B所赋数值×2+1时，砍除节点A和节点B之间的边；反之，则保留节点A和节点B之间的边；
④形成团簇：对步骤①中所建树上的每个节点都进行步骤③所述的操作，最终将步骤①中所建树分割成多棵新的树，每棵新树就是一个团簇。
7.根据权利要求1所述的分析方法，其特征在于，所述步骤（7）中的组装，对于全长测序，则根据末端重叠区域进行拼接。
8.根据权利要求1所述的分析方法，其特征在于，所述步骤（7）中的组装，对于非全长序列，采用比对imgt数据库中参考序列的方法进行拼接。
9.根据权利要求1所述的分析方法，其特征在于，所述步骤（7）组装的过程中，基因序列存在缺失现象时，根据imgt数据库中的参考序列对缺失部分进行填充。
10.根据权利要求1所述的分析方法，其特征在于，所述步骤（9）中有功能的基因是指免疫基因中的核酸长度是3的倍数且不含终止密码子的CDR3区基因。
基于分子标记的免疫组库生物信息分析方法
技术领域
[0001] 本发明属于分子生物信息分析处理系统领域，具体地，涉及一种基于分子标记的免疫组库生物信息分析方法。
背景技术
[0002] 高通量测序技术(High-throughput sequencing,HTS)，又称为深度测序技术，是对传统Sanger测序(称为一代测序技术)革命性的改变，能够一次对几十万到几百万条核酸分子进行序列测定，使得对一个物种的转录组和基因组进行细致全貌的分析成为可能。高通量测序技术的发展极大的促进了精准医学的发展，落地了较多高通量测序临床应用，如无创产前基因检测(NIPT)等。
[0003] 免疫组库是指某个个体在任何特定时间点其循环系统中所有功能多样性B淋巴细胞和T淋巴细胞的总和。T淋巴细胞和B淋巴细胞分别通过其表面的细胞受体(TCR或BCR)来识别和结合抗原，进而发挥功能清除病原体或肿瘤细胞等。一个T或B淋巴细胞只表达一种TCR或BCR，每条TCR或BCR都由可变区和恒定区组成，不同克隆T、B细胞的恒定区可相同，但可变区不同，人体T、B淋巴细胞总数约1012，因而具有复杂的识别抗原受体的多样性。
[0004] 免疫组库测序(Immune Repertoire sequencing(IR-SEQ))是以T/B淋巴细胞为研究目标，用多重PCR技术扩增决定B淋巴细胞受体(BCR)或T淋巴细胞受体(TCR)多样性的互补决定区(CDR3区)，再结合高通量测序技术，全面评估免疫系统的多样性，深入挖掘免疫组库与疾病的关系。
[0005] 免疫组库测序作为一种新的高通量测序技术，近年来一直处于科研前沿，特别是随着免疫治疗的兴起与临床落地，极大的推进了免疫组库测序技术的发展。免疫治疗目前拥有较大的市场前景，随着相关产品的批准与落地，也极大的刺激了免疫治疗的研发，免疫组库测序作为免疫治疗研发和预后监控的关键一环，其市场前景也十分巨大，提升免疫组库测序数据的准确性能够极大的促进免疫治疗的研发与临床预后监控效果。免疫组库测序的应用场景不只是免疫治疗，在其他方面也有较多的运用，如抗体研发等等，其应用市场规模巨大，应用场景多样，对其进行研究拥有较大的意义。
[0006] 但目前免疫组库测序有一些难点，如PCR和测序造成的错误无法较好的纠正，会极大的影响后续分析免疫组库的多样性，而免疫组库的多样性是一些临床场景引用的基础，相关的分析流程和算法也无法满足特定场合的免疫组库临床产品的需求。因此，对免疫组库测序技术的研发，解决相关技术难题，对于免疫治疗和癌症预后监控具有重大社会及经济意义。
[0007] 《MiXCR:software for comprehensive adaptive immunity profiling》一文介绍了Mixcr软件在免疫测序中的应用和相对于原有软件的特点。虽然Mixcr软件目前在非分子标记免疫组库测序中用的比较普遍，而这种测序方法会导致分析结果中测序错误和PCR错误较多，没法进行较好的纠正，导致分析结果出现一定的偏差。此外，基于MIGEC软件的免疫组库测序分析也是一种较为常见的方法，但是，其针对的实验类型太窄，泛用性不好，测序结果也存在一定的偏差。
[0008] 中国专利公布CN107122626A公开了一种二代测序DNA突变检测的生物信息学分析的方法和系统，包括生物信息分析模块，用于提供生物分析流程基本组成单元，完成生物信息分析基本功能；中间数据转换模块用于对生物信息分析模块产生的数据进行格式转换，提供符合要求的生物分析数据源和结果；运行环境配置模块用于配置不同生物分析流程运行时所有输入文件、输出文件、配置文件、数据文件、临时文件、日志记录、脚本及应用程序的相对路径或绝对路径以及运行相关环境变量。该发明只针对DNA突变检测这一个具体生物信息分析方法进行设计，并不能很好地用于免疫组库生物信息分析。
[0009] 目前急需一种用于免疫组库生物信息分析且能够纠正过程中出现的错误的分析方法。
发明内容
[0010] 本发明提供了一种基于分子标记的免疫组库生物信息分析方法，该方法采用独特的UMI(unique molecular identifiers)纠错算法，提高了基于分子标记的免疫组库测序生物信息分析的准确性，而且适用范围广。
[0011] 本发明公开了基于分子标记的免疫组库生物信息分析方法，包括以下步骤：
[0012] (1)构建双末端测序所需的测试序列：在单链免疫基因上引入扩增引物和UMI标记序列，进行PCR扩增，得到测试序列；
[0013] (2)剔除不完整的测试序列：根据是否带有UMI和引物分离步骤(1)所得的测试序列，保留含有UMI和引物的测试序列，去除不带引物或UMI的测试序列；
[0014] (3)测序质量控制：根据测序质量值对步骤(2)中保留下来的测试序列进行过滤；
[0015] (4)UMI自身校正：根据不同UMI序列之间的汉明距离和每个UMI所标记的免疫基因序列的种类数，将步骤(3)所得的测试序列分成不同的团簇，团簇的数目就是纠正后UMI的种类数；
[0016] (5)每个团簇中免疫基因的数目统计：纠正前各个UMI所标记的免疫基因序列的种类数之和就是该UMI所在团簇的免疫基因序列的种类数；
[0017] (6)每个团簇中免疫基因序列的校正：对同一个团簇中的免疫基因序列采用多序列对比软件muscle进行相互之间的序列对比，如果某个位置的一致性碱基的比例大于0.6，则该位点为该一致性碱基，反之用N代替，得出每个团簇中免疫基因的序列；
[0018] (7)免疫基因序列组装：将各个团簇上的免疫基因进行组装；
[0019] (8)统计免疫基因的真实表达量：对步骤(7)中组装好的免疫基因序列进行过滤，除去没有组装序列，并分析不同团簇中免疫基因的相似性，如果一致，则合并，并记录相关的数目信息，即为该基因的真实表达量；
[0020] (9)统计与报告：采用igblast进行数据库比对注释，分析注释结果分别统计有功能的和没有功能的基因类型。
[0021] 所述步骤(1)中的UMI是结构为UNNNNUNNNNUNNNNU的四个U碱基框架，两两U碱基之12
间包含四个随机的碱基，该UMI理论种类有4 种。
[0022] 所述步骤(1)中的免疫基因是免疫球蛋白M基因、免疫球蛋白G基因、免疫球蛋白A基因、免疫球蛋白D基因、免疫球蛋白E基因中的一种或多种。
[0023] 所述步骤(1)中在且仅在单链免疫基因的一端添加引物和UMI序列，且单链免疫基因3’端和5’端添加相同的UMI序列，一条3’端带有引物和UMI的免疫基因序列和一条5’端带引物和UMI的相同序列的免疫基因序列称为一对测试序列。
[0024] 所述步骤(3)对整个测试序列进行测序，通过过滤，保留测序质量值在20以上的测试序列。
[0025] 所述步骤(4)中UMI自身校正过程包括以下步骤：
[0026] ①建树：将不同的UMI作为不同的节点，连接汉明距离为1的UMI节点，建成多棵有向无环树；
[0027] ②赋值：对步骤①中建成的有向无环树的节点进行赋值，所赋数值为该UMI所标记的免疫基因序列的种类数；
[0028] ③砍树：当节点A(任意节点)所赋数值大于与之相连的节点B(任意节点)所赋数值×2+1时，砍除节点A和节点B之间的边；反之，则保留节点A和节点B之间的边；
[0029] ④形成团簇：对步骤①中所建树上的每个节点都进行步骤③所述的操作，最终将步骤①中所建树分割成多棵新的树，每棵新树就是一个团簇。
[0030] 所述步骤(7)中的组装，对于全长测序，则根据末端重叠区域进行拼接；对于非全长序列，采用比对imgt数据库中参考序列的方法进行拼接。
[0031] 所述步骤(7)组装的过程中，基因序列存在缺失现象时，根据imgt数据库中的参考序列对缺失部分进行填充。
[0032] 所述步骤(9)中有功能的基因是指免疫基因中的核酸长度是3的倍数且不含终止密码子的CDR3区基因。
[0033] 与现有技术相比，本发明的有益效果为：
[0034] (1)本发明提供了一种新的免疫组库生物信息分析方法，促进了分子标记的免疫组库测序技术的发展，对于免疫治疗和癌症预后监控具有重大社会及经济意义。
[0035] (2)本发明提供的方法适用范围广，而且可以同时对多种类型免疫球蛋白基因进行建库。
[0036] (3)本发明基于分子标记对免疫组库进行测序，通过对16bpUMI序列的校正和同一UMI序列对应的免疫基因序列之间的相互校正，提高了所测数据的准确度。
附图说明
[0037] 图1为本发明免疫组库生物信息分析的操作流程图。
具体实施方式
[0038] 实施例1UMI自身校正之前的测试序列的构建和筛选
[0039] (1)构建双末端测序所需的测试序列：取结构相同的单链免疫球蛋白M基因并分成数量相等的两份，在其中一份的每条免疫球蛋白M基因的3’端添加引物和结构为UAAAGUCCAGUGCAAU的UMI序列，在另一份的每条免疫球蛋白M基因的5’端添加引物和结构为UAAAGUCCAGUGCAAU的UMI序列，其中一条3’端带有引物和UMI的单链免疫球蛋白M基因和一条5’端带有引物和UMI的单链免疫球蛋白M基因称为一对测试序列；取结构相同的单链免疫球蛋白A基因并分成数量相等的两份，在其中一份的每条免疫球蛋白A基因的3’端添加引物和结构为UGGCAUAAGCUAGCAU的UMI序列，在另一份的每条免疫球蛋白A基因的5’端添加引物和结构为UGGCAUAAGCUAGCAU的UMI序列，其中一条3’端带有引物和UMI的单链免疫球蛋白A基因和一条5’端带有引物和UMI的单链免疫球蛋白A基因称为一对测试序列:将上述经过标记的免疫球蛋白基因混合后进行PCR扩增；
[0040] (2)剔除不完整的测试序列：根据是否带有UMI和引物分离步骤(1)所得的测试序列，保留含有UMI和引物的测试序列，去除不带引物或UMI的测试序列；
[0041] (3)测序质量控制：根据测序质量值对步骤(2)中保留下来的测试序列进行过滤，保留测序质量在20以上的测试序列。
[0042] 表1中实验数据显示：经过PCR扩增得的文库中有部分基因不带有引物或UMI，通过筛选将这部分基因进行进行剔除，有利于提高后续过程的检测效率。
[0043] 表1
[0044]
[0045]
[0046] 实施例2 UMI自身校正和免疫基因数目统计
[0047] (4)采用以下步骤，对步骤(3)所得的测试序列进行UMI自身校正：
[0048] ①建树：将不同的UMI作为不同的节点，连接汉明距离为1的UMI节点，建成多棵有向无环树；
[0049] ②赋值：对步骤①中建成的有向无环树的节点进行赋值，所赋数值为该UMI所标记的免疫基因序列的种类数；
[0050] ③砍树：当节点A(任意节点)所赋数值大于与之相连的节点B(任意节点)所赋数值×2+1时，砍除节点A和节点B之间的边；反之，则保留节点A和节点B之间的边；
[0051] ④形成团簇：对步骤①中所建树上的每个节点都进行步骤③所述的操作，最终将步骤①中所建树分割成多棵新的树，每棵新树就是一个团簇，团簇的数目就是校正后UMI的数目。
[0052] (5)每个团簇中免疫基因的数目统计：各个节点UMI所标记的免疫基因序列的种类数之和就是该节点所在团簇的免疫基因序列的种类数。
[0053] 实施例3团簇内部免疫基因序列的校正和组装
[0054] (6)每个团簇中免疫基因序列的校正：对同一个团簇中的免疫基因序列采用多序列对比软件muscle进行相互之间的序列对比，如果某个位置的一致性碱基的比例大于0.6，则该位点为该一致性碱基，反之用N代替，得出每个团簇中免疫基因的序列；
[0055] (7)免疫基因序列组装：将各个团簇上的免疫基因进行组装，对于全长测序，则根据末端重叠区域进行拼接，对于非全长序列，采用比对imgt数据库中参考序列的方法进行拼接，基因序列存在缺失现象时，根据imgt数据库中的参考序列对缺失部分进行填充。
[0056] 实施例4统计与报告
[0057] (8)统计免疫基因的真实表达量：对步骤(7)中组装好的免疫基因序列进行过滤，除去没有组装序列，并分析不同团簇中免疫基因的相似性，如果一致，则合并，并记录相关的数目信息，即为该基因的真实表达量；
[0058] (9)统计与报告：采用igblast进行数据库比对注释，分析注释结果分别统计有功能的和没有功能的基因类型。
[0059] 本发明以多种免疫基因为测序对象，按照上述的实验方法进行标记、扩增、筛选和校正，得到了表2中20个基因文库所示的实验数据，通过比较可以发现，校正后UMI的种类数，明显小于校正前UMI的种类数，且UMI的校正率高达70％，有效降低PCR过程中出现的UMI序列错误对免疫基因测序准确度的影响。
[0060] 表2
[0061]
[0062]
