+ fixed a bug in the C++ interface
[qpalma.git] / standalone / palma / genomic.py
1 #
2 # This program is free software; you can redistribute it and/or modify
3 # it under the terms of the GNU General Public License as published by
4 # the Free Software Foundation; either version 2 of the License, or
5 # (at your option) any later version.
6 #
7 # Written (W) 2006-2007 Soeren Sonnenburg
8 # Written (W) 2006-2007 Mikio Braun
9 # Copyright (C) 2007 Fraunhofer Institute FIRST and Max-Planck-Society
10 #
11
12 import time
13 from string import maketrans
14
15 """ this function is 100% compatible to the matlab function, thus it is one based (!)
16 use one_based=False if needed, then however the interval is [start,stop) (excluding stop)
17 """
18 def load_genomic(chromosome, strand, start, stop, genome, one_based=True):
19 fname = '/fml/ag-raetsch/share/databases/genomes/' + genome + '/' + chromosome[3:] + '.flat'
20 f=file(fname)
21 if one_based:
22 f.seek(start-1)
23 str=f.read(stop-start+1)
24 else:
25 f.seek(start)
26 str=f.read(stop-start)
27
28 if strand=='-':
29 return reverse_complement(str)
30 elif strand=='+':
31 return str
32 else:
33 print 'strand must be + or -'
34 raise KeyError
35
36 """ read a table browser ascii output file (http://genome.ucsc.edu/cgi-bin/hgTables) """
37 def read_table_browser(f):
38 table=dict();
39 for l in f.readlines():
40 if not l.startswith('#'):
41 (name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds,proteinID,alignID)=l.split('\t')
42 exonStarts=[ int(i) for i in exonStarts.split(',')[:-1] ]
43 exonEnds=[ int(i) for i in exonEnds.split(',')[:-1] ]
44
45 table[name]={ 'chrom': chrom, 'strand': strand, 'txStart': int(txStart), 'txEnd': int(txEnd),
46 'cdsStart': int(cdsStart), 'cdsEnd': int(cdsEnd), 'exonCount': int(exonCount), 'exonStarts': exonStarts,
47 'exonEnds': exonEnds, 'proteinID': proteinID, 'alignID': alignID[:-1] }
48
49 return table
50
51 """ get promoter region """
52 def get_promoter_region(chromosome, strand, gene_start, gene_end, genome, length):
53
54 if strand == '+':
55 return load_genomic(chromosome, strand, gene_start, gene_start+length, genome, one_based=False)
56 elif strand == '-':
57 return load_genomic(chromosome, strand, gene_end, gene_end+length, genome, one_based=False)
58 else:
59 print 'unknown strand'
60 return None
61
62 """ reverse + complement a DNA sequence (only letters ACGT are translated!)
63 FIXME won't work with all the rest like y... """
64 def reverse_complement(str):
65 t=maketrans('acgtACGT','tgcaTGCA')
66 return str[len(str)::-1].translate(t)
67
68 """ works only with .fa files that contain a single entry """
69 def read_single_fasta(fname):
70 str=file(fname).read()
71 str=str[str.index('\n')+1:].replace('\n','')
72 return str
73
74 """ writes only single enty .fa files """
75 def write_single_fasta(fname, name, str, linelen=60):
76 header= '>' + name + '\n'
77 f=file(fname,'a')
78 f.write(header)
79 for i in xrange(0,len(str),linelen):
80 f.write(str[i:i+linelen]+'\n')
81 f.close()
82
83 """ read fasta as dictionary """
84 def read_fasta(f, fasta=dict(), id_only=True):
85 import numpy
86 str=f.read()
87 idx = str.find('>')
88 #if idx==-1:
89 if 0: # no support for contig list files -> would need extra blastn workaround
90 # file name list?
91 files = str.split('\n') ;
92 for s in files:
93 if len(s)==0: continue
94 print s.strip() + '\n'
95 fasta = read_fasta(file(s.strip()), fasta)
96 else:
97 # real fasta file?
98 sequences = str.split('>')
99 for s in sequences:
100 if len(s)==0: continue
101 header = s[0:s.index('\n')]
102 if id_only:
103 header = header.split(' ')[0] ;
104 header = header.split('\t')[0] ;
105 seq = s[s.index('\n')+1:].replace('\n','').upper()
106 #print 'has_key', fasta.has_key(header),header
107 assert(not fasta.has_key(header))
108 fasta[header]=seq ;
109
110 return fasta
111
112 """ write dictionary fasta """
113 def write_fasta(f, d, linelen=60):
114 for k in sorted(d):
115 f.write('>%s\n' % k);
116 s = d[k]
117 for i in xrange(0, len(s), linelen):
118 f.write(s[i:i+linelen] + '\n')
119
120 def write_gff(f, (source, version), (seqtype, seqname), descrlist, skipheader=False):
121 """ writes a gff version 2 file
122 descrlist is a list of dictionaries, each of which contain these fields:
123 <seqname> <source> <feature> <start> <end> <score> <strand> <frame> [attributes] [comments]
124 """
125
126 if not skipheader:
127 f.write('##gff-version 2\n')
128 f.write('##source-version %s %s\n' % (source, version) )
129
130 t=time.localtime()
131 f.write("##date %d-%d-%d %d:%d:%d\n" % t[0:6])
132
133 f.write('##Type %s %s\n' % (seqtype, seqname) )
134
135 for d in descrlist:
136 f.write('%s\t%s\t%s\t%d\t%d\t%f\t%s\t%d' % (d['seqname'], d['source'],
137 d['feature'], d['start'], d['end'],
138 d['score'], d['strand'], d['frame']))
139 if d.has_key('attributes'):
140 f.write('\t' + d['attributes'])
141 if d.has_key('comments'):
142 f.write('\t' + d['comments'])
143 f.write('\n')
144
145
146 if __name__ == '__main__':
147 import sys,os
148
149 table=read_table_browser(file('/fml/ag-raetsch/home/sonne/addnet/tfbs/share/data/wt1_bibliosphere_table_browser_hg17.txt'))
150 print table.keys()
151 print table[table.keys()[0]]
152 d = { 'ahoernchen' : 'ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT',
153 'bhoernchen' : 'GATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACAGATTACA' }
154
155 write_fasta(sys.stdout, d)
156 write_fasta(file('/tmp/test.fa','w'), d)
157
158 d2 = read_fasta(file('/tmp/test.fa'))
159 os.unlink('/tmp/test.fa')
160
161 print d
162 print d2
163 print d == d2
164
165 p=load_genomic('chr5', '+', 100000, 100100,'hg17')
166 n=load_genomic('chr1', '-', 3000000, 3001000,'mm7')
167 write_single_fasta('bla.fa','bla', 'ACGT')
168 n2=read_single_fasta('bla.fa')