基因组组装数学建模

基因组组装

摘要

基因组测序是生物信息学的核心，有着极其重要的应用价值。新的测序技术大量涌现，产生的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需要过滤掉。