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 return genomicSeq, currentAcc, currentDon
171 # build reverse complement if on negative strand
172 if strand == '-':
173 fn = 'chr%d.dna.flat' % chr
174 filename = os.path.join(dna_flat_files,fn)
175 cmd = 'wc -c %s | cut -f1 -d \' \'' % filename
176 #print cmd
177 import subprocess
178 obj = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
179 out,err = obj.communicate()
181 if err != '':
182 print 'Error occurred while trying to obtain file size'
183 end = int(out)
185 #print 'size is %d' % end
187 intervalBegin = genomicSeq_start-100
188 intervalEnd = genomicSeq_stop+100
190 #print 'before getSpliceScores'
191 #print intervalBegin ,intervalEnd
193 total_size = end
194 currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
196 currentAcc = currentAcc[100:-98]
197 currentAcc = currentAcc[1:]
198 currentDon = currentDon[100:-100]
200 length = len(genomicSeq)
201 currentAcc = currentAcc[:length]
203 currentDon = currentDon+[-inf]*(length-len(currentDon))
205 ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and genomicSeq[p-1]=='a' and genomicSeq[p]=='g']
206 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')]
208 assert ag_tuple_pos == [p for p,e in enumerate(currentAcc) if e != -inf and p > 1], pdb.set_trace()
209 assert gt_tuple_pos == [p for p,e in enumerate(currentDon) if e != -inf and p > 0], pdb.set_trace()
210 assert len(currentAcc) == len(currentDon)
212 return genomicSeq, currentAcc, currentDon