1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
4 import os
5 import pdb
6 import random
7 import re
8 import sys
10 from numpy.matlib import inf
12 from Genefinding import *
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.
22 est = cgt[ac][tg]a
26 rev_comp = t[ac][tg]acg
27 """
29 bpos = seq.find('[')
30 rc = lambda x: {'a':'t','c':'g','g':'c','t':'a'}[x]
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]
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)
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))
51 ret_val = '%s[%s]%s'%(reverse_complement(seq[bpos+4:]),inside_brackets,first_part)
53 return ret_val
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.
61 est = cgt[ac][tg]aa
65 result = cgtcgaa
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 """
71 new_est = ''
72 e = 0
74 while True:
75 if e >= len(est):
76 break
78 if est[e] == '[':
79 new_est += est[e+2]
80 e += 4
81 else:
82 new_est += est[e]
83 e += 1
85 return "".join(new_est).lower()
89 """
90 This function takes a dna sequence and a read sequence and returns the
91 bracket format of the match/mismatches i.e.
93 dna : aaa
95 is written in bracket notation: aa[ac]
96 """
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 """
107 size = intervalEnd-intervalBegin
108 assert size > 1, 'Error (getSpliceScores): interval size is less than 2!'
110 acc = size*[0.0]
111 don = size*[0.0]
113 interval_matrix = createIntArrayFromList([intervalBegin,intervalEnd])
114 pos_size = new_intp()
115 intp_assign(pos_size,1)
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)
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)
125 return acc, don
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 """
134 assert genomicSeq_start < genomicSeq_stop
136 chrom = 'chr%d' % chr
138 genomicSeq = genomicSeq.lower()
140 # check the obtained dna sequence
141 assert genomicSeq != '', 'load_genomic returned empty sequence!'
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)
147 if strand == '+':
148 intervalBegin = genomicSeq_start-100
149 intervalEnd = genomicSeq_stop+100
151 currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd)
153 currentAcc = currentAcc[100:-98]
154 currentAcc = currentAcc[1:]
155 currentDon = currentDon[100:-100]
157 length = len(genomicSeq)
158 currentAcc = currentAcc[:length]
160 currentDon = currentDon+[-inf]*(length-len(currentDon))
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')]
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)
169 pdb.set_trace()
171 return genomicSeq, currentAcc, currentDon
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()
182 if err != '':
183 print 'Error occurred while trying to obtain file size'
184 end = int(out)
186 intervalBegin = genomicSeq_start#-100
187 intervalEnd = genomicSeq_stop#+100
189 total_size = end
190 currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
192 #currentAcc = currentAcc[100:-98]
193 currentAcc = currentAcc[1:]
194 #currentDon = currentDon[100:-100]
196 length = len(genomicSeq)
197 currentAcc = currentAcc[:length]
198 currentAcc = currentAcc+[-inf]*(length-len(currentAcc))
200 currentDon = currentDon+[-inf]*(length-len(currentDon))
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')]
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()
209 return genomicSeq, currentAcc, currentDon