augix.com | 里昂印象: 用python编写序列联配程序

April 02, 2005

用python编写序列联配程序

ACGGTGCGGACCAAATCCTTAGACT
ACAGTGCAAACCAAATCACACCCCTTATT

对上面两条序列的进行联配(alignement)，相同碱基配上加1分，不同碱基配对记作-1，一个gap记作-1，自左向右逐个完成配对，得到最大得分即为最佳联配。用python编写程序完成。老师写的程序只得到了一个匹配结果：

ACGGTGCGGACCAAAT______CCTTAGACT
ACAGTGCAAACCAAATCACACCCCTT__ATT

事实上，不同联配方式，得分可能均为最大值。我的程序比老师的短，并且考虑到了这一因素。我的结果是：

score=7
ACGGTGCGGACCAAAT____CC__TTAGACT
ACAGTGCAAACCAAATCACACCCCTTAT__T
//
ACGGTGCGGACCAAAT____CC__TTAGACT
ACAGTGCAAACCAAATCACACCCCTTA_T_T
//
ACGGTGCGGACCAAAT____CC__TTAGACT
ACAGTGCAAACCAAATCACACCCCTTA__TT
//
ACGGTGCGGACCAAAT____CC__TTAGACT
ACAGTGCAAACCAAATCACACCCCTT__ATT
......

我的代码是：

#!/usr/bin/python
seq1='1'+'ACGGTGCGGACCAAATCCTTAGACT'
seq2='2'+'ACAGTGCAAACCAAATCACACCCCTTATT'
flag={};score={};flag[0,0]='';score[0,0]=0;newseq1={};newseq2={};position={}
def fonc(i,j):
	flag[i,j]='';newseq1[i,j]=[];newseq2[i,j]=[]
	if j==0: flag[i,j],score[i,j]='-',score[i-1,j]-1;return
	if i==0: flag[i,j],score[i,j]='|',score[i,j-1]-1;return
	if seq1[i]==seq2[j]: s=1
	else: s=-1
	AB,AG,GB=score[i-1,j-1]+s,score[i-1,j]-1,score[i,j-1]-1
	if AB>=AG and AB>=GB: flag[i,j]+='\\'; score[i,j]=AB
	if GB>=AB and GB>=AG: flag[i,j]+='|'; score[i,j]=GB
	if AG>=AB and AG>=GB: flag[i,j]+='-'; score[i,j]=AG
for j in range(len(seq2)):
	for i in range(len(seq1)):
		if i!=0 or j!=0: fonc(i,j)
print 'score='+str(score[i,j])
newseq1[i,j]=[''];newseq2[i,j]=[''];position[i,j]='true';newseq1[0,0]=[];newseq2[0,0]=[]
while i!=0 or j!=0:
	for i,j in position.keys():
		if i==0 and j==0: break
		del position[i,j]
		if '\\' in flag[i,j]:			
			newseq1[i-1,j-1]+=[seq1[i]+x for x in newseq1[i,j]]
			newseq2[i-1,j-1]+=[seq2[j]+x for x in newseq2[i,j]]
			position[i-1,j-1]='true'
		if '-' in flag[i,j]:
			newseq1[i-1,j]+=[seq1[i]+x for x in newseq1[i,j]]
			newseq2[i-1,j]+=['_'+x for x in newseq2[i,j]]
			position[i-1,j]='true'
		if '|' in flag[i,j]:
			newseq1[i,j-1]+=['_'+x for x in newseq1[i,j]]
			newseq2[i,j-1]+=[seq2[j]+x for x in newseq2[i,j]]
			position[i,j-1]='true'
for i in range(len(newseq1[0,0])):
	print newseq1[0,0][i]+'\n'+newseq2[0,0][i]+'\n//'

Posted by 西瓜 at April 2, 2005 01:21 PM

Trackback Pings

TrackBack URL for this entry:
http://www.augix.com/cgi-bin/MT/mt-tb.cgi/69

Listed below are links to weblogs that reference 用python编写序列联配程序:

» phentermine from phentermine
[Read More]

Tracked on March 7, 2006 01:59 PM

» home mortgage from
home mortgage [Read More]

Tracked on March 17, 2006 07:36 PM

» buy ephedra from buy ephedra
[Read More]

Tracked on March 19, 2006 12:45 AM

» buy phentermine from buy phentermine
buy phentermine [Read More]

Tracked on March 19, 2006 07:59 AM

» legal ephedra from legal ephedra
[Read More]

Tracked on March 19, 2006 11:36 PM

» mortgage from mortgage
[Read More]

Tracked on March 21, 2006 12:16 PM

» mortgage from mortgage
[Read More]

Tracked on March 21, 2006 04:44 PM

» mortgage from mortgage
[Read More]

Tracked on March 21, 2006 04:56 PM

» mortgage from mortgage
[Read More]

Tracked on March 21, 2006 05:06 PM

» http://xenical.avelavok.com/detox-hgh-phentermine-quit-smoking-xenical.htm / detox hgh phentermine quit smoking xenical from
[Read More]

Tracked on April 21, 2006 12:58 AM

» http://xenical.avelavok.com/detox-hgh-phentermine-quit-smoking-xenical.htm / detox hgh phentermine quit smoking xenical from
[Read More]