对免疫组库高通量测序样本间序列污染进行过滤的方法
阅读说明：本技术 对免疫组库高通量测序样本间序列污染进行过滤的方法 (Method for filtering sequence contamination among immune repertoire high-throughput sequencing samples ) 是由 张伟 罗礼华 刘晓 于 2019-03-28 设计创作，主要内容包括：本发明公开了一种对免疫组库高通量测序样本间序列污染进行过滤的方法。本发明方法包括对免疫组库高通量测序所得有效数据进行泳道内样本间低频过滤、泳道间样本间低频过滤和核苷酸序列多样性过滤的步骤。在免疫组库(TCR&BCR)建库和测序方法中测序污染一直是个无可避免的问题,影响着数据的可信性。但是,在免疫组库的计算机分析研究领域,系统且通用性的污染过滤方法一直缺失。本发明填补了这块空白,针对测序过程中可能的污染源进行准确的过滤,保证了后续数据分析的准确性。(The invention discloses a method for filtering sequence pollution among immune repertoire high-throughput sequencing samples. The method comprises the step of carrying out intra-lane sample low-frequency filtration, inter-lane sample low-frequency filtration and nucleotide sequence diversity filtration on effective data obtained by high-throughput sequencing of the immune repertoire. Sequencing contamination has always been an inevitable problem in the construction of immune repertoires (TCR & BCR) and sequencing methods, affecting the credibility of the data. However, in the field of computer analysis research of immune repertoires, systematic and versatile contamination filtration methods have been lacking. The method fills the blank, accurately filters possible pollution sources in the sequencing process, and ensures the accuracy of subsequent data analysis.)
技术领域
本发明涉及生物信息学领域，特别涉及一种对免疫组库高通量测序样本间序列污染进行过滤的方法。
背景技术
免疫组库(immune repertoire，IR)，指特定时间内某机体中功能多样性T淋巴细胞和B淋巴细胞的总和。TCR(T cell receptor)是位于T细胞表面的受体。BCR(B cellreceptor)是位于B细胞表面的免疫球蛋白，用于识别抗原、接受抗原的刺激从而启动免疫应答。TCR和BCR都是由两条链组成(重链和轻链，或者α和β链)，分别含450～550个核苷酸或211～217个氨基酸残基。每条链上，有些区域多态性很高，称为可变区(variable region)，其中多样性最高，并能与抗原表位空间互补，叫互补决定区(complementaritydetermining region,CDR)。TCR的α链和BCR重链由V、D、J基因簇编码，TCR的β链和BCR的轻链由V、J基因簇编码。数量众多的V(D)J基因串联排列在同一条染色体，相互之间被内含子分开。在T、B淋巴细胞的发育过程中，V(D)J基因会发生基因重排。除了V(D)J基因随机组合产生大量的多样性，在VD或DJ连接过程中随机***或删除一些核苷酸进一步丰富了受体的多样性。理论上，对于人类来说，将会有超过1018个独特的TCR和2×1012个独特的BCR，从而构成容量巨大的抗原识别受体库，即免疫组库。
为了捕获如此高多样性的免疫组库，一般是利用特异性引物，对TCR和BCR的基因区进行特异性扩增，再进行高通量测序。目前常用的有三种实验捕获方法，即多重PCR(Polymerase Chain Reaction，多聚酶链式反应)，5’RACE(Rapid amplification of cDNAends)技术，随机合成标签序列的UID(unique molecular identifier)技术。近十年内，利用免疫组库技术，已经有很多方面的科学研究和应用，包括于白血病的致病克隆检测以及治疗后的免疫恢复监控、肿瘤的免疫微环境和免疫治疗、疫苗前后的免疫反应及不同疫苗效果的评价、单克隆抗体的快速筛选和对艾滋病感染HIV病毒的中和抗体鉴定等。特别地，免疫组库技术在白血病的治疗恢复监控上，有着巨大的优势，有更好的灵敏度和更系统的免疫评价。近几年，肿瘤免疫微环境研究的快速发展，免疫组库技术也起到很重要的作用，利用TCR作为T淋巴细胞的识别标记，可以准确对淋巴细胞的进化和分化进行分析；同时，免疫治疗中，TCR也是起到决定性的作用。
免疫组库领域，之前的研究都只是用少量的样本，由于少量样本在建库和测序过程中，污染率会很低，并且就算有污染的序列，对项目的研究分析并不会造成影响。当然，少量样本也因为考虑到成本以及对该领域的研究认识局限。最近一年，才开始进行大规模的免疫组库测序，这种污染问题最近才被发现，所以，还没有用于过滤样本间污染序列的分析方法和思路。
免疫组库对TCR和BCR进行实验捕获和扩增时，通常会在一个96孔板进行PCR扩增，一个孔放一个样本。对于高通量测序来说，由于免疫组库样本需要的测序量并不大，所以，一般是有多个样本混合(pooling)在一个测序泳道(lane)，比如48个样本或者96个样本在一个lane上机测序，对于这些样本，一般会加一段标签序列进行区分。这种多样本一起进行扩增和测序的过程，会引入样本之间的相互污染，即一个样本中的序列，最后有少量也出现在另一个样本中。这种序列的污染，对于TCR和BCR的鉴定的多样性，以及分析样本间的共有克隆，造成了极大的干扰。例如，一个TCR本来只在一个样本中出现，但由于样本间的相互污染，最后发现在10个样本中都出现了这个TCR克隆，那么就会错误的认为，这个克隆在这组样本中有富集，或者最后被错误的认为是一个抗原相关性的克隆。
发明内容
针对上述问题，本发明的目的是提供一种利用信息分析的方法，去除免疫组库高通量测序样本间的污染序列。
第一方面，本发明要求保护一种对免疫组库高通量测序样本间序列污染进行过滤的方法。
本发明所提供的对免疫组库高通量测序样本间序列污染进行过滤的方法，可包括如下步骤：
(A)对免疫组库高通量测序所得有效数据进行泳道内样本间低频过滤，输出过滤后数据；
(B)对步骤(A)输出的过滤后数据进行泳道间样本间低频过滤，输出过滤后数据；
(C)对步骤(B)输出的过滤后数据进行核苷酸序列多样性过滤，输出最终有效数据。
步骤(A)中，所述免疫组库高通量测序所得有效数据可按照包括如下步骤的方法获得：对免疫组库高通量测序所得下机数据进行基本处理和常规信息分析(可以用常规的免疫组库分析软件，如IMonitor，使用参数为：-a-b-A1-A2-o-n-t-Rs)，从而获得所述有效数据。
简单来说，包括基础数据处理，V(D)J分配，序列结构分析和数据统计及可视化。首先，过滤低质量的读长(reads)，并合并配对的reads(使序列更长)；其次，将配对后的序列与已有的V/D/J参考序列对比，根据比对得分分配对应的V/D/J基因；第三步，纠正PCR和测序错误，确定序列结构和翻译，过滤掉没有比对到V、J序列以及没有CDR3区域的序列；最后，对有效的序列进行数据统计和图形展示。
在多样本的测序中，频率较高的克隆(“克隆”一词为免疫学专业词语。定义：一种TCR(T细胞受体)或者BCR(B细胞受体)核苷酸序列。下文出现“克隆”含义与此相同)容易污染到同一泳道的其他样本中，因此高频率克隆会在其他多个样本中以相对较低的频率出现。因此，可在步骤(A)中按照包括如下步骤的方法对所述有效数据进行泳道内样本间低频过滤：
(A1)合并同泳道内所有样本的克隆，并统计计算每个克隆在每个样本中的频率；
(A2)如果在同一泳道内的某两个样本间，均存在某一克隆A(核苷酸序列一致)，且所述克隆A在所述两个样本中的频率比大于某个阈值α，则过滤掉所述两个样本中低频样本中的所述克隆A。例如，同一泳道内的样本1和样本2均存在克隆A(核苷酸序列一致)，且所述克隆A在所述样本1中的频率与所述克隆A在所述样本2中的频率的比值大于所述阈值α(且所述克隆A在所述样本1中的频率大于所述克隆A在所述样本2中的频率)，则过滤掉所述样本2中的所述克隆A。
进一步地，所述阈值α可根据同表型样本集(同表型如这些样本是同一种疾病或者同是健康人)在泳道内和泳道间的测序结果，统计克隆频率分布，选择对同一克隆在某一泳道内的频率和在其他泳道内的频率有区分度的值，即为阈值α。
更进一步地，在本发明中所述阈值α具体为2000:1。
对于一批同表型的样本，任意一个克隆，从随机上来说，出现在每个泳道的样本数比例是差不多，如果某个泳道发现更多样本出现某个克隆，则很可能这个泳道有污染情况。因此，可在步骤(B)中按照包括如下步骤的方法对步骤(A)输出的过滤后数据进行泳道间样本间低频过滤：
(B1)合并所有样本的克隆，并统计计算在每个泳道内每个克隆出现的样本数比例；所述样本数比例为某一克隆B在某一泳道a内出现的样本数和所述泳道a中总样本数之比；
(B2)如果某一克隆C在某一泳道b的所述样本数比例高于阈值β，则按照步骤(B3)对所述泳道b中的所述克隆C进行过滤；
(B3)按照所述克隆C在所述泳道b的每个样本中出现的频率由低到高，将所述泳道b中的各样本进行排序，从频率最低(不为0)的样本开始逐一过滤掉样本中的所述克隆C，直至所述克隆C在所述泳道b的所述样本数比例小于或等于所述阈值β。
进一步地，所述阈值β可为所述克隆C在除所述泳道b外的其他泳道的所述样本数比例平均值的5倍。
根据密码子的简并性，同一条氨基酸序列可以由多条不同的核苷酸序列翻译而成。因此，在步骤(C)中可按照包括如下步骤的方法对步骤(B)输出的过滤后数据进行核苷酸序列多样性过滤：
(C1)将所有样本的所有克隆的核苷酸序列均翻译成对应的氨基酸序列，统计在至少N个样本中均出现的公共氨基酸序列；
(C2)如果在均出现某一氨基酸序列M的所有样本中，所述氨基酸序列M均由同一核苷酸序列m翻译而成，则认为所述核苷酸序列m为污染序列，滤去所有样本中的所述核苷酸序列m。
进一步地，所述N个样本可为8-12个样本，如10个样本。
在步骤(C)之后还可包括如下步骤(D)：
(D)针对每一个样本，分别统计所述步骤(A)、所述步骤(B)和所述步骤(C)中每一步被过滤掉的核苷酸序列占总核苷酸序列的百分比；如果某一步骤中某一样本X被过滤掉的核苷酸序列占总核苷酸序列的百分比高于阈值γ，则认为所述样本X的数据被污染太多，不能继续使用，过滤掉所述样本X的所有数据。
进一步地，所述阈值γ可为20％。
第二方面，本发明要求保护一种用于对免疫组库高通量测序样本间序列污染进行过滤的系统。
本发明所提供的对免疫组库高通量测序样本间序列污染进行过滤的系统，可为系统I或系统II；
所述系统I包括装置A、装置B和装置C；
所述系统II包括装置A、装置B、装置C和装置D；
所述装置A能够实现前述第一方面中所述的步骤(A)；所述装置B能够实现前述第一方面中所述的步骤(B)；所述装置C能够实现前述第一方面中所述的步骤(C)；所述装置D能够实现前述第一方面中所述的步骤(D)。
根据需要，所述系统还可包括高通量测序仪和/或能够对免疫组库高通量测序所得下机数据进行基本处理和常规信息分析(可以用常规的免疫组库分析软件，如IMonitor)以获得所述步骤(A)中所述有效数据的仪器。
第三方面，本发明要求保护第二方面中所述系统在对免疫组库高通量测序样本间序列污染进行过滤中的应用。
在上述三方面中，所述样本可为血液、组织中提取的DNA或者RNA。所述高通量测序为多样本混合(pooling)高通量测序，测序平台不限，包括illumina、BGIseq等，测序类型不限，可为单端测序，也可为双端测序。在本发明的一个实施例中具体为单端侧200bp。
将所述样本通过免疫组库实验捕获技术(如多重PCR，5’RACE技术，UID技术等)对T细胞的TCR或者B细胞的BCR进行捕获扩增。这种捕获可以是扩增出全部TCR和BCR的序列，或者仅仅捕获多样性最高的可变区(CDR3区)。然后对扩增产物进行多样本混合(pooling)高通量测序。最后，得到免疫组库高通量测序的下机数据。
在免疫组库(TCR&BCR)建库和测序方法中测序污染一直是个无可避免的问题，影响着数据的可信性。但是，在免疫组库的计算机分析研究领域，系统且通用性的污染过滤方法一直缺失。本发明填补了这块空白，针对测序过程中可能的污染源进行准确的过滤，保证了后续数据分析的准确性。
附图说明
图1为免疫组库样本间序列污染过滤流程图。
图2为克隆在群体中过滤前后分布对比。左图为过滤前，右图为过滤后，横坐标为样本，纵坐标为克隆，黑色表示克隆在样本中存在，白色表示克隆在样本中不存在。
具体实施方式
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。
下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、过滤免疫组库高通量测序样本间的污染序列
本实施例中所采用的免疫组库样本间序列污染过滤流程图如图1所示。
一、实验扩增TCR
从224个正常人的外周血，提取DNA，在96孔板上进行多重PCR捕获并扩增T细胞受体(TCR)，每个孔为一个样本。对这些样本进行建库，然后将每96个样本混合(pooling)到一条泳道上，用BGI-SEQ500测序仪进行高通量测序，测序类型为单端测200bp。
二、免疫组库数据的基本信息分析
使用IMonitor分析下机的测序数据，使用参数为：-a-b-A1-A2-o-n-t-Rs)，经过基础数据处理，V(D)J分配，序列结构分析得到比对正确的CDR3序列数据，对各样本的这些序列进行结构数据统计。表1为原始测序数据经过IMonitor软件分析后的数据统计。以样本S1为例，原始的测序下机数据序列有6264323条，滤去测序接头污染的序列占原始序列比为0.78％，滤去测序质量值低的序列后，98.15％的序列能被保留下来，能比对到V基因片段、D基因片段、J基因片段、同时VJ基因的片段分别为94.54％、59.95％、94.45％、92.19％，测序错误校正后保留99.33％的序列，最后有效数据占原始数据的69.99％，最终S1样本经过IMonitor基本信息分析后得到的有效数据为4384249条。
表1部分样本的基本数据分析结果
样本S1
样本S2
样本S3
样本S4
样本S5
样本S6
原始序列(条)
6264323
22861463
5803853
9737697
8369471
6284164
测序接头污染过滤(％)
0.78
0.86
0.83
0.81
0.76
0.84
测序质量值高的序列比(％)
98.15
98.57
98.38
98.62
98.63
98.2
比对的V基因片段(％)
94.54
98.1
96.29
96.55
95.83
96.92
比对的D基因片段(％)
59.93
13.61
59.59
53.07
64.01
70.15
比对的J基因片段(％)
94.45
98.88
95.52
95.73
95.2
96.48
同时比对V和J基因片段(％)
92.19
97.36
93.72
93.96
93.19
95.06
序列的测序错误纠错率(％)
99.33
99.52
99.27
99.22
99.23
99.27
有效序列比例(％)
69.99
82.23
75.83
74.45
69.05
77.8
有效序列(条)
4384249
18798416
4400856
7249256
5779127
4888830
三、泳道(lane)内样本间的低频过滤
在上述的基本信息分析后，每个克隆能被定义出来，以及计算出了每个克隆在样本中的频率。首先合并同一个泳道内所有样本的克隆数据，利用克隆在不同样本中的频率比较(频率比>2000:1)，过滤掉频率比较低的克隆(即如果在同一泳道内的某两个样本间，均存在某一克隆，且该克隆在所述两个样本中的频率比大于2000:1，则过滤掉所述两个样本中低频样本中的该克隆)。表2中的6个样本中显示，过滤后，绝大部分序列都被保留下来了，只有0.36％～4.65％的序列被认为是污染的序列。
表2部分样本在泳道内样本间的低频过滤的结果
四、泳道间样本间低频过滤
合并所有样本的克隆以及记录克隆在每个样本中的频率，计算每个克隆在每个泳道的样本数比列(样本数比例为某一克隆在某一泳道内出现的样本数和该泳道中总样本数之比)，再比较不同泳道的样本数比例。如果某一克隆在某个泳道的样本数比例大于该克隆在其他泳道的样本数比例的平均值的5倍，则对该泳道中的该克隆进行过滤。从过滤结果看，只有很少数序列被认为是污染的序列，根据表3中6个样本的统计，只有0.01％～0.08％的序列被过滤掉。
表3部分样本在泳道间样本间低频过滤的结果
五、核苷酸序列多样性过滤
将所有样本的所有克隆的核苷酸序列均翻译成对应的氨基酸序列，统计在至少10个样本中均出现的公共氨基酸序列。如果在均出现某一氨基酸序列的所有样本中，该氨基酸序列所对应的都只有一条核苷酸序列，并且在所有样本中一致，则认为该核苷酸序列是污染的序列，过滤掉样本中出现该克隆的序列。从过滤的结果看，表4的6个样本，被过滤掉的序列百分比为0.08％～0.88％。
表4部分样本的核苷酸序列多样性过滤结果
六、过滤后对结果的影响
对于大样本量的免疫组库数据，在多个样本共有的克隆是非常重要的，很可能是特定表型或者疾病相关的克隆。如图2所示，过滤前，发现有一堆克隆在一些样本中有聚集，但经过过滤后，这些克隆聚集现象已经不存在了，说明之前聚集的克隆是样本间的污染所致。因为从以往研究结果来看，由于免疫组库多样性很大，克隆在不同样本中共享的比例非常少，很少有克隆聚集的现象(参考文献：Glanville,J.,Kuo,T.C.,von Budingen,H.C.,Guey,L.,Berka,J.,Sundar,P.D.,Huerta,G.,Mehta,G.R.,Oksenberg,J.R.,Hauser,S.L.et al.(2011)Naive antibody gene-segment frequencies are heritable andunaltered by chronic lymphocyte ablation.Proceedings of the National Academyof Sciences of the United States of America,108,20066-20071.)。