脚本之家,脚本语言编程技术及教程分享平台!
分类导航

Python|VBS|Ruby|Lua|perl|VBA|Golang|PowerShell|Erlang|autoit|Dos|bat|

服务器之家 - 脚本之家 - Python - 使用python实现BLAST

使用python实现BLAST

2021-01-15 00:30脚本之家 Python

最近在自学python,又用python实现了一下BLAST。 这次更新了打分函数如下,空位罚分改为-5,但不区分gap open 和 gap extend。

最近在自学python,又用python实现了一下BLAST。

这次更新了打分函数如下,空位罚分改为-5,但不区分gap open 和 gap extend。

使用python实现BLAST

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
'''''
@author: JiuYu
'''
 
def score(a,b):#scoring function
  score=0
  lst=['AC','GT','CA','TG']
  if a==b:
    score +=2
  elif a+b in lst:
    score += -5
  else:
    score += -7
  return score
 
def BLAST(seq1,seq2):#Basic Local Alignment Search Tool
  l1 = len(seq1)
  l2 = len(seq2)
  GAP =-5   #-5 for any gap
  scores =[]
  point =[]
   
  for j in range(l2+1):
    if j == 0:
      line1=[0]
      line2=[0]
      for i in range(1,l1+1):
        line1.append(GAP*i)
        line2.append(2)
    else:
      line1=[]
      line2=[]
      line1.append(GAP*j)
      line2.append(3)
    scores.append(line1)
    point.append(line2)
   
  #fill the blank of scores and point
  for j in range(1,l2+1):
    letter2 = seq2[j-1]
    for i in range(1,l1+1):
      letter1 = seq1[i-1]
      diagonal_score = score(letter1, letter2) + scores[j-1][i-1]
      left_score = GAP + scores[j][i-1]
      up_score = GAP + scores[j-1][i]
      max_score = max(diagonal_score, left_score, up_score)
      scores[j].append(max_score)
       
      if scores[j][i] == diagonal_score:
        point[j].append(1)
      elif scores[j][i] == left_score:
        point[j].append(2)
      else:
        point[j].append(3)
         
  #trace back
  alignment1=''
  alignment2=''
  i = l2
  j = l1
  print 'scores =',scores[i][j]
  while True:
    if point[i][j] == 0:
      break
    elif point[i][j] == 1:
      alignment1 += seq1[j-1]
      alignment2 += seq2[i-1]
      i -= 1
      j -= 1
    elif point[i][j] == 2:
      alignment1 += seq1[j-1]
      alignment2 += '-'
      j -= 1
    else:
      alignment1 += '-'
      alignment2 += seq2[i-1]
      i -= 1
       
  #reverse alignment
  alignment1 = alignment1[::-1]
  alignment2 = alignment2[::-1]
  print 'The best alignment:'
  print alignment1
  print alignment2
 
seq1=raw_input('Please input your first sequences:\n')
seq2=raw_input('input second sequences:\n')
BLAST(seq1, seq2)

运行结果:

使用python实现BLAST

无疑python对字符串的处理更加强大,语言也更加简单,优雅。比如最后逆序输出alignment,java我是单独写了一个逆序函数,而python只用一个语句就可以完成相同任务。

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持服务器之家。

延伸 · 阅读

精彩推荐