CN109817280B

CN109817280B - 一种测序数据组装方法

Info

Publication number: CN109817280B
Application number: CN201910146649.7A
Authority: CN
Inventors: 马丰收; 张艺; 何飞; 刘洋
Original assignee: Genergy Bio Technology Shanghai Co ltd
Current assignee: Genergy Bio Technology Shanghai Co ltd
Priority date: 2016-04-06
Filing date: 2016-04-06
Publication date: 2023-04-14
Anticipated expiration: 2036-04-06
Also published as: CN106021978B; CN106021978A; CN109817280A

Abstract

本发明涉及一种测序数据组装方法，使用光学图谱平台Irys得到基因组装文件；同时，拿到NGS的scaffold文件：fai文件；数据预处理：通过设定阈值，过滤掉可信度低的比对结果，合并cmap文件，排序，计算N50；组装效果统计：统计BioNano与NGS比对结果，包括BioNano的contig与NGS的scaffold长度、个数以及总量；根据BioNano的contig与NGS的scaffold之间的网络拓扑关系，分类分析组装出的新的contig长度与scaffold长度。可以辅助基因组组装，明显提高物种的基因组装效果。

Description

一种测序数据组装方法

技术领域

本发明涉及生物信息学，特别涉及主要应用于辅助de novo测序数据组装与结构变异的检测。

背景技术

基因组de novo测序即基因组从头测序，是指对基因组序列未知或没有近缘物种基因组的某个物种的全基因组序列的测序。然后利用生物信息学手段对测序序列进行拼接、组装和注释，从而获得该物种完整的基因组序列图谱。目前，最常用的方法是二代测序(Next Generation Sequencing,NGS)，但是NGS方法会使大量重复元件和结构变异的信息丢失，所以组装完整的基因组图谱就成为一个难题。

BioNano Genomics公司拓展了纳米通道技术，并把它发展为一个灵活的具有高分辨率和极长单分子DNA的光学图谱平台Irys，使得常规和精确检测基因组结构变异和完成基因组组装成为可能(Shelton J M,Coleman M C,Herndon N,et al.Tools andpipelines for BioNano data:molecule assembly pipeline and FASTA superscaffolding tool[J].BMC genomics,2015,16(1):1.)。通过捕捉高分辨率的极长的DNA分子，Irys系统提供一个基因组图谱，可以观察任何生物的基因组包括人类(Pendleton M,Sebra R,Pang A W C,et al.Assembly and diploid architecture of an individualhuman genome via single-molecule technologies[J].Nature methods,2015.)。该图谱可以洞察引起表型变异的新的基因组结构变异(Cao H,Hastie A R,Cao D,et al.Rapiddetection of structural variation in a human genome using nanochannel-basedgenome mapping technology[J].GigaScience,2014,3(1):1-11.)；利用全面的基因组视角更快更完整的实现从头组装；利用单分子成像打开真实的基因组生物学与一系列应用；避免扩增和剪切带来的错误，在单分子水平解剖复杂的混合物，而不是在一个平均测量值上面丢失罕见的变异；高性价比和高通量的Irys系统更全面和更准确的刻画了整个基因组范围内的基因组事件；更全面的了解整个基因组的事件。

发明内容

本发明的目的在于提供一种基于光学图谱平台Irys的一种denovo测序数据组装方法，计算得到更长、更准确的contig与scaffold。

为解决上述技术问题，本发明的实施方式提供了一种基于光学图谱平台Irys的一种denovo测序数据组装方法，步骤为：

第一步，使用光学图谱平台Irys得到基因组装文件：xmap文件、q.cmap文件、r.cmap文件；同时，拿到NGS的scaffold文件：fai文件；

第二步，数据预处理：通过设定阈值，过滤掉低可信度的比对结果，合并cmap文件，排序，计算N50；

第三步，组装效果统计：统计BioNano与NGS比对结果，包括BioNano的contig与NGS的scaffold长度、个数以及总量；

第四步，根据BioNano的contig与NGS的scaffold之间的网络拓扑关系，分类分析组装出的新的contig长度与scaffold长度。

第二步数据预处理具体步骤为：i)筛选fai文件：按照xmap文件中QryContigID提取fai文件中的query名称，作为xmap的第11列；ii)筛选xmap文件：设定阈值，筛选出所有符合条件的比对文件；iii)分别对q.cmap文件、r.cmap文件去重，然后与xmap文件合并，并按照BioNano的contig ID排序，并计算N50。

第四步分类分析是指采用无向拓扑图显示两种方法中参与比对contig的大小与网络关系，BioNano的contig与NGS的scaffold按照一一对应、一对多、多对一的关系进行分类分析；一一对应，计算contig；一对多，计算scaffold的长度。

所述的一一对应，计算contig具体分为四类：i)BioNano的contig片段包含在NGS的scaffold中，表明组装出来的contig长度没有NGS的scaffold长度长；故此时gap为0，组装的contig记为NGS的scaffold长度Contig_length，NGS的scaffold长度为QL

Contig_length＝QL；

ii)NGS的scaffold片段包含在BioNano的contig中，说明两个酶切位点之间的距离，也就是记录的contig之间存在gap，则计算gap长度；假设两个酶切位点分别为RS、RE，记录的contig长度为RL，则gap长度，记为Gap_size，则

Gap_size＝RS+RL-RE；

iii)BioNano的contig长度与NGS的scaffold长度相等，此时gap为0，定义此时组装的contig长度记为Contig_length为两者的平均值，假设NGS的scaffold起始位置为QS、QE，长度为QL，则有

RE-RS-RL＝QE-QS-QL，

Contig_length＝(RL+QL)/2；

iiii)BioNano的contig长度与NGS的scaffold分为正方向比对与负方向比对，分别用“+”与“-”区分；正方向比对，此时的contig长度为

Contig_length＝RL-RE+QE

或者

Contig_length＝RS+QL-QS

Gap为RL-RE或者RS；

BioNano的contig长度与NGS的scaffold负方向比对，此时，contig长度为

Contig_length＝RL-RE+QL-QE

或者

Contig_length＝RS+QS

Gap为QE或者QL-QS。

所述的一对多，计算scaffold的长度是指提取一条BioNano的contig与两条NGS的scaffold比对上的contig，并按照拼接方向，分为“++”、“+-”、“-+”、“--”4种情况组装出更长的contig，也就是scaffold；

当拼接方向为“++”时，定义计算公式如下：

scaffold_length＝(RS2-RS1)+QS1+(QL2-QS2)；

当拼接方向为“+-”时，定义此时组装出来的新的scaffold长度如下：

scaffold_length＝(RS2-RS1)+QS1+(QS2-QE2)；

当拼接方向为“-+”时，新的scaffold长度为：

scaffold_length＝(RS2-RS1)+(QL1-QS1)+(QL2-QS2)；

当拼接方向为“--”时，有

scaffold_length＝(RS2-RS1)+(QL2-QS1)+QE2；

RS1：BioNano方法中contig1的开始位置；

RE1：BioNano方法中contig1的末端位置；

RS2：BioNano方法中contig2的开始位置；

RE2：BioNano方法中contig2的末端位置；

QS1：NGS方法中contig1的开始位置；

QS2：NGS方法中contig1的末端位置；

QE1：NGS方法中contig2的开始位置；

QE2：NGS方法中contig2的末端位置；

RL：BioNano方法中contig1+contig2长度；

QL1：NGS方法中contig1的长度；

QL2：NGS方法中contig2的长度。

本发明在光学图谱平台Irys的基础上，通过分析BioNano的contig与NGS的scaffold之间的网络拓扑关系，提出一种新的de novo测序数据组装方法，计算得到更长、更准确的contig与scaffold。

本发明的基于光学图谱平台Irys的统计方法，可以辅助基因组组装，明显提高物种的基因组装效果。

附图说明

图1是BioNano的contig与NGS的scaffold网络关系拓扑图。圆点代表BioNano的contig，三角形代表NGS的scaffold。

图2是BioNano的contig片段包含于NGS的scaffold。

图3是NGS的scaffold片段包含于BioNano的contig。

图4是BioNano的contig长度与NGS的scaffold长度相等。

图5是BioNano的contig长度与NGS的scaffold正方向比对。

图6是BioNano的contig长度与NGS的scaffold负方向比对。

图7一条BioNano的contig与两条NGS的scaffold“++”方向比对。

图8一条BioNano的contig与两条NGS的scaffold“+-”方向比对。

图9一条BioNano的contig与两条NGS的scaffold“-+”方向比对。

图10一条BioNano的contig与两条NGS的scaffold“--”方向比对。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明的各实施方式进行详细的阐述。然而，本领域的普通技术人员可以理解，在本发明各实施方式中，为了使读者更好地理解本申请而提出了许多技术细节。但是，即使没有这些技术细节和基于以下各实施方式的种种变化和修改，也可以实现本申请各权利要求所要求保护的技术方案。

首先介绍本发明中用到的几个概念：

1.read:测序过程中，一个DNA分子先经过克隆形成若干个拷贝，然后这些拷贝被打碎成若干条短的，可以直接测序的片段，每一条片段称作一个"read"，测序仪产生的即是read的集合。

2.Contig：拼接软件基于reads之间的overlap区，拼接获得的序列称为Contig(重叠群)。

3.Scaffold：基因组de novo测序，通过reads拼接获得Contigs后，往往还需要构建454Paired-end库或Illumina Mate-pair库，以获得一定大小片段(如3Kb、6Kb、10Kb、20Kb)两端的序列。基于这些序列，可以确定一些Contig之间的顺序关系，这些先后顺序已知的Contigs组成Scaffold。

4.Contig N50：Reads拼接后会获得一些不同长度的Contigs。将所有的Contig长度相加，能获得一个Contig总长度。然后将所有的Contigs按照从长到短进行排序，如获得Contig 1，Contig 2，Contig 3…………Contig 25。将Contig按照这个顺序依次相加，当相加的长度达到Contig总长度的一半时，最后一个加上的Contig长度即为Contig N50。

5.Scaffold N50：Scaffold N50与Contig N50的定义类似。Contigs拼接组装获得一些不同长度的Scaffolds。将所有的Scaffold长度相加，能获得一个Scaffold总长度。然后将所有的Scaffolds按照从长到短进行排序，如获得Scaffold 1，Scaffold 2，Scaffold3…………Scaffold 25。将Scaffold按照这个顺序依次相加，当相加的长度达到Scaffold总长度的一半时，最后一个加上的Scaffold长度即为Scaffold N50。

本发明的目的是在光学图谱平台Irys分析基础上，提出一种新的de novo测序数据组装方法，更快的组装出更长、更准确的contig与scaffold。

本发明提出的方法，包括以下步骤：

1Irys结果数据的简要说明

将物种的基因测序数据导入Irys软件后，会得到以下组装的质量报告与cmap数据文件。报告

测序数据质量报告

测序后过滤得到的数据质量报告，包含大于100kb、150kb、180kb、250kb、500kb以及150-200kb的reads个数、数据量、bin的质量分数等信息。

比对结果报告

拼接组装质量报告，比对统计结果包括10-500kb、100-500kb、150-500kb、200-500kb、250-500kb、以及大于500kb的bin的个数、数据量、质量分数、置信度以及真假阳性率等。

数据

(一)cmap文件

1.q.cmap文件

q.cmap文件包含Query的比对ID(Qry.CMapId)、Query的contig长度(Qry.ContigLength)、Query的酶切位点个数(Qry.NumSites)、Query的酶切位点标号(Qry.SiteID)、Query的酶切位点位置(Qry.Position)、Query的标准差(Qry.StdDev)与Query的酶切位点覆盖度(Qry.Coverage)等信息。

2.r.cmap文件

r.cmap文件中包含reference的cmap编号(Ref.CMapId)、reference的contig长度(Ref.ContigLength)、reference的酶切位点个数(Ref.NumSites)reference的酶切位点标号(Ref.SiteID)、reference的酶切位点位置(Ref.Position)、reference的酶切位点标准差(Ref.StdDev)以及reference的酶切位点覆盖度(Ref.Coverage)等。

3.xmap文件

此文件中每列所含信息命名以及如下：

比对事件的ID(XmapEntryID)、query的contig ID(QryContigID)、referencecontig ID(RefcontigID)、比对在query的起始位置(QryStartPos)、比对在query的终点位置(QryEndPos)、比对在reference的起始位置(RefStartPos)、比对在reference的终点位置(RefEndPos)、比对发生的方向(Orientation)、比对的可信度(Confidence)以及比对的详细match信息，M表示match，I表示insertion，D表示deletion(HitEnum)。

(二)NGS的scaffold文件(fai文件)

含有杂乱无章的scaffold信息，例如比对scaffold的名称(Qry.name)、scaffold长度(Qry.length)等，数据量为353Mb,总长度10,542bp，N50为721kb。

2数据预处理及统计

2.1筛选fai文件

按照xmap文件中QryContigID提取fai文件中的query名称，作为xmap的第11列。

2.2筛选xmap文件

为了过滤掉低可信度的比对结果，需要设定阈值，筛选出所有符合条件的比对文件。

这里，假设阈值为10，置信度超过(包含)10的比对，被认为是可靠的。

2.3合并

分别对q.cmap文件、r.cmap文件去重，然后与xmap文件合并，并按照BioNano的contig ID排序，方便统计比对序列的N50。

2.4比较N50

统计BioNano的contig与NGS的scaffold长度、个数以及总量，如表1所示。

表1 BioNano与NGS比对结果统计

由表1可以看出，BioNano比对上的contig比率为44.60％，要远大于NGS得到的scaffold；而总体参与比对的contig比率与NGS相差不大。

3BioNano的contig与NGS的scaffold网络拓扑图

表1此表显示了BioNano的contig与NGS的scaffold的总体分布信息，但是并不能表明每个contig的顺序、位置信息。如何组装出更长的contig才是本发明最关心的问题。首先，用无向拓扑图显示两种方法中参与比对contig的大小与网络关系，如图1所示。从图中可以看到几乎每个NGS的scaffold都被覆盖到，并且两者之间的关系十分复杂，存在BioNano的contig与NGS的scaffold一一对应、一对多、多对一等关系。

3.1一一对应，计算contig

首先，考虑最简单的情况，contig与NGS的scaffold一一对应。按照这种关系从图3中提取所有符合条件的contig，将其可以分为4类：

(一)BioNano的contig片段包含在NGS的scaffold中，如图2所示；

此种情况表明组装出来的contig长度还没有NGS的scaffold长度长；故此时gap为0，组装的contig记为NGS的scaffold长度(记为QL)。

Contig_length＝QL

(二)NGS的scaffold片段包含在BioNano的contig中，如图3所示；

这时，说明两个酶切位点之间的距离，也就是记录的contig之间存在gap(也就是两端空白)，并且可以计算gap长度。例如，假设两个酶切位点分别为RS、RE，记录的contig长度为RL，则gap长度(记为Gap_size)为

Gap_size＝RS+RL-RE

(三)BioNano的contig长度与NGS的scaffold长度相等，如图4所示；

显然，此时gap为0，定义此时组装的contig长度(记为Contig_length)为两者的平均值。例如，假设NGS的scaffold起始位置为QS、QE，长度为QL，则有

RE-RS-RL＝QE-QS-QL

Contig_length＝(RL+QL)/2

(四)其他情况

一般分为正方向比对与负方向比对，分别用“+”与“-”区分。

BioNano的contig长度与NGS的scaffold正方向比对，如图5所示。

此时的contig长度为

Contig_length＝RL-RE+QE

或者

Contig_length＝RS+QL-QS

Gap为RL-RE或者RS。

b)BioNano的contig长度与NGS的scaffold负方向比对，如图6所示。

此时，contig长度为

Contig_length＝RL-RE+QL-QE

或者

Contig_length＝RS+QS

Gap为QE或者QL-QS。

按照上述四种情况计算出组装的contig长度，可以统计出最终组装的contig总长度、片段数、平均长度以及平均gap长度、contig50等来评估组装效果的好坏。统计结果如表2所示。

表2全基因组拼接contig统计表

组装指标	一一对应
		总数据量(Mb)	31.54
contig数	35
		平均contig(Mb)	0.9
平均Gap大小(Kb)	250.552

4.3.2一对多，计算scaffold的长度：

同样地，也可以从图1中提取一条BioNano的contig与两条NGS的scaffold比对上的contig，并按照拼接方向，分为“++”、“+-”、“-+”、“--”等4种情况组装出更长的contig，也就是scaffold。首先，定义下列公式中将要出现的参数含义：

RS1：BioNano方法中contig1的开始位置；

RE1：BioNano方法中contig1的末端位置；

RS2：BioNano方法中contig2的开始位置；

RE2：BioNano方法中contig2的末端位置；

QS1：NGS方法中contig1的开始位置；

QS2：NGS方法中contig1的末端位置；

QE1：NGS方法中contig2的开始位置；

QE2：NGS方法中contig2的末端位置；

RL：BioNano方法中contig1+contig2长度；

QL1：NGS方法中contig1的长度；

QL2：NGS方法中contig2的长度；

当拼接方向为“++”时，如图7所示，(第一条虚线表示BioNano的contig，下方两条虚线表示NGS的scaffold，下同)。

从图9中可以看出，BioNano的contig长度正好与NGS的两条scaffold中间重叠，组成一条新的scaffold，长度即为3条片段的长度减去重叠部分。这里，定义计算公式如下：

scaffold_length＝(RS2-RS1)+QS1+(QL2-QS2)

同理，可以定义如图8-10所示的新的scaffold长度。当拼接方向为“+-”时(如图8所示)，定义此时组装出来的新的scaffold长度如下：

scaffold_length＝(RS2-RS1)+QS1+(QS2-QE2)；

当拼接方向为“-+”时，新的scaffold长度为：

scaffold_length＝(RS2-RS1)+(QL1-QS1)+(QL2-QS2)；

当拼接方向为“--”时，有

scaffold_length＝(RS2-RS1)+(QL2-QS1)+QE2。

按照上述四种情况计算出组装的scaffold长度，统计scaffold总长度、数量、scaffold50等来评估组装效果的好坏。统计结果如表3所示。

表3全基因组拼接scaffold统计表

4.4本方法在其他的物种的应用

基于光学图谱平台Irys的统计方法，可以辅助基因组组装，明显提高物种的基因组装效果。以拟南芥为例，Tair10～120Mb，数据量48G～400x(>＝150kb)。组装结果如表4所示。

表4拟南芥基因组装结果统计表

N50从1.526Mb提高到4.229Mb，增加了177％；Scaffold+Contig数量从688缩减到595，减少了13.5％；最长的scaffold从5.3Mb提高到8.89Mb，增加了67.7％。

为了更好的说明本方法的有效性与实用性，分别对不同的物种使用，并得到如表5的统计结果。

表5基于光学图谱平台Irys的基因组装方法的优势在其他物种上的验证

项目	前后对比
		昆虫，数据量50x覆盖	N50长度提升80.5％，最长scaffold提升82.5％
动物，数据量100x覆盖	N50长度提升151％，最长scaffold提升94.7％
		植物，数据量100x覆盖	N50长度提升85.7％，最长scaffold提升60.4％

上面各种方法的步骤划分，只是为了描述清楚，实现时可以合并为一个步骤或者对某些步骤进行拆分，分解为多个步骤，只要包含相同的逻辑关系，都在本专利的保护范围内；对算法中或者流程中添加无关紧要的修改或者引入无关紧要的设计，但不改变其算法和流程的核心设计都在该专利的保护范围内。

本领域的普通技术人员可以理解，上述各实施方式是实现本发明的具体实施例，而在实际应用中，可以在形式上和细节上对其作各种改变，而不偏离本发明的精神和范围。

Claims

1.一种测序数据组装方法，其特征在于，步骤为：

第二步，数据预处理：通过设定阈值，所述阈值为10，过滤掉低可信度的比对结果，合并cmap文件，排序，计算N50；

第四步，根据BioNano的contig与NGS的scaffold之间的网络拓扑关系，分类分析组装出的新的contig长度与scaffold长度；

其中，第二步数据预处理具体步骤为：i)筛选fai文件：按照xmap文件中QryContigID提取fai文件中的query名称，作为xmap的第11列；ii)筛选xmap文件：设定阈值，筛选出所有符合条件的比对文件；iii)分别对q.cmap文件、r.cmap文件去重，然后与xmap文件合并，并按照BioNano的contig ID排序，并计算N50；

其中，第四步分类分析是指采用无向拓扑图显示两种方法中参与比对contig的大小与网络关系，BioNano的contig与NGS的scaffold按照一一对应、一对多、多对一的关系进行分类分析；一一对应，计算组装的contig长度；一对多，计算组装的scaffold的长度；

所述的一一对应，计算组装的contig长度具体分为四类：

i)BioNano的contig片段包含在NGS的scaffold中，表明组装出来的contig长度没有NGS的scaffold长度长；故此时gap为0，组装的contig长度Contig_length记为NGS的scaffold长度，NGS的scaffold长度为QL；

Contig_length＝QL；

ii)NGS的scaffold片段包含在BioNano的contig中，说明记录的contig之间存在gap，则计算gap长度；设两个酶切位点分别为RS、RE，记录的contig长度为RL，则gap长度为Gap_size，则：

Gap_size＝RS+RL-RE；

iii)BioNano的contig长度与NGS的scaffold长度相等，此时gap为0，定义此时组装的contig长度记为Contig_length为两者的平均值，设NGS的scaffold起始和结束位置分别为QS、QE，长度为QL，则：

RE-RS-RL＝QE-QS-QL，

Contig_length＝(RL+QL)/2；

iiii)BioNano的contig长度与NGS的scaffold分为正方向比对与负方向比对，分别用“+”与“-”区分；正方向比对，此时组装的contig长度Contig_length为：

Contig_length＝RL-RE+QE

或者

Contig_length＝RS+QL-QS

Gap为RL-RE或者RS；

BioNano的contig长度与NGS的scaffold负方向比对，此时，组装的contig长度Contig_length为：

Contig_length＝RL-RE+QL-QE

或者

Contig_length＝RS+QS

Gap为QE或者QL-QS；

所述的一对多，计算组装的scaffold长度是指：

提取一条BioNano的contig与两条NGS的scaffold比对上的contig，并按照拼接方向，分为“++”、“+-”、“-+”、“--”4种情况组装出更长的contig，即为组装的scaffold；

当拼接方向为“++”时，组装的scaffold长度scaffold_length为：

scaffold_length＝(RS2-RS1)+QS1+(QL2-QS2)；

当拼接方向为“+-”时，组装的scaffold长度scaffold_length为：

scaffold_length＝(RS2-RS1)+QS1+(QS2-QE2)；

当拼接方向为“-+”时，组装的scaffold长度scaffold_length为：

scaffold_length＝(RS2-RS1)+(QL1-QS1)+(QL2-QS2)；

当拼接方向为“--”时，组装的scaffold长度scaffold_length为：

scaffold_length＝(RS2-RS1)+(QL2-QS1)+QE2；

RS1：BioNano方法中contig1的开始位置；

RE1：BioNano方法中contig1的末端位置；

RS2：BioNano方法中contig2的开始位置；

RE2：BioNano方法中contig2的末端位置；

QS1：NGS方法中contig1的开始位置；

QS2：NGS方法中contig1的末端位置；

QE1：NGS方法中contig2的开始位置；

QE2：NGS方法中contig2的末端位置；

QL1：NGS方法中contig1的长度；

QL2：NGS方法中contig2的长度。