+ fixed other bug (missing n last entries of the acceptor scores)
[qpalma.git] / qpalma / sequence_utils.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 import os
5 import pdb
6 import random
7 import re
8 import sys
9
10 from numpy.matlib import inf
11
12 from Genefinding import *
13 from genome_utils import load_genomic
14
15
16 def reverse_complement(seq):
17 """
18 This function takes a read in plain or bracket notation and returns the
19 reverse complement of it.
20 I.e.
21
22 est = cgt[ac][tg]a
23
24 leads to
25
26 rev_comp = t[ac][tg]acg
27 """
28
29 bpos = seq.find('[')
30 rc = lambda x: {'a':'t','c':'g','g':'c','t':'a'}[x]
31
32 # check first whether seq contains no brackets at all
33 if bpos == -1:
34 ret_val = map(rc,seq)
35 ret_val.reverse()
36 ret_val = "".join(ret_val)
37 else:
38 brc = lambda x: {'a':'t','c':'g','g':'c','t':'a','[':'[',']':']'}[x]
39
40 # first_part is the part of the seq up to the first occurrence of a
41 # bracket
42 first_part = seq[:bpos]
43 first_part = map(rc,first_part)
44 first_part.reverse()
45 first_part = "".join(first_part)
46
47 # inside brackets has to be complemented but NOT reversed
48 inside_brackets = seq[bpos+1:bpos+3]
49 inside_brackets = "".join(map(rc,inside_brackets))
50
51 ret_val = '%s[%s]%s'%(reverse_complement(seq[bpos+4:]),inside_brackets,first_part)
52
53 return ret_val
54
55
56 def unbracket_est(est):
57 """
58 This function takes a read in bracket notation and restores the read sequence from it.
59 I.e.
60
61 est = cgt[ac][tg]aa
62
63 leads to
64
65 result = cgtcgaa
66
67 so the second entry within brackets is the base on the read whereas the first
68 entry is the base from the dna.
69 """
70
71 new_est = ''
72 e = 0
73
74 while True:
75 if e >= len(est):
76 break
77
78 if est[e] == '[':
79 new_est += est[e+2]
80 e += 4
81 else:
82 new_est += est[e]
83 e += 1
84
85 return "".join(new_est).lower()
86
87
88 def create_bracket_seq(dna_seq,read_seq):
89 """
90 This function takes a dna sequence and a read sequence and returns the
91 bracket format of the match/mismatches i.e.
92
93 dna : aaa
94 read: aac
95 is written in bracket notation: aa[ac]
96 """
97 assert len(dna_seq) == len(read_seq)
98 return "".join(map(lambda x,y: ['[%s%s]'%(x,y),x][x==y],dna_seq,read_seq))
99
100
101 def getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size=0):
102 """
103 Now we want to use interval_query to get the predicted splice scores trained
104 on the TAIR sequence and annotation.
105 """
106
107 size = intervalEnd-intervalBegin
108 assert size > 1, 'Error (getSpliceScores): interval size is less than 2!'
109
110 acc = size*[0.0]
111 don = size*[0.0]
112
113 interval_matrix = createIntArrayFromList([intervalBegin,intervalEnd])
114 pos_size = new_intp()
115 intp_assign(pos_size,1)
116
117 # fetch acceptor scores
118 sscore_filename = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/acc/contig_%d%s'
119 acc = doIntervalQuery(chr,strand,intervalBegin,intervalEnd,sscore_filename,total_size)
120
121 # fetch donor scores
122 sscore_filename = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/don/contig_%d%s'
123 don = doIntervalQuery(chr,strand,intervalBegin,intervalEnd,sscore_filename,total_size)
124
125 return acc, don
126
127
128 def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_files):
129 """
130 This function expects an interval, chromosome and strand information and
131 returns then the genomic sequence of this interval and the associated scores.
132 """
133
134 assert genomicSeq_start < genomicSeq_stop
135
136 chrom = 'chr%d' % chr
137 genomicSeq = load_genomic(chrom,strand,genomicSeq_start-1,genomicSeq_stop,dna_flat_files,one_based=False)
138 genomicSeq = genomicSeq.lower()
139
140 # check the obtained dna sequence
141 assert genomicSeq != '', 'load_genomic returned empty sequence!'
142
143 # all entries other than a c g t are set to n
144 no_base = re.compile('[^acgt]')
145 genomicSeq = no_base.sub('n',genomicSeq)
146
147 if strand == '+':
148 intervalBegin = genomicSeq_start-100
149 intervalEnd = genomicSeq_stop+100
150
151 currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd)
152
153 currentAcc = currentAcc[100:-98]
154 currentAcc = currentAcc[1:]
155 currentDon = currentDon[100:-100]
156
157 length = len(genomicSeq)
158 currentAcc = currentAcc[:length]
159
160 currentDon = currentDon+[-inf]*(length-len(currentDon))
161
162 ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and genomicSeq[p-1]=='a' and genomicSeq[p]=='g' ]
163 gt_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>0 and p<len(genomicSeq)-1 and e=='g' and (genomicSeq[p+1]=='t' or genomicSeq[p+1]=='c')]
164
165 assert ag_tuple_pos == [p for p,e in enumerate(currentAcc) if e != -inf and p > 1], pdb.set_trace()
166 assert gt_tuple_pos == [p for p,e in enumerate(currentDon) if e != -inf and p > 0], pdb.set_trace()
167 assert len(currentAcc) == len(currentDon)
168
169 pdb.set_trace()
170
171 return genomicSeq, currentAcc, currentDon
172
173 # build reverse complement if on negative strand
174 if strand == '-':
175 fn = 'chr%d.dna.flat' % chr
176 filename = os.path.join(dna_flat_files,fn)
177 cmd = 'wc -c %s | cut -f1 -d \' \'' % filename
178 import subprocess
179 obj = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
180 out,err = obj.communicate()
181
182 if err != '':
183 print 'Error occurred while trying to obtain file size'
184 end = int(out)
185
186 intervalBegin = genomicSeq_start#-100
187 intervalEnd = genomicSeq_stop#+100
188
189 total_size = end
190 currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
191
192 #currentAcc = currentAcc[100:-98]
193 currentAcc = currentAcc[1:]
194 #currentDon = currentDon[100:-100]
195
196 length = len(genomicSeq)
197 currentAcc = currentAcc[:length]
198 currentAcc = currentAcc+[-inf]*(length-len(currentAcc))
199
200 currentDon = currentDon+[-inf]*(length-len(currentDon))
201
202 ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and genomicSeq[p-1]=='a' and genomicSeq[p]=='g']
203 gt_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>0 and p<len(genomicSeq)-1 and e=='g' and (genomicSeq[p+1]=='t' or genomicSeq[p+1]=='c')]
204
205 assert ag_tuple_pos == [p for p,e in enumerate(currentAcc) if e != -inf and p > 1], pdb.set_trace()
206 assert gt_tuple_pos == [p for p,e in enumerate(currentDon) if e != -inf and p > 0], pdb.set_trace()
207 assert len(currentAcc) == len(currentDon), pdb.set_trace()
208
209 return genomicSeq, currentAcc, currentDon