vcNtd2G6

1#This software is a free software. 
2#Thus, it is licensed under GNU General Public License.
3#Python implementation to Nussinov Folding Algorithm 
4#for Homework 7 of Bioinformatics class.
5#Forrest Bao, Nov. 29 <http://fsbao.net> <forrest.bao aT gmail.com>
6
7import sys,string
8from numpy import *
9from matplotlib import *
10
11#read sequences
12#sequences are stored in many lines
13f=open(sys.argv[1], 'r')
14seq=[]
15for line in f.readlines():
16    seq.append(line.strip('\n'));
17
18def delta(l,m):
19    delta=0;
20    if l=='A' and (m=='U' or m=='T'):
21        return 2;
22    elif (l=='U' or l=='T') and m=='A':
23        return 2;
24    elif l=='G' and m=='C':
25        return 3;
26    elif l=='C' and m=='G':
27        return 3;
28    else:
29        return 0;
30
31def buildDP(seq):
32    L=len(seq);
33    s=zeros((L,L));
34    for n in range(1,L):
35        for j in range(n,L):
36            i=j-n;
37            if abs(i-j)>3:
38                case1=s[i+1,j-1]+delta(seq[i],seq[j]);
39                case2=s[i+1,j];
40                case3=s[i,j-1];
41                if i+3<=j:
42                    tmp=[];
43                    for k in range(i+1,j):
44                        tmp.append(s[i,k]+s[k+1,j]);
45                    case4=max(tmp);
46                    s[i,j]=max(case1,case2,case3,case4);
47                else:
48                    s[i,j]=max(case1,case2,case3);
49            else:        
50                case2=s[i+1,j];
51                case3=s[i,j-1];
52                if i+3<=j:
53                    tmp=[];
54                    for k in range(i+1,j):
55                        tmp.append(s[i,k]+s[k+1,j]);
56                    case4=max(tmp);
57                    s[i,j]=max(case2,case3,case4);
58                else:
59                    s[i,j]=max(case2,case3);
60    print("Score: ", s[0,L-1])                
61    return s;
62
63def traceback(s,seq,i,j,pair):
64    if i<j:
65        if s[i,j]==s[i+1,j]:
66            traceback(s,seq,i+1,j,pair);
67        elif s[i,j]==s[i,j-1]:
68            traceback(s,seq,i,j-1,pair);
69        elif s[i,j]==s[i+1,j-1]+delta(seq[i],seq[j]):
70            pair.append([i,j,str(seq[i]),str(seq[j])]);
71            traceback(s,seq,i+1,j-1,pair);
72        else:
73            for k in range(i+1,j):
74                if s[i,j]==s[i,k]+s[k+1,j]:
75                    traceback(s,seq,i,k,pair);
76                    traceback(s,seq,k+1,j,pair);
77                    break;
78    return pair;
79
80for q in range(0,len(seq)):
81    pair=traceback(buildDP(seq[q]),seq[q],0,len(seq[q])-1,[])
82    print("max # of folding pairs: ",len(pair));
83    for x in range(0,len(pair)):
84        print('%d %d %s==%s' % (pair[x][0],pair[x][1],pair[x][2],pair[x][3]));
85    print("---");