基因组测序的数学模型分解

基因组组装

摘要

基因组测序是生物信息学的核心，有着极其重要的应用价值。新的测序技术大量涌现，产生的reads长度更短，数量更多，覆盖率更大，能直接读取的碱基对序列长度远小于基因组长度。所以测序之前DNA分子要经过复制若干份、随机打断成短片段。要获取整个DNA片段，需要把这些片段利用重合部分信息组织连接。如何在保证组装序列的连续性、完整性和准确性的同时设计耗时短、内存小的组装算法是本题的关键。

本文建立改进后OLC算法模型。该模型首先使用了特定的编码规定，通过C++程序对庞大的数据先后进行十进制和二进制的处理，不改变数据准确性的前提下尽可能减小内存和缩短计算机操作时间，并引入解决碱基识别错误问题的一般思路消除初始reads中的碱基错误。然后通过深度优先算法，设定适当的阈值，找出具有重叠关系的碱基片段并形成一有向赋权图，其中点是碱基片段，边代表具有重叠关系，权值代表片段重叠的多少，将问题转化为图论中寻找最大赋权通路的问题，从而对OLC算法进行改进，采用图论的方法更直观和更具操作性的解决DNA的拼接问题，从而对OLC算法进行改进。最后再根据OLC算法对Hamilton 路径进行拼接，生成共有序列，通过多序列比对等方法，获得最终的基因组序列。

关键词：基因组测序 OLC算法深度优先算法Hamilton路径

一问题的重述

1.1 问题背景

快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。

1.2 问题提出

确定基因组碱基对序列的过程称为测序。目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为

ATACCTT GCTAGCGT

GCTAGCGT AGGTCTGA

则有可能基因组序列中包含有ATACCTT GCTAGCGT AGGTCTGA这一段。

由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。

利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长（reads）。基因组复制份数约为50–100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。具体解决问题如下：

(1)建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。

(2)现有一个全长约为120,000个碱基对的细菌人工染色体，采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度约为70×，即基因组每个位置平均被测到约70次。试利

用你的算法和程序进行组装，并使之具有良好的组装效果。

二问题分析

2.1 问题一分析

鉴于现代测序技术的不完备性，各种基因组测序技术还有待改进和发展的空间，本文尝试性的建立数学模型，一方面对经典的OLC（overlap-layout-consensus）算法进行改进和发展，另一方面对现代测序技术提供参考和见解。对于基因组测序问题，本文采用图论的方法更直观和更具操作性的解决DNA的拼接问题。为了较好地解决测序中可能出现的个别碱基对识别错误，本文首先引入解决碱基识别错误问题的一般思路。鉴于OLC技术需要对碱基片段进行两两配对寻找重叠的碱基片段所造成的时间度复杂问题。本模型使用了特定的编码规定，通过C++程序对庞大的数据先后进行十进制和二进制的处理，使得不改变数据准确性的前提下大大降低了内存和缩短计算机操作时间。本模型首先通过深度优先算法，设定适当的阈值，找出具有重叠关系的碱基片段并形成一有向赋权图。其中点是碱基片段，边代表具有重叠关系，权值代表片段重叠的多少。这样问题将转化为图论中寻找最大赋权通路的问题。

2.2 问题二分析

基于问题一建立的模型，代入数据进行验算。

三模型假设

(1)假设测序过程中没有其他因素的干扰；

(2)假设题目所给定的序列相对位置的碱基全部遵循GU-AC法则；

(3)假设题目中所有的序列都是正常可判别的序列，没有出现序列的基因突变等情况；

(4)假设一个完整基因组，打断成500bp的片段是随机的；

(5)假设基因组每个位置被测到的几率是等可能的；

(6)所有片段上的碱基都已经被识别出来，不存在未知碱基。

四符号说明符号意义

reads 利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长

contig由reads经过一定算法拼接产生3kb～10Mb以内的一些基因组片段k-mer长度为k的一段DNA片段

quality 每一个reads都含有一个质量值，该值能反映该reads的正确率。质量值越高，reads的正确率越高

五模型建立及求解

5.1 数据预处理

5.1.1 数据简化处理

由于基因组进行编码的时候信息量非常的巨大，而且本文采用的数学模型需要对待定的所有reads进行两两的配对，以此确定无向图。若采用字符串的存储方式，显然会造成内存空间的大量消耗，甚至内存耗尽。为此，必须寻找其他的的存储方式，以达到降低内存空间消耗的目的。算法采用一套编码规则，将字符 A 编码为 00，字符 T 编码为 11，字符 G 编码为 01，字符 C 编码为 10。为便于研究，将二进制数转再化为对应的十进制数，这样就能大大的减少数据庞大给计算机运行和计算带来的难度，如图1

图1 编码规则

5.1.2 消除初始reads中的碱基错误

(1)收集的大量资料表明，测序数据中会有许多全A或者基本上全是A的reads，这些数据很可能是Solexa测序过程中的人工数据，需要去除。方法为：设定A的含量阈值为0.9，过滤掉含量大于等于0.9的reads。

(2)测序数据中含有一些未知的碱基，通常用“N”或“.”表示，其对拼接有不利的影响，因此含有未知碱基的read需要过滤掉。

5.1.3 序列片段中错误碱基的修正

本模型建立在传统测序技术中的OLC（overlap-layout-consensus）算法的改进上，由于现代测序技术并不完美，在测序前要通过 PCR 手段对待测片段进行扩增，从而增加了测序的错误率。在测序模型建立之前，为了降低PCR手段扩增带来的错误。有必要对 reads 数据进行预处理，修正 reads 中测序错误的碱基从而提高 DNA 序列拼接的效果。以下将引用常用的一种修正序列片段中错误碱基的方法。

由于基因组中每个位置进行测序的次数可能不止一次,每个位置的碱基在测序得到的序列片段集合中出现的期望次数为序列片段集合的覆盖率,因此在序列片段集合中可能存在多条在某一区域重叠的序列片段,如图 2 所示。

图2 序列片段集合中可能存在多条在某一区域重叠的序列片段

基于这个事实,当某个公共序列 U 达到一定的长度,并且序列片段集合中包含该公共序列的序列片段达到一定的数目时,我们可以认为该公共序列 U 是从基因组G 的某一个区域测序得到的,并且序列片段集合中所有包含该公共序列 U 的序列片段都是从该区域附近的某一个位置开始测序得到的。我们可以对紧跟在满足上述条件的公共序列后面的序列进行多序列比对,以此来修正序列片段中的错误碱基。图 3 是修正序列片段中错误碱基的一个简单、直观的例子,我们可以看到,通过这种方法第二条序列片段的倒数第四个碱基 C 被改为 G,最后一条序列片段第 19 个位置缺失的碱基 G 也被补上了。

图3 修正序列片段中错误碱基过程

5.2 基于OLC 策略及改进的深度优先算法对问题一模型的建立

针对 Sanger 测序技术产生的长度较长、错误率较低的序列片段,人们进行了广泛的研究，其中大部分技术都是采用基于 Hamilton 路径的算法实现的。本文基于哈密顿路径问题建立数学模型，使得传统的OLC 测序算法达到更优。

以全部待拼接的reads 为节点，给定一个适当的阈值λ，则用节点间的连线代表reads 点之间有重叠部分，且这个重叠部分大于阈值λ。那么就把DNA 测序问题转化为一般图论问题。

对于可定图），，（W E V G =，V 、E 和W 分别代表图的顶点、边和边上的权的

集合。其中，W 表示重叠部分）

（λ≥i W ，以待定reads 为始点，寻找一条通路，使得有且只有一次经过尽可能多的点并使得权值最大，即哈密顿通路。此时DNA 测序问题将转化为图论中对于给定图求赋权值最大的所有哈密顿通路问题。其中哈密度通路的条数为contig 条数，权值最大的哈密顿通路为最长contig 。如图4所示：

图4 重叠关系图

该算法的核心是构建重叠关系图对于处理 Sanger 数据或者 454、Ion Torrent 数据具有优势。主要包括 2个步骤：

（1）处理

本模型首先需要对待定的所有reads进行两两的配对，当两对reads的重叠部分超过某个设定的阈值的时候，说明这一对reads有联系。针对该问题，我们采用改进的深度优先算法把有联系的reads点连接起来，从而得到一个复杂的有向赋权图G。首先介绍改进的深度优先算法的基本思路：

1）把一个具体的问题抽象成了一个图论的模型有向图状态对应着结点，状态之间的关系（或者说决策方案）对应着边；

2）从当前的某个节点开始历遍所有的点，去掉所有低于阈值的路，构成一个新的有向赋权图；

3）在各个阶段尝试方案时，采取的是穷举的思想。

根据该算法，我们定义每两条reads重叠部分的碱基数量为权，两个reads之间重叠越多则两个节点之间的权越大。

（2）拼接

该步骤是将第一步中全局比对得到的覆盖信息组装并构建一张重叠关系图。根

据节点处数的大小，可以判断该链接是否为可靠链接。计算机根据全图的节点，计算 Hamiltonian通路。所有通路上的reads串联就构成了一条完整的链。

1）首先取任一条reads为contig，接着寻找与该reads的两端含有重叠区域的reads，则可能存在无数条这种reads，那么我们需要先设定一个阈，当重叠区域的碱基数量超过阈值时才能将其视为满足条件；

2）排列reads，确定reads之间的相对位置，建立overlap图，然后分析overlap，获得历遍整个图的最佳近似路径，找到Hamilton 路径；

3）生成共有序列，通过多序列比对等方法，获得最终的基因组序列。

六模型评价

6.1优点：

(1)本模型的算法容易推广到实际的基因组组装中，具有一定的实际应用价值；

(2)利用并行技术提高数据处理效率。新一代测序平台产生的序列片段数据量庞大,要将本文的数据处理方法用于实际测序工程需要解决速度和内存两方面的问题。利用并行技术,实现算法在集群环境的并行化将是不错的解决方案。

6.2缺点：

(1)忽略了碱基存在的内环境因素及其生化结构的影响；

(2)在实际中，基因组组装是一个复杂的数学问题，存在着大量的不确定性；

(3)一些新的想法缺乏足够的理论根据，所以有些问题的解决带有一定的主观性

(4)该模型处理重复片段的能力较弱。当有reads拼接失败时，意味着小于reads 长度的重复片段被检验出来，但无法处理大于reads长度的重复片段。

6.3模型的改进

对于处理重复片段弱的情况，尝试使用覆盖度分析，建立相应对的概率模型，必要时使用配对reads信息，来发现大于reads长度的重复片段。

七模型推广

运用本模型可以在个人计算机上解决相对不大的数据，也可以在内存较大的计算机上解决数据量较大的基因组组装问题，同时此种算法模型还可以作为基础模型，

解决一些降低重复数据的生产生活问题。

八参考文献

[1]陈传艺.针对新一代测序技术的序列拼接算法研究[D].福建农林大学,2012.

[2]赵洁,赵志军.新一代测序技术及其应用[J].白求恩军医学院学报,2012,04:344-345.

[3]韩东涛.基于概率模型的基因组从头测序算法研究[D].哈尔滨工业大学,2012.

[4]蔡毅,骆志刚.DNA序列拼接算法分析及并行化探讨[A].中国电子学会信息论分会、北京邮电大学研究生院.2007北京地区高校研究生学术交流会通信与信息技术会议论文集（上册）[C].中国电子学会信息论分会、北京邮电大学研究生院:,2008:6.

[5]张书翠. 基于高通量测序的Klebsiella pneumoniae基因组拼接的研究[D].上海师范大学,2013.

九附录

1.提取数据及处理数据的C++程序

程序1

#include

#include

#include

using namespace std;

void main()

{

ifstream in("F://McMc_BAC_1.fq.gz.clean.dup.clean");//打开文件1

if(!in){

cout<<"File open error!"<

return;

}

ifstream fin("F://McMc_BAC_2.fq.gz.clean.dup.clean");//打开文件2

if(!fin){

cout<<"File open error!"<

return;

}

ofstream f("F://read1碱基序列集.txt");//打开写入文件，若文件不存在则新建if(!f){

cout<<"File write error!"<

return;

}

int i=0;string s;

while(getline(in,s))//着行读取文件1数据并存于s中，直至数据全部读取

{i=i+1;

if((i+2)%4==0)

f<

}

ofstream of("F://read2碱基序列集.txt");//打开写入文件，若文件不存在则新建if(!of){

cout<<"File write error!"<

return;

}

while(getline(fin,s))//着行读取文件2数据并存于s中，直至数据全部读取

{i=i+1;

if((i+2)%4==0)

of<

}

}

2.OLC算法的MATLAB程序

程序1

%%

% 求那些数据相同

% 其中：

% j表示有n-j个连续碱基相同

% i，k表示第i个与第j个数据相同

% 假定k=4

% 计算相同项的位置

clear,clc

fid=fopen('E:\123.txt','r');

data=fread(fid);

[m,n]=size(data);

data=reshape(data,[90,85]);

data=data';

K=84;

[m,n]=size(data);

n=n-2;

d1=[];

d2=[];

d3=[];

for j=1:K-1

for i=1:m

for k=i+1:m

index1=[];

for t=j+1:n

index=abs(data(i,t)-data(k,t-j));

index1=[index1,index];

end

index2=sum(index1);

if index2==0

d1=[d1,j];

d2=[d2,i];

d3=[d3,k];

end

end

end

end

d=[d1;d2;d3];

c=d';

xlswrite('E:\dnanew.xlsx',c);

fclose(fid);

程序2

%%

% 计算最佳路径

clear,clc

c=xlsread('E:\dnanew.xlsx');

c=c';

j=c(1,:);

i=c(2,:);

k=c(3,:);

n=length(j);

a=[];

b0=0;

a0=[];

for t=1:n

a=i(t);

b=88;

for tt=t:n

if k(t)==i(tt)

k(t)=k(tt);

i(t)=i(tt);

a=[a,i(tt)];

b=b+j(tt);

end

end

if b0<=b

b0=b;

a0=a;

end

end

a0,b0

l=a0';

xlswrite('E:\a0.xlsx',l); % 将所得的数据储存在a0.xlsx中

程序3

%%

% 此程序解决基因链组装

clear,clc

fid=fopen('E:\123.txt','r'); % 打开E盘123.txt文件，其中储存基因链组data=fread(fid); % 将fid转换成数据（更具ASCII码）

[m,n]=size(data);

data=reshape(data,[90,85]); % 重塑数据成90*85的矩阵（多了回车和换行字符）

data=data';

c=xlsread('E:\dnanew.xlsx'); % 读取前期形成的dnanew.xlsx，其中储存路径c=c';

j=c(1,:);

a=xlsread('E:\a0.xlsx'); % 读取前期形成的a0.xlsx，其中储存具体路径a=a';

dt=data(a(1),1:88);

for i=2:length(a)

k=89-j(a(i));

dt=[dt,data(i,k:88)];

end

dt;

dlmwrite('E:\final.txt',char(dt),''); % 将获得的长链储存在final.txt中fclose(fid);