基因问题算法.doc

基因的拼接与组装模型摘要本文研究的是基因的拼接与组装问题，建立了逐段延伸拼接模型、分离对位组装模型对测序结果进行拼接与组装，并利用序列相似度、先到先得原理对测序出现的误差和重复序列进行处理。对问题一：本文建立了逐段延伸拼接算法模型与分离对位组装算法模型来对读长序列拼接与组装。由于测序结果有两组读长序列，本文只对其中一组进行拼接，再用另一组与拼接结果对位组装。首先,从一组读长中任取一条序列作为种子，从剩余序列中检索与其尾部相同的序列，并拼接到序列末端成为新的粽子，再次检索与其尾部相同的序列，继续拼接，直到序列长度达到规定值，再从剩余序列选出一条序列作为种子进行拼接，直至序列全部成为种子为止；其次，用生物学知识将拼接得的序列的反向互补序列找出，将其逐个分离成与读长序列等长的序列，并在另一组读长中检索出所有与分离得的序列相同的序列，并放入对应位置，将重叠部分拼接；最后，将得到的碱基对片段拼接成一条完整的基因，既可完成序列的拼接与组装。对问题二：基于问题一的算法模型，对问题二中给出的read1和read2的读长序列进行拼接与组装。先将read1中的序列读入Matlab中，利用逐段延伸算法，将其拼接成为若干条长度为500的序列；再将拼接得到的每条长度为500的序列的反向互补序列求出，按照分离对位算法将其分离成为长度为88的序列413（=500-88+1）条，并记录下每条分离序列的位置；然后，在read2中检索出所有与分离得到的系列相同的序列，并将其放入对应的位置，将重叠部分拼接。再将剩下的序列按照相似度大小放入空余的位置；最后，将长度为500的片段拼接为一条完基因，既可完成组装。得到的基因中含有碱基对个数117786对。关键词：逐段延伸算法、分离对位算法、相似度一、问题的重述随着新一代测序技术的快速发展，对测序得到的读长序列的拼接与组装也变得非常困难。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。问题一：本文建立了数学模型，运用了逐段拼接算法和逐段分离组装算法，并编制相应程序，将读长序列组装成基因组。本文的算法和程序能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。问题二：本文对现有的一个全长约为120,000个碱基对的细菌人工染色体（BAC），根据测得的读长数据read1和read2，利用本文的算法和程序进行组装，并有良好的组装效果。二问题分析2.1对问题一的分析 由于新一代测序技术测出的读长序列不仅长度短、数量巨大，而且测序结果可能存在错误，这给拼接与组装带来了巨大的困难。因此，在拼接和组装过程中，算法不仅要能解决拼接速度和组装速度的问题以外，还要能有效的处理测序中出现的错误与重复序列，而传统的OLC算法和贪婪图法在处理这么巨大的数据是速度很慢，并且由于将数据用图像存储，占用内存空间较大。因此，本文在建立算法时优先考虑内存占用问题与拼接速度问题。 本文建立的逐段延伸拼接算法与分离对位组装算法避免了用图像储存数据，能节省一定的内存空间，而且拼接与组装速度相对较快，能较好的解决内存与速度的问题。 而在拼接过程中，可能检索出多条符合条件的序列或者检索到测序错误的序列，本文遵循先到先得原理与序列相似度原理，即先检索到的序列优先拼接，越相似的先拼接。2.2对问题二的分析 由于题中给出了DNA两条链上的读长序列数据，根据数据格式，本文利用Matlab只读取了第二行的数据，并将其转化为对应的ascii码矩阵，即转化为数值型矩阵，以便于比较和求相似度。然后根据问题一的逐段延伸拼接算法将read1中的读长序列拼接成为长度为500bp左右的长序列，并保存；再将这一系列的长序列的反向互补序列求出，将条逐个分离成为长度为88的序列，并记其序列位置为每个序列末端到长序列末端之间的碱基个数，最后一个分离的序列的位置为0；然后，从read2中检索出所有与分离的序列相同的序列放在对应的位置，将重叠部分拼接，再求出剩下序列与空余的分离序列的相似度大小，按照相似度大小将其放在空余的部分，直至填满所有的长序列，即可完成组装；最后，将这些组装好的长序列拼接成为一条完整的基因。3、 符号说明C1i表示第一组读长序列的第i条读长转化的矩阵Cij表示第一组读长序列的第i个读长转化的矩阵中的第j个碱基C2i表示第二组读长序列的第i条读长转化的矩阵c2ij表示第二组读长序列的第i条读长转化的矩阵中的第j个碱基T1表示第一组读长序列转化为的矩阵序列的集合P1表示对第一组读长序列进行逐段延伸拼接后的矩阵序列的集合T2表示第二组读长序列转化为的矩阵序列的集合P2表示对第二组读长序列进行对位后的矩阵序列的集合N表示每一组读长序列的总个数L表示一条读长序列的长度k表示两个读长序列拼接时每一个读长序列中未重叠的碱基个数4、 模型假设1、 假设读长序列中的碱基的方向与测序方向一致，即read1中碱基方向为从左到右，read2中方向为右到左。2、 假设在拼接与组装过程中，碱基不会发生变异。3、 假设只有在检索不到完全符合符合条件的序列时才考虑测序出的误差。4、 假设两序列的相似度大于0.9则可以拼接。5、 假设k最大值等于6。5、 模型的建立与求解5.1模型建立前的准备 为了便于比较两个序列之间的关系，本文首先利用Matlab将所得的读长数据分别读入成为矩阵，并且将其中的碱基分别转化为其对应的ascii码，即转化为数值型矩阵。在拼接过程中，本文只对第一组读长数据进行拼接。5.2问题一的算法与模型5.2.1逐段延伸算法1.算法步骤：（1） 建立两个集合P1与T1，其中P1表示种子序列所在集合（初始时为空集），T1为读长序列所在集合（初始时为第一组读长序列的全体）；（2） 从T1中任意取出一条序列C1i作为种子序列放入P1中；（3） 在T1剩下的序列中检索前缀序列与C1i末端的L-k个碱基相同的序列，若检索到，转入下一步；若未检索到，则增大k的值后重复该步；若增大k值到最大后仍未找到，则还原k的值，检索与u最相似的一条序列，转下一步；（4） ；将检索到的序列末端的k个碱基拼接到C1i的末端成为长度为L=L+k的新的种子序列C1i，并将找到的序列从T集合中删除；（5） 判断u的长度是否达到要求，若未达到，重复第（3）步；若达到要求，则转第（2）步；直到T集合为空集时，停止。2.模型建立：（1） 定义相似度与的差距系数 相似度s=1-d（2） 比对过程任取T1放入P1中用与比较，如果则如果则求s=1-d，若则可以拼接，否则再取=按上述拼接。如果则求，若则可以拼接，否则继续查找重复上述过程。若T1中S都不大于0.9或者则更换种子序列，任取T1，重复上述过程直到T1=或者终止。5.2.2分离对位算法1.算法步骤：(1) 将P1中的种子序列利用生物学知识转化为其对应的反向互补序列，并且放入P2中；(2) 将P2集合中的每条序列逐项分离成与读长序列的长度L相等的m（m=u-L+1）条序列，并记下这m条序列末端距离u末端的长度，即为这m条序列的位置；(3) 在T2中检索出与m序列相同的序列，若找到，并放入对应的位置，并且将重叠部分拼接；若未找到，则检索相似度最大的序列，放入该位置；2.模型建立：组装过程，将p中所有转为反向互补序列存入中，即A（65）T（84），C（67）G（71）记为将分离并记的位置为（i,j），在中取，如果，k=1,2,3,则令若，则计算s=1-d，若s0.9则的位置为，再取，若同上操作，直至为空集综上可速，将位置中的第以分量相同的拼接即为。5.3问题二的算法与模型基于问题一的算法和模型，首先将read1和read2中的数据转化为对应的ascii码的数值型矩阵，即然后，利用逐段延伸算法和模型，read1的读长序列进行拼接，得到长度为500bp的序列 条，即再由分离对位算法将500bp的长序列的反向互补序列求出，并将每条500bp的长序列分离成88bp的序列413条，在read2中检索并放在对应的位置，形成长为500的片段，即在将这些片段拼接成为一条完整的基因，即该拼接得到的基因中含有碱基对117786对，相识度误差为6、 模型评价6.1模型的优点由于模型避免了用图像储存数据，所以占用内存空间较小，并且拼接与组装速度较快，具有一定的有效性。本文选用的基因拼接和组装方法既能快速拼接基因长读序列和有效的组装基因长读序列，又能有效的处理测序中出现的错误与重复序列。6.2模型中的不足虽然本文已引入相似度、先到先得原理对拼接和组装过程中可能出现的测序错误进行处理，但是效果仍然不是很好。6.2改进方法在拼接和组装之前应该先对读长序列结果进行筛选，将无效结果丢弃。但是出于筛选条件难以控制，太弱会保留无效序列，太强则会丢弃有效序列。因此本文没有这么做。7、 参考文献1 韩明,王家宝,李林，数学实验，同济大学出版，2012年2 司守奎,孙玺菁，数学建模算法与应用，北京：国防工业出版社，2012年3 沃森，杨焕明，基因的分子生物学， 科学出版社，2009年4 杨清，余丽芸 ，分子生物学与基因工程实验技术，中国农业大学，2014年5 李海英，杨峰山，邵淑丽，现代分子生物学与基因工程，化学工业出版社，2008年 6严海英 ，基因工程与分子生物学实验教程，武汉大学出版社，2009年8、 附录 相关程序：1.%a=fopen(C:UsersAdministratorDesktop2014BMcMc_BAC_1.fq.gz.clean.dup.clean);A,count = fscanf(a,%s); b1=find(A=+);cc=ATTGGGAAACAAGACTTGCAGAAGAACCGAAATACAATTCAAGGAAGAAGCGTGCTGAATTTTCTAAAGAGGCAACCCCGAAAGTCGC;n1=size(cc,2);for i=1:size(b1,2)endC1=;for i=1:size(b1,2)C1=C1;A(b1(i)-n1):(b1(i)-1);end a=fopen(C:UsersAdministratorDesktop2014BMcMc_BAC_2.fq.gz.clean.dup.clean);B,count = fscanf(a,%s); b2=find(B=+);cc=ATTGGGAAACAAGACTTGCAGAAGAACCGAAATACAATTCAAGGAAGAAGCGTGCTGAATTTTCTAAAGAGGCAACCCCGAAAGTCGC;n1=size(cc,2);for i=1:size(b2,2)endC2=;for i=1:size(b2,2)C2=C2;B(b2(i)-n1):(b2(i)-1);end%JY JY=;for i=1:46847JY(:,:,i)=C1(i,:);C2(i,:);end%AA AA(:,:,1)=C1(1,:);C2(1,:); AA(:,:,2)=C1(2,:);C2(2,:); for i=3:46847 AA(:,:,i)=C1(i,:);C2(i,:);end2.function varargout = wind_rose(D,F,varargin)%WIND_ROSE Wind rose of direction and intensity% % Syntax:% HANDLES = WIND_ROSE(D,I,VARARGIN)% Inputs:% D Directions% I Intensities% VARARGIN:% -dtype, type of input directions D, standard or meteo,% if meteo, the conversion dnew=mod(-90-D,360) is done;% if not meteo, standard is used (default)% -n, number of D subdivisons% -di, intensities subdivisons, default is automatic% -ci, percentage circles to draw, default is automatic% -labtitle, main title% -lablegend, legend title% -cmap, colormap jet% -colors, to use instead of colormap, for each di% -quad, Quadrant to show percentages 1% -ri, empty internal radius, relative to size of higher% percentage 1/30% -legtype, legend type: 1, continuous, 2, separated boxes 2% -bcolor, full rectangle border color none% -lcolor, line colors for axes and circles k% -percbg, percentage labels bg w% -ax, to place wind rose on pervious axes, the input for ax% must be theax x y width, where theax is the previous% axes, x and y are the location and width is the wind% rose width relative to theax width (default=1/5)% -parent, by default a new axes is created unless parent is% given, ex, parent may be a subplot% -iflip, flip the intensities as they go outward radially, ie,% highest values are placed nearest the origin 0 1% -inorm, normalize intensities, means all angles will have 100% -incout, if 0, data outside di limits will not be used 0 1% Output:% HANDLES Handles of all lines, fills, texts% Examle:% d=0:10:350;% D=;% V=;% for i=1:length(d)% n=d(i)/10;% D=D ones(1,n)*d(i);% V=V 1:n;% end% figure% wind_rose(D,V)% figure% wind_rose(D,V,iflip,1)% figure% wind_rose(D,V,ci,1 2 7,dtype,meteo)% % place it on a previous axes:% ax=axes;% plot(lon,lat)% wind_rose(D,V,ax,ax x y 1/3)% MMA 26-11-2007, % IEO, Instituto Espaol de Oceanografa% La Corua, Espaa % 10-12-2007 - Added varargin ci and n (nAngles removed as input)% 17-12-2007 - Added varargin ax, colors% 22-02-2008 - Added varargin dtype% 08-05-2008 - Bug fix (bar at dir=0 could be incorrect in some cases)% 14-05-2008 - Added varargin iflip% 16-06-2008 - Added varargin parent% 10-06-2009 - Added varargin incout handles=; % varargin options:dtype=standard;nAngles=16;ri=1/30;quad=1;legType=2;percBg=w;titStr=;legStr=;cmap=jet;colors=;Ag=; % intensity subdivs.ci=; % percentage circleslineColors=k;borderColor=none;onAxes=false;iflip=0;inorm=0;parent=0;IncHiLow=1; % include values higher and lower that the limits of Ag. vin=varargin;for i=1:length(vin) if isequal(vini,dtype) dtype=vini+1; elseif isequal(vini,n) nAngles=vini+1; elseif isequal(vini,ri) ri=vini+1; elseif isequal(vini,quad) quad=vini+1; elseif isequal(vini,legtype) legType=vini+1; elseif isequal(vini,percbg) percBg=vini+1; elseif isequal(vini,labtitle) titStr=vini+1; elseif isequal(vini,lablegend) legStr=vini+1; elseif isequal(vini,cmap) cmap=vini+1; elseif isequal(vini,colors) colors=vini+1; elseif isequal(vini,di) Ag=vini+1; elseif isequal(vini,ci) ci=vini+1; elseif isequal(vini,lcolor) lineColors=vini+1; elseif isequal(vini,bcolor) borderColor=vini+1; elseif isequal(vini,ax) ax=vini+1; try onAxes=ax(1); onAxesX=ax(2); onAxesY=ax(3); onAxesR=ax(4); catch disp(: cannot place wind rose on axes, bad argument for ax) return end elseif isequal(vini,iflip) iflip=vini+1; elseif isequal(vini,inorm) inorm=vini+1; elseif isequal(vini,parent) parent=vini+1; elseif isequal(vini,incout) IncHiLow=vini+1; endend % other options:% size of the full rectangle:rs=1.2;rl=1.7; % directions conversion:if isequal(dtype,meteo) D=mod(-90-D,360);end % angles subdivisons:D=mod(D,360);Ay=linspace(0,360,nAngles+1)-0.5*360/nAngles; % calc instensity subdivisions:if isempty(Ag) % gen Ag: f=figure(visible,off); plot(F); axis tight; yl=get(gca,ytick); close(f) dyl=diff(yl); dyl=dyl(1); if min(F)yl(end), yl=yl yl(end)+dyl; end Ag=yl;end for i=1:length(Ay)-1 if i=1 I=find( (D=Ay(i) & D=Ay(end); else I=find(D=Ay(i) & D=Ag(j) & b=Ag(j) & bAg(j+1); end E(i,j)=length(J); end if IncHiLow E(i,1)=length(find(bAg(end-1); endendb=sum(E,2)/length(D)*100; % normalize data:if inorm n=sum(E,2); for i=1:length(n) E(i,:)=E(i,:)/n(i); end b=100*ones(size(b);end % check if has values higher or lower than the Ag limitshasH=length(find(FAg(end);hasL=length(find(FAg(1); % calc number of percentage circles to draw:if isempty(ci) if inorm ci=25 50 75; g=120; ncircles=3; else dcircles=1 2 5 10 15 20 25 30 50; ncircles=3; d=abs(1./(dcircles/max(b)-ncircles); i=find(d=min(d); d=dcircles(i(1); if d*ncircles0, handles(end+1)=fill(x,y,corjcor); end r1=r2; end endendaxis equalaxis off % uistack has problems in some matlab versions, so:%uistack(labs,top)%uistack(circles,top)ch=get(wrAx,children);if inorm % only bring circles up in inorm case. for i=1:length(circles) ch(ch=circles(i)=; ch=circles(i); ch; endendfor i=1:length(labs) ch(ch=labs(i)=; ch=labs(i); ch;endset(wrAx,children,ch); % N S E W labels:bg=none;args=BackgroundColor,bg,FontSize,8;h(1)=text(-g-ri, 0,WEST, VerticalAlignment,top, HorizontalAlignment,left, args:);h(2)=text( g+ri, 0,EAST, VerticalAlignment,top, HorizontalAlignment,right,args:);h(3)=text( 0,-g-ri,SOUTH,VerticalAlignment,bottom,HorizontalAlignment,left, args:);h(4)=text( 0, g+ri,NORTH,VerticalAlignment,top, HorizontalAlignment,left, args:);handles=handles h; % scale legend:L=(g*rl-g-ri)/7;h=(g+ri)/10;dy=h/3; x0=g+ri+(g*rl-g-ri)/7;x1=x0+L;y0=-g-ri; if legType=1 % contimuous. for j=1:length(Ag)-1 lab=num2str(Ag(j); if j=1 & hasL & IncHiLow lab=; end y1=y0+h; handles(end+1)=fill(x0 x1 x1 x0,y0 y0 y1 y1,corj); handles(end+1)=text(x1+L/4,y0,lab,VerticalAlignment,middle,fontsize,8); y0=y1; end if (hasH & IncHiLow) handles(end+1)=text(x1+L/4,y0,num2str(Ag(end),VerticalAlignment,middle,fontsize,8); endelseif legType=2 % separated boxes. for j=1:length(Ag)-1 lab=num2str(Ag(j) - num2str(Ag(j+1); if j=1 & hasL & IncHiLow lab=,num2str(Ag(j); end y1=y0+h; handles(end+1)=fill(x0 x1 x1 x0,y0+dy y0+dy y1 y1,corj); handles(end+1)=text(x1+L/4,(y0+dy+y1)/2,lab,VerticalAlignment,middle,fontsize,8); y0=y1; end end % title and legend label:x=mean(-g*rs,g*rl);y=mean(g+ri,g*rs);handles(end+1)=text(x,y,titStr,HorizontalAlignment,center); x=x0;y=y1+dy;handles(end+1)=text(x,y,legStr,HorizontalAlignment,left,VerticalAlignment,bottom); if onAxes place_wr(onAxes,wrAx,onAxesX,onAxesY,onAxesR);end if nargout=1 varargout1=handles;end function place_wr(ax,ax2,x,y,width)if nargin 5 width=1/5;enduax=get(ax,units);pax=get(ax,position);set(ax,