确定混合测序数据中读段的样本源的方法及装置
技术领域
本发明涉及混合数据处理技术,特别是混合测序数据中数据来源的样本的确定方法和装置。
背景技术
Sanger测序是确定基因分型的金标准,飞行时间质谱检测能够实现定点检测基因分型,比如深圳华大基因推出的一款产品针对四个耳聋常见突变基因的20个位点进行质谱检测,这20个位点在我国耳聋人群的致病因素中占据主要作用,还有全外显子组测序,三种方法都具有各自的局限性,比如Sanger和质谱法通量低、成本高,而全外显子组测序则不能有效利用全部测序数据。
先天性耳聋是一类常见疾病,在我国新生儿中的发病率高于1‰,其中60%以上是遗传因素导致的。因此,除了常规的医学诊断方法,通过测定相关基因的基因分型、判断是否发生基因突变,可以辅助医生诊断新生儿是否患有耳聋。
根据国内研究人员针对我国人群中耳聋基因突变进行分子流行病学调查的结果,GJB2、GJB3、SLC26A4和12sRNA的突变最为常见,在人群中的突变比例高达40%,在这四个基因上的突变位点是导致遗传性耳聋发生的常见突变。
发明内容
本发明一方面提供了一种确定混合测序数据中读段的样本源的方法,混合测序数据由多个读段组成,该方法包括:A.利用多个标签分别标记多个核酸样本,使得每个核酸样本至少带有一条或多条标签以区分其它核酸样本,记录所述每个核酸样本与其所带的标签的对应关系;B.混合所述多个核酸样本,对混合核酸样本进行核酸序列测定,获得混合测序数据;C.将所述混合测序数据与参考序列比对,获得比对结果,从所述比对结果中筛选出与所述参考序列不完全匹配而且不匹配部分的长度不小于比A中的标签小1bp的长度的读段。
本发明另一方面提供了一种确定混合测序数据中读段的样本源的装置,混合测序数据由多个读段组成,该装置包括:样本标记单元,用以实现利用多个标签分别标记多个核酸样本,使得每个核酸样本至少带有一条或多条标签以区分其它核酸样本,记录每个核酸样本与其所带的标签的对应关系;混合测序单元,与样本标记单元相连,用以混合获自样本标记单元的标签标记过的多个核酸样本,以及对混合核酸样本进行核酸序列测定,获得混合测序数据;比对筛选单元,与混合测序单元相连,用以实现将混合测序数据与参考序列比对,获得比对结果,以及从比对结果中筛选出与参考序列不完全匹配而且不匹配部分的长度不小于比所述标签小1bp的长度的读段;归类单元,与样本标记单元和比对筛选单元相连,用以实现利用从比对筛选单元中筛选出的读段的信息和所述样本标记单元的对应关系,确定混合测序数据中读段源于的核酸样本。
利用本发明一方面提供的确定混合测序数据中数据的样本源的方法或装置,将多个样本核酸的混合测序后的混合数据正确对应到样本源,使得不浪费测序通量,特别是适合于每个样本数据量需求相对低而测序通量相对高的平台。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解,其中:
图1是本发明的一个具体实施方式中的文库构建示意图;
图2是本发明的一个具体实施方式中的确定混合测序数据中读段的样本源的装置示意图。
具体实施方式
根据本发明的一个实施方式,提供了一种确定混合测序数据中读段的样本源的方法,混合测序数据是由多个读段组成,所说方法包括:
A.利用多个标签分别标记多个核酸样本,使得每个核酸样本至少带有一条或多条标签以区分其它核酸样本,记录所述每个核酸样本与其所带的标签的对应关系;
B.混合所述多个核酸样本,对混合核酸样本进行核酸序列测定,获得混合测序数据;
C.将所述混合测序数据与参考序列比对,获得比对结果,从所述比对结果中筛选出与所述参考序列不完全匹配而且不匹配部分的长度不小于比A中的标签小1bp的长度的读段;
D.依据C中筛选出的读段的信息和A中的对应关系,确定所述混合测序数据中读段源于的核酸样本。
根据本发明的一个具体实施方式,A中标签长度为5~12bp。A中的标签可以选自SEQ ID NO:27~124所示的序列。SEQ ID NO:27~124序列见表1,这组标签,是发明人考虑序列长度、碱基组成、碱基位置比例、与其它标签碱基的关系设计大量序列,多次试验筛选获得的,这组标签的部分或者全部可以置于同一反应体系中而又相互之间不干扰影响,而且不干扰常规体系内的其它反应物或反应,比如不影响文库构建中的各反应体系及反应,测序芯片上的固定序列等。
表1
根据本发明的一个具体实施方式,A中利用标签标记核酸样本是通过标签引物扩增所述核酸样本的至少一部分核酸来实现的,所说的标签引物由位于5‘端的标签连接引物序列构成。所用的标签可以自己设计也可以依据已知公开的,比如购买Illumina或Lifetechnologies公司设计验证的包含可混合一起使用的标签的试剂盒。标签引物中的标签还可选自SEQ ID NO:27~124所示的序列,所说的标签引物中的引物序列可以选自SEQ IDNO:1和2,SEQ ID NO:3和4,SEQ ID NO:5和6,SEQ ID NO:7和8,SEQ ID NO:9和10,SEQ IDNO:11和12,SEQ ID NO:13和14,SEQ ID NO:15和16,SEQ ID NO:17和18,SEQ ID NO:19和20,SEQ ID NO:21和22,SEQ ID NO:23和24以及SEQ ID NO:25和26所示的13对序列中的至少1对。序列SEQ ID NO:1-26如表2,这组构成标签引物的引物,可以扩增各样本中耳聋相关基因区域,并使每个样本的扩增产物中带有相应的一个或多个有顺序关系的标签,这样建立样本与标签的对应关系,依据该对应关系,能够依据本发明的这一实施方式的办法使混合的样本核酸序列数据,通常是数目庞大的混合数据,对应至正确的样本源,分析各样本核酸信息。利用本发明的这一方法,使多个样本的数据能够混合在一起处理而且最终能利用该标记对应关系区分混合数据,归类其至各样本。
表2
引物编号 |
引物序列 |
F1 |
TCTTTTCCAGAGCAAACCGC(SEQ ID NO:1) |
F2 |
ACGTGCATGGCCACTAGGAG(SEQ ID NO:3) |
F3 |
TGCAGCTGATCTTCGTGTCC(SEQ ID NO:5) |
F4 |
ATGGTGAGTACGATGCAGAC(SEQ ID NO:7) |
F5 |
GCCTTTGGTGTGCTAAAGAC(SEQ ID NO:9) |
F6 |
GGGTTCCAGGAAATTACTTTG(SEQ ID NO:11) |
F7 |
AAATGATCGGTTTAGACAC(SEQ ID NO:13) |
F8 |
AGGATCGTTGTCATCCAGTC(SEQ ID NO:15) |
F9 |
TAGGGCCTATTCCTGATTGG(SEQ ID NO:17) |
F10 |
CCAAAGCTCCAAATGTATA(SEQ ID NO:19) |
F11 |
AGAAAAGCTGGAGCAATGCG(SEQ ID NO:21) |
F12 |
ACACACAATAGCTAAGACCC(SEQ ID NO:23) |
F13 |
GAGTGCTTAGTTGAACAGGG(SEQ ID NO:25) |
R1 |
GGGTGTTGCAGACAAAGTCG(SEQ ID NO:2) |
R2 |
TTGTGGCTGCAAAGGAGGTG(SEQ ID NO:4) |
R3 |
ACCACAGGGAGCCTTCGATG(SEQ ID NO:6) |
R4 |
CAAGCTCATCATTGAGTTCC(SEQ ID NO:8) |
R5 |
GGAGAAGTGTTAAACTCCTG(SEQ ID NO:10) |
R6 |
ACAGCTAGAGTCCTGATTGC(SEQ ID NO:12) |
R7 |
TTTCCAGGTTGGCTCCATAT(SEQ ID NO:14) |
R8 |
AAGGCTGTTGTTCCTACCTG(SEQ ID NO:16) |
R9 |
CCAGTCCTATTTTCTATGGC(SEQ ID NO:18) |
R10 |
GTGGATTGGAACTCTGAGC(SEQ ID NO:20) |
R11 |
GATACATCTGTAGAAAGGTTG(SEQ ID NO:22) |
R12 |
GATTACAGAACAGGCTCCTC(SEQ ID NO:24) |
R13 |
AAGCTACACTCTGGTTCGTC(SEQ ID NO:26) |
根据本发明的一个具体实施方式,所说的标签引物中的引物序列还可选自SEQ IDNO:1和2,SEQ ID NO:3和4,SEQ ID NO:5和6,SEQ ID NO:7和8,SEQ ID NO:9和10,SEQ IDNO:11和12,SEQ ID NO:13和14,SEQ ID NO:15和16,SEQ ID NO:17和18,SEQ ID NO:19和20,SEQ ID NO:21和22,SEQ ID NO:23和24以及SEQ ID NO:25和26所示的13对序列中的至少2对。
根据本发明的一个具体实施方式,所说的标签引物中的引物序列选自SEQ ID NO:1和2,SEQ ID NO:3和4,SEQ ID NO:5和6,SEQ ID NO:7和8,SEQ ID NO:9和10,SEQ ID NO:11和12,SEQ ID NO:13和14,SEQ ID NO:15和16,SEQ ID NO:17和18,SEQ ID NO:19和20,SEQ ID NO:21和22,SEQ ID NO:23和24以及SEQ ID NO:25和26所示的13对序列中的至少5对。
根据本发明的一个具体实施方式,所说的标签引物中的引物序列选自SEQ ID NO:1和2,SEQ ID NO:3和4,SEQ ID NO:5和6,SEQ ID NO:7和8,SEQ ID NO:9和10,SEQ ID NO:11和12,SEQ ID NO:13和14,SEQ ID NO:15和16,SEQ ID NO:17和18,SEQ ID NO:19和20,SEQ ID NO:21和22,SEQ ID NO:23和24以及SEQ ID NO:25和26所示的13对序列中的至少10对。
根据本发明的一个具体实施方式,所说的标签引物中的引物序列为SEQ ID NO:1和2,SEQ ID NO:3和4,SEQ ID NO:5和6,SEQ ID NO:7和8,SEQ ID NO:9和10,SEQ ID NO:11和12,SEQ ID NO:13和14,SEQ ID NO:15和16,SEQ ID NO:17和18,SEQ ID NO:19和20,SEQID NO:21和22,SEQ ID NO:23和24以及SEQ ID NO:25和26所示的13对序列。
根据本发明的一个具体实施方式,在获得混合测序数据之后,去除混合测序数据中长度不小于50bp的读段,去除短的读段有利于提高整体数据质量。
根据本发明的一个具体实施方式,C进一步包括依据不匹配部分在读段中的位置对筛选出的读段进行分类,获得第一读段和第二读段,第一读段中的读段的两个末端都与参考序列不匹配,第二读段中的读段的两个末端中的一个与参考序列不匹配。在本发明的一个具体实施方式中,参考序列是由目标序列构成的,将完全比对上特定参考序列的完全不包含标签信息的读段剔除。
根据本发明的一个具体实施方式,比对第一读段中的每个读段中的与参考序列都不匹配的两个末端,去除两个末端不互相匹配的以及两个末端互相匹配的长度小于比A中标签小1bp的长度的读段。比如,标签的长度是7bp,读段两末端互相匹配的长度至少为6bp,该读段才会得以保留用于归类,这是因为我们这组标签能够容纳1bp的错配,即使用该组标签的一部分或全部时,对比判断得的一个标签序列与它自己原本序列有1bp错配,仍旧认为是同一标签,仍旧能将该标签所在读段正确归到该标签标记的样本,这样利于提高数据利用率。
根据本发明的一个具体实施方式,B中核酸序列测序是在半导体芯片测序平台上进行的,B中核酸序列测定包括混合核酸样本的测序文库的构建。测序可以从现有测序平台中选择,可选测序平台包括但不限于CG(Complete Genomics)、Illumina/Solexa、ABI/SOLiD、Roche 454和单分子测序平台,依据所选测序平台进行单端或双端测序文库的制备。
根据本发明的另一实施方式,提供一种确定混合测序数据中读段的样本源的装置,混合测序数据是由多个读段组成的,所说装置1000如图2所示,包括:样本标记单元100,用以实现利用多个标签分别标记多个核酸样本,使得每个核酸样本至少带有一条或多条标签以区分其它核酸样本,记录每个核酸样本与其所带的标签的对应关系;混合测序单元200,与样本标记单元100相连,用以混合获自样本标记单元的标签标记过的多个核酸样本,以及对混合核酸样本进行核酸序列测定,获得混合测序数据;比对筛选单元300,与混合测序单元200相连,用以实现将混合测序数据与参考序列比对,获得比对结果,以及从比对结果中筛选出与参考序列不完全匹配而且不匹配部分的长度不小于比标签小1bp的长度的读段;归类单元400,与样本标记单元100和比对筛选单元300相连,用以实现利用从比对筛选单元中筛选出的读段的信息和样本标记单元的对应关系,确定混合测序数据中读段源于的核酸样本。利用本发明的装置,可以实施完成本发明的方法。
根据本发明的一个具体实施方式,所说的装置的各单元还包含子单元,混合测序单元200还包括短读段去除子单元,可用于设置长度要求以去除混合测序数据中过短的读段;比对筛选单元300包括分类比对子单元,该子单元能够依据不匹配部分在读段中的位置对筛选出的读段进行分类,获得第一读段和第二读段,第一读段中的读段的两个末端都与参考序列不匹配,第二读段中的读段的两个末端中的一个与参考序列不匹配,以及比对第一读段中的每个读段中的与参考序列都不匹配的两个末端,去除两个末端不互相匹配的以及两个末端互相匹配的长度小于比A中标签小1bp的长度的读段。比如,标签的长度是7bp,读段两末端互相匹配的长度至少为6bp,该读段才会得以保留传至下个单元。
利用本发明提供的方法和/或装置,能够对多个样本的混合测序数据正确对应至各自的样本,提高了数据的利用率,利用对应到各自样本的数据检测各样本的变异情况,与Sanger测序检测各样本的变异结果完全一致。将各个检测数据需求量相对低、比如远低于测序平台通量的样本核酸混合建库、测序获得混合测序数据,能否进行正确归类将混合数据中的读段以及能够将多少的读段对应到相应的样本是进行后续检测的前提。本发明的方法中,还提供了自己设计的序列包括标签、标签和引物构成标签引物,能够区分大于等于标签数的样本,比如利用一对上下游引物序列都带同一个标签的引物对标记一个样本,使该样本的核酸扩增产物带一个特定标签序列以区别其它样本,再比如利用一对上下游引物序列带不同标签的引物对来标记一个样本,使该样本的核酸扩增产物带两个标签序列,这样只要所带的两个标签的任一个与其它样本的不同就能够区分其它样本核酸,再再比如利用多对引物标记一个样本,各对引物的上下游序列都带同一个标签或不同标签,各对引物之间带一样或不同标签,使该样本核酸扩增产物带一个或多个特定排列的标签,其中任一标签或者标签位置的不同足以区分不同样本。本发明提供的这组标签或者标签引物,是发明人通过多次试验、创造性劳动,设计考虑筛查各序列本身碱基组成、各种碱基比例以及序列间关系比如标签与标签之间、标签与引物之间、引物与引物之间、标签连接到引物后不同标签引物之间等等的关系,使这些序列的全部或者任一部分都能够在同一个反应体系中使用发挥作用。
下面详细描述本发明的实施例,所述实施例的示例在附图中示出,其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的,仅用于解释本发明,而不能理解为对本发明的限制。需要说明的是在本文中所使用的术语“第一”、“第二”等仅用于方便描述目的,而不能理解为指示或暗示相对重要性,也不能理解为之间有先后顺序关系。在本发明的描述中,除非另有说明,“多个”的含义是两个或两个以上。
除另有交待,以下实施例中涉及的试剂及仪器,都是常规市售产品,比如购自生命技术公司(life technologies)等。
实施例一混合样本文库构建、测序
96个人血液样本获自天津妇幼保健院。
将96个样本的核酸分别置于96孔板各孔中,在96孔板中对所有样本进行PCR标记反应:在提取血液样品中的DNA之后,将各个样本的DNA模板分别加入不同的孔内,并在每个孔内加入13对标签引物,标签引物是通过将标签序列连接到引物5’端预先制备的,记录样本与PCR引入的标签(barcode或index)的对应关系,然后置于PCR仪中进行PCR扩增,获得靶序列扩增产物,该多重PCR反应体系中各缓冲液组分及其量或比例的配置可参照已知的多重PCR反应体系,比如参考Hayden MJ,Nguyen TM,Waterman A,ChalmersKJ.2008.Multiplex-ready PCR:a new method for multiplexed SSR and SNPgenotyping.BMC Genomics.;doi:10.1186/1471-2164-9-80中的多重PCR反应体系配置进行,合成标签引物中的标签及引物分别选自表1和表2。然后根据标记混合扩增产物及纯化,再依据所采用的测序平台的文库构建说明进行混合文库构建,这边是根据Pronton的文库构建手册,包括末端修复和接头连接,这样,一个文库包含96个样本,接着利用AgilentBioanalyzer 2100检测文库片段大小和浓度,得到合格上机文库,接着按照Ion Proton半导体微芯片测序平台的上机操作流程对文库进行测序。
该实施例中由于扩增区域小,一个样本的数据量远低于Proton测序仪中一个通道(lane)的数据量,可参考了Illumina的多重(Multiplex)测序技术(Meyer M,KircherM.2010.Illumina Sequencing Library Preparation for Highly Multiplexed TargetCapture and Sequencing.Cold Spring Harb.Protoc.;doi:10.1101/pdb.prot5448),将其移植到Proton平台,比如在文库构建时利用连接标签接头、或者扩增引进新的标签来标记多个混合文库,可混合多个混合文库上机测序。图1示意的是利用本发明的标签序列扩增标记各个样本,混合各样本核酸建库,获得一个或多个文库,若是构建多个混合样本文库,在接头连接或者接头连接后的扩增中引入新的标签,这边引入的标签可以利用市售的公开可一起使用的标签,也可选自表1中除标记样本的标签序列以区分多个一起上机的混合文库。
实施例二混合测序数据的归类
获得下机数据之后,依据读段(reads)长度筛掉过短的reads,比如过滤掉小于50bp的reads,接着依据标签与样本的对应关系对reads进行归类。
构建特殊参考序列,这边的参考序列是由靶序列(目标扩增区域)构成的,参考靶序列区域的截取确定可依据实施例一中的13对引物比对到的参考基因组的位置来截取扩增区域确定。对混合数据进行归类,可以利用传统的方法,截取reads 5’端的7bp,即截取的长度等于标签的长度,将截取的这7bp与标记样本的标签序列对应,对应上的就将该reads归至该样本。
发明人发现从该半导体未芯片测序平台的下机数据直接与设计的扩增区域(捕获区域)序列构成的参考序列进行比对,由于该试验中文库的大小为200bp左右,扣掉两端接头即插入片段大约160bp,而测序获得的reads长度范围为50~300bp,大部分为120~180bp,大量的reads的长度接近或者大于等于插入片段的长度,由于该建库过程中没有打断扩增产物,相当于大量的reads的长度接近或者大于等于该由扩增区域构成的参考序列的长度,reads长度越长,比如测通的reads占的比例高,利于提高数据利用率。观测reads在这个参考序列上的比对位置和统计,发明人发现当比对上的起始位置位于该reads的第8个碱基时,该reads在其5’端具有一个完整的长度恰好为7bp的barcode;同样的,当比对的终止位置距离read 3’端的距离为7bp时,该read在3’端具有一个完整的barcode。经过统计,有64%的reads在5’端和3’端都具有完整的barcode,只在5’端有完整barcode的reads占到14%,只在3’端有完整barcode的reads占到12%。按照传统方法,只在3’端有完整barcode的reads是被放弃的。
这边先说明一定义:soft clip reads指reads比对回到基因组时,一条reads被切成两段,匹配到不同的区域,这样的reads叫做soft clip reads,一般是由于基因组发生某一段的缺失,或转录组的剪接,在测序过程中,这类reads横跨了缺失位点或剪接位点,在这边soft clip reads指只一部分比对到特殊参考序列的reads,未比对上的那段称为softclip。
基于这个发现,发明人提出另一套混合数据的归类方法,可以使归类对应的准确率提高至少12%,该归类对应reads至相应样本源的方法,也是首先利用lifetechnologies的比对软件tmap以默认参数设置将读段(reads)比对到以靶序列(扩增区域)构成的参考序列上,根据比对结果中的soft clip reads信息,对于每条soft clip reads,当其5’端或3’端有一个长度为7bp的soft clip发生时,认定这条read在5’端或3’端有完整的barcode序列。若read在5’端和3’端都具有完整的barcode序列,比较二者是否相同,不同则舍去该reads;若reads只在5’端或3’端具有barcode序列,则该序列就是该reads的barcode序列。基于之前已记录的barcode与样本的对应关系,可以得知该barcode序列对应的样本,将该reads的序列和质量值中截去barcode对应的部分,剩余的部分归入到该样本的数据中。在这一步,数据利用率能够高达90%。
实施例三样本核酸变异检测
对混合测序数据正确归类后,针对每个样本单独进行变异检测分析,使用到的软件是life technologies提供的比对软件tmap和变异检测软件tvc。这些软件都可以在proton测序仪的随机服务器上运行,其默认参数都是针对人类基因组设置的,可以不用调整。参考tmap和Torrent Variant Caller(tvc)的说明文档,利用扩增区域信息和已正确归类的reads信息,就能完成比对和变异检测的分析。其中,变异检测的核心原理是利用贝叶斯推断(Shoemaker JS,Painter IS,Weir BS.1999.Bayesian statistics in genetics:Aguide for the uninitiated.Trends Genet15:354–358)。该算法结合比对结果和read上的碱基质量值信息,能够得出该点的基因分型,从而判断该基因是否发生了突变。
通过利用13对标签引物同时扩增,能够实现针对与耳聋相关的GJB2、GJB3、SLC26A4和12sRNA四个基因进行变异检测。包括上述96例样本的总共300例样本的检测结果,与Sanger测序检测结果的一致性达到了100%,本方法的部分样本的检测结果如表3所示,表3中因为人是二倍体,所以检测结果用两个字母来表示,如果是SNP,就直接以碱基来表示,比如样本“14HL078963”的点突变“9G>A”的检测结果是“G.G”,说明该位点的一对等位碱基均未发生突变,即该位点没有发生SNP,针对插入缺失(indel)突变类型的检测,用R代表未突变,用V代表突变,Sanger测序检测这部分样本的突变,结果同表3。同时,本方法所需时间一般为2~3天,相比较与质谱法缩短了2天,而相对于全外显子组测序,本方法的测序成本仅为其1%,并且,本方法一次上机可以同时检测500份样本甚至更多,极大地提升了检测通量。
表3