« back to all changes in this revision
Viewing changes to src/interPopula/HapMap/__init__.py
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Tiago Antao
- Date: 2010-04-07 13:26:50 UTC
- Revision ID: tiagoantao@gmail.com-20100407132650-dv1yrnw3rvbtypdo
initial import
- files added:
- epydoc
- scripts
- src
- src/interPopula
- src/interPopula/HapMap
- src/interPopula/TableBrowser
- tests
- tests/adhoc
- tests/unit
added
removed
src/interPopula/HapMap/__init__.py
1
from ftplib import FTP
2
from sys import exit
3
import gzip
4
import os
5
import sqlite3
6
from interPopula import Config
7
8
def prepareDataDir():
9
hapMapDir = Config.dataDir + os.sep + 'HapMap'
10
#print hapMapDir
11
try:
12
os.mkdir(hapMapDir)
13
except OSError:
14
pass #Probably OK, dir already exists
15
dbConn = sqlite3.connect(hapMapDir + os.sep + 'db')
16
17
try: #indiv
18
dbConn.execute('''CREATE TABLE individual (
19
fid VARCHAR(15),
20
iid VARCHAR(15),
21
dad VARCHAR(15),
22
mom VARCHAR(15),
23
gender INTEGER,
24
pop VARCHAR(15)
25
)''')
26
except sqlite3.OperationalError, e:
27
print e
28
pass #tables already exist...OK
29
30
try: #FDB
31
dbConn.execute('''CREATE TABLE freq_snp (
32
rs_id INTEGER,
33
population VARCHAR(20),
34
chromosome VARCHAR(2),
35
position INTEGER,
36
strand CHARACTER(1),
37
nuc1 CHARACTER(1),
38
nuc2 CHARACTER(1),
39
homo1 INTEGER,
40
hetero INTEGER,
41
homo2 INTEGER
42
)''')
43
dbConn.execute('CREATE INDEX fq_snp_rs ON freq_snp(rs_id)')
44
dbConn.execute('CREATE INDEX fq_snp_pop ON freq_snp(population)')
45
dbConn.execute('CREATE INDEX fq_snp_chr ON freq_snp(chromosome)')
46
except sqlite3.OperationalError, e:
47
print e
48
pass #tables already exist...OK
49
50
try:#Genome
51
dbConn.execute('''CREATE TABLE geno_snp (
52
rs_id INTEGER,
53
population VARCHAR(20),
54
chromosome VARCHAR(2),
55
iid VARCHAR(15),
56
nuc1 CHARACTER(1),
57
nuc2 CHARACTER(1)
58
)''')
59
dbConn.execute('CREATE INDEX gn_snp_rs ON geno_snp(rs_id)')
60
dbConn.execute('CREATE INDEX gn_snp_pop ON geno_snp(population)')
61
dbConn.execute('CREATE INDEX gn_snp_chr ON geno_snp(chromosome)')
62
dbConn.execute('CREATE INDEX gn_snp_cs_pop ON geno_snp(chromosome, population)')
63
except sqlite3.OperationalError, e:
64
print e
65
pass #tables already exist...OK
66
67
return dbConn
68
69
70
class Indiv:
71
"""Database of map between indivs and pops and pedigree
72
"""
73
74
def __init__(self):
75
self.dbConn = prepareDataDir()
76
77
def get(self, forceLoad = False):
78
'''Get data.
79
80
This will retrieve data from HapMap and load it on the database.
81
Load might fail, repeating might be a possibility.
82
'''
83
hapMapDir = Config.dataDir + os.sep + 'HapMap'
84
if forceLoad or (len(self.getPops()) == 0):
85
tempHMFile = hapMapDir + os.sep + 'tmp_i.txt'
86
try:
87
os.remove(tempHMFile)
88
except OSError:
89
pass #Doesn't exist -> OK
90
ftp = FTP('ftp.ncbi.nlm.nih.gov')
91
ftpDir = 'hapmap/samples_individuals'
92
ftp.login()
93
flist = ftp.nlst(ftpDir)
94
#print 'FL', flist
95
for fname in flist:
96
#XXX: if we dont fetch nothing we should halt and report
97
if fname.find("relation") > -1:
98
print 'Fetching' , fname
99
ftp.retrbinary('RETR ' + fname,
100
open(tempHMFile, 'wb').write)
101
break
102
ftp.close()
103
104
self.dbConn.execute("DELETE FROM individual")
105
f = open(tempHMFile, 'r')
106
l = f.readline() #we jump the first
107
108
l = f.readline()
109
while l<>'':
110
toks = l.rstrip().split('\t')
111
#print toks
112
fid = toks[0]
113
iid = toks[1]
114
dad = toks[2]
115
mom = toks[3]
116
gender = int(toks[4])
117
#pheno not loaded
118
pop = toks[6]
119
self.dbConn.execute('''
120
INSERT INTO individual (
121
fid, iid, dad, mom, gender, pop)
122
VALUES (?,?,?,?,?,?)''',
123
(fid, iid, dad, mom, gender, pop))
124
l = f.readline()
125
f.close()
126
self.dbConn.commit()
127
128
def cleanDB(self):
129
'''Cleans the database. Caution...
130
'''
131
self.dbConn.execute('DELETE FROM individual')
132
133
def nukeDB(self):
134
'''Removes tables from the database. Caution...
135
'''
136
self.dbConn.execute('DROP TABLE individual')
137
138
def getPops(self):
139
"""Get the list of available pops."""
140
c = self.dbConn.cursor()
141
c.execute('''
142
SELECT DISTINCT pop
143
FROM individual''')
144
p_list = []
145
for rs in c:
146
p_list.append(rs)
147
return p_list
148
149
def getIndivsForPop(self, pop):
150
'''Gets a list of individuals for a pop
151
format: fid, iid, dad, mom, gender
152
'''
153
c = self.dbConn.cursor()
154
c.execute('''
155
SELECT fid, iid, dad, mom, gender
156
FROM individual
157
WHERE pop = ?''', (pop))
158
i_list = []
159
for rs in c:
160
i_list.append(rs)
161
return i_list
162
163
def getOfsForPop(self, pop):
164
'''Gets a list of offspring for a pop
165
format: iid
166
'''
167
c = self.dbConn.cursor()
168
c.execute('''
169
SELECT iid
170
FROM individual
171
WHERE dad IS NOT NULL
172
AND pop = ?''', (pop))
173
i_list = []
174
for rs in c:
175
i_list.append(rs)
176
return i_list
177
178
179
def close(self):
180
'''Closes the database.
181
'''
182
self.dbConn.execute('VACUUM')
183
self.dbConn.close()
184
185
186
187
class FDB:
188
"""Database of frequency/genotype HapMap data.
189
Note: Local database
190
"""
191
192
def __init__(self):
193
self.dbConn = prepareDataDir()
194
195
def getPopsForChr(self, chr):
196
'''Check if database has a certain chromosome. Returns list of pops.
197
'''
198
pops = []
199
c = self.dbConn.cursor()
200
c.execute('''
201
SELECT DISTINCT population
202
FROM freq_snp
203
WHERE chromosome = ?''', chr)
204
for pop in c:
205
pops.append(pop)
206
return pops
207
208
def getChrsForPop(self, pop):
209
'''Check if database has a certain population. Returns list of chrs.
210
'''
211
chrs = []
212
c = self.dbConn.cursor()
213
c.execute('''
214
SELECT chromosome
215
FROM freq_snp
216
WHERE population = ?''', pop)
217
for chr in c:
218
chrs.append(chr)
219
return chrs
220
221
def hasChrPop(self, chr, pop):
222
'''Check if database has a certain chromosome and population.
223
'''
224
c = self.dbConn.cursor()
225
#c.execute('select * from freq_snp')
226
#for line in c:
227
# print line
228
c.execute('''
229
SELECT count(*)
230
FROM freq_snp
231
WHERE chromosome = ?
232
AND population = ?''', (chr, pop))
233
cnt = iter(c).next()
234
#print chr, pop, cnt
235
if cnt == (0,):
236
return False
237
else:
238
return True
239
240
def requireChrPop(self, chr, pop, forceLoad = False):
241
"""Require a chromosome and population
242
"""
243
self.requireChrPopFreq(chr, pop, forceLoad)
244
self.requireChrPopGeno(chr, pop, forceLoad)
245
246
def requireChrPopFreq(self, chr, pop, forceLoad = False):
247
'''Requires a chromosome and population freq.
248
249
This will retrieve data from HapMap and load it on the database.
250
Hogs network and takes time and disk space!
251
Load might fail, repeating might be a possibility.
252
'''
253
hapMapDir = Config.dataDir + os.sep + 'HapMap'
254
if forceLoad or (not self.hasChrPop(chr, pop)):
255
tempHMFile = hapMapDir + os.sep + 'tmp_hp.gz'
256
try:
257
os.remove(tempHMFile)
258
except OSError:
259
pass #Doesn't exist -> OK
260
print 'Will load', chr, pop, 'NOTE: This will take a lot of time'
261
ftp = FTP('ftp.hapmap.org')
262
ftpDir = 'hapmap/frequencies/latest/forward/non-redundant'
263
ftp.login()
264
flist = ftp.nlst(ftpDir)
265
#XXX: We are loading fwd_strand, non redundant... report on paper
266
#print 'FL', flist
267
for fname in flist:
268
#XXX: if we dont fetch nothing we should halt and report
269
if fname.find('_chr' + chr + '_' + pop) > -1:
270
print 'Fetching' , fname
271
ftp.retrbinary('RETR ' + fname,
272
open(tempHMFile, 'wb').write)
273
ftp.close()
274
gz = gzip.open(tempHMFile, 'rb')
275
l = gz.readline() #we jump the first
276
l = gz.readline()
277
print 'Injecting data'
278
self.dbConn.execute('''
279
DELETE
280
FROM freq_snp
281
WHERE chromosome = ?
282
AND population = ?''', (chr, pop))
283
while l<>'':
284
toks = l.rstrip().split(' ')
285
#print toks
286
rs_id = int(toks[0][2:])
287
pos = int(toks[2])
288
strand = toks[3]
289
nuc1 = toks[13][0]
290
nuc2 = toks[13][2]
291
homo1 = int(toks[12])
292
hetero = int(toks[15])
293
homo2 = int(toks[18])
294
self.dbConn.execute('''
295
INSERT INTO freq_snp (
296
rs_id, population, chromosome, position, strand,
297
nuc1, nuc2, homo1, hetero, homo2)
298
VALUES (?,?,?,?,?,?,?,?,?,?)''',
299
(rs_id, pop, chr, pos, strand,
300
nuc1, nuc2, homo1, hetero, homo2))
301
l = gz.readline()
302
gz.close()
303
self.dbConn.commit()
304
305
def requireChrPopGeno(self, chr, pop, forceLoad = False):
306
'''Requires a chromosome and population geno.
307
308
This will retrieve data from HapMap and load it on the database.
309
Hogs network and takes time and disk space!
310
Load might fail, repeating might be a possibility.
311
'''
312
hapMapDir = Config.dataDir + os.sep + 'HapMap'
313
if forceLoad or (not self.hasChrPopGeno(chr, pop)):
314
tempHMFile = hapMapDir + os.sep + 'tmp_gn.gz'
315
try:
316
os.remove(tempHMFile)
317
except OSError:
318
pass #Doesn't exist -> OK
319
print 'Will load geno ', chr, pop, 'NOTE: will take a lot of time'
320
ftp = FTP('ftp.hapmap.org')
321
ftpDir = 'hapmap/genotypes/latest/forward/non-redundant'
322
ftp.login()
323
flist = ftp.nlst(ftpDir)
324
#XXX: We are loading fwd_strand, non redundant... report on paper
325
#print 'FL', flist
326
for fname in flist:
327
#XXX: if we dont fetch nothing we should halt and report
328
if fname.find('_chr' + chr + '_' + pop) > -1:
329
print 'Fetching' , fname
330
ftp.retrbinary('RETR ' + fname,
331
open(tempHMFile, 'wb').write)
332
ftp.close()
333
gz = gzip.open(tempHMFile, 'rb')
334
l = gz.readline() #we get indivs
335
indivs = l.rstrip().split(' ')[11:]
336
l = gz.readline()
337
print 'Injecting data'
338
self.dbConn.execute('''
339
DELETE
340
FROM geno_snp
341
WHERE chromosome = ?
342
AND population = ?''', (chr, pop))
343
while l<>'':
344
toks = l.rstrip().split(' ')
345
#print toks
346
rs_id = int(toks[0][2:])
347
#print rs_id
348
for i in range(len(indivs)):
349
nuc1 = toks[11+i][0]
350
nuc2 = toks[11+i][1]
351
self.dbConn.execute('''
352
INSERT INTO geno_snp (
353
rs_id, population, chromosome, iid, nuc1, nuc2)
354
VALUES (?,?,?,?,?,?)''',
355
(rs_id, pop, chr, indivs[i], nuc1, nuc2))
356
l = gz.readline()
357
gz.close()
358
self.dbConn.commit()
359
360
361
362
def requirePop(self, pop):
363
'''Requires a population, all chromosomes.
364
365
This will retrieve data from HapMap and load it on the database.
366
Hogs network and takes time and disk space!
367
Load might fail, repeating might be a possibility.
368
'''
369
for chr in range(1, 23):
370
self.requireChrPop(str(chr), pop)
371
self.requireChrPop('X', pop)
372
self.requireChrPop('Y', pop)
373
374
def requireChr(self, chr, pop_list):
375
'''Requires a chromosome for specified populations.
376
377
This will retrieve data from HapMap and load it on the database.
378
379
We let the user decide which populations to load, mainly because
380
we don't know between (JPT and CHB) or JPT+CHB.
381
382
Hogs network and takes time and disk space!
383
384
Load might fail, repeating might be a possibility.
385
'''
386
for pop in pop_list:
387
self.requireChrPop(chr, pop)
388
389
def cleanDB(self):
390
'''Cleans the database. Caution...
391
'''
392
self.dbConn.execute('DELETE FROM freq_snp')
393
self.dbConn.execute('DELETE FROM geno_snp')
394
395
def nukeDB(self):
396
'''Removes tables from the database. Caution...
397
'''
398
self.dbConn.execute('DROP TABLE freq_snp')
399
self.dbConn.execute('DROP TABLE geno_snp')
400
401
def getRSsForInterval(self, chr, begin, end):
402
'''Gets a list of rs_ids for an interval (limits included)
403
'''
404
c = self.dbConn.cursor()
405
c.execute('''
406
SELECT DISTINCT rs_id
407
FROM freq_snp
408
WHERE chromosome = ?
409
AND position >= ?
410
AND position <= ?''', (chr, begin, end))
411
rs_list = []
412
for rs in c:
413
rs_list.append(rs[0])
414
return rs_list
415
416
def getSNPs(self, rs_id):
417
'''Returns a list of SNPs for an rs_id (1 per pop with that SNP).
418
419
@param rs_id: rsId.
420
@return: List of SNPs.
421
'''
422
c = self.dbConn.execute('''
423
SELECT *
424
FROM freq_snp
425
WHERE rs_id = ?''', (rs_id,))
426
snp_list = []
427
for snp in c:
428
snp_list.append(snp)
429
return snp_list
430
431
def close(self):
432
'''Closes the database.
433
'''
434
self.dbConn.execute('VACUUM')
435
self.dbConn.close()
436
437
def hasChrPopGeno(self, chr, pop):
438
'''Check if database has a certain chromosome and population (geno).
439
'''
440
c = self.dbConn.cursor()
441
c.execute('''
442
SELECT count(*)
443
FROM geno_snp
444
WHERE chromosome = ?
445
AND population = ?''', (chr, pop))
446
cnt = iter(c).next()
447
#print chr, pop, cnt
448
if cnt == (0,):
449
return False
450
else:
451
return True
452
453
def getGenepop(self, f, rsIds, pops, allowOfs = True, doMale = True):
454
"""Generate genepop file
455
allowOfs - allow offspring
456
doMale - if chr X, do Male or Female. Male is haploid
457
458
TBD: Y
459
"""
460
f.write("Data for " + str(pops) + "\n")
461
c = self.dbConn.cursor()
462
rsType = {}
463
rsIds.sort()
464
#print rsIds
465
for rsId in rsIds:
466
c.execute('''
467
SELECT chromosome, position
468
FROM freq_snp
469
WHERE rs_id = ?''', (rsId,))
470
chro, pos = iter(c).next()
471
rsType[rsId] = chro
472
f.write("%s/%s/%d\n" % (rsId, chro, pos))
473
for pop in pops:
474
f.write("Pop\n")
475
c = self.dbConn.cursor()
476
if allowOfs:
477
c.execute("select iid, gender from individual where pop=?", (pop))
478
else:
479
c.execute("select iid, gender from individual where pop=? and dad=0 and mom=0",
480
(pop,))
481
for iid, gender in c:
482
chro = rsType[rsIds[0]] #wrong
483
if chro != "X" or (gender==1 and doMale) or (gender==2 and not doMale): f.write("%s," % (iid))
484
for rsId in rsIds:
485
c = self.dbConn.cursor()
486
c.execute("select nuc1, nuc2 from geno_snp where iid=? and rs_id=?",(iid, rsId))
487
try:
488
n1, n2 = iter(c).next()
489
except StopIteration, e:
490
n1, n2 = ("N", "N")
491
chro = rsType[rsId]
492
if chro != "X":
493
f.write(" %s%s" % (self.gp(n1), self.gp(n2)))
494
elif doMale and gender==1:
495
f.write(" %s" % (self.gp(n1)))
496
elif (not doMale) and gender==2:
497
f.write(" %s%s" % (self.gp(n1), self.gp(n2)))
498
if chro != "X" or (gender==1 and doMale) or (gender==2 and not doMale): f.write("\n")
499
500
def gp(self, st):
501
if st =="N": return "00"
502
if st =="A": return "01"
503
if st =="C": return "02"
504
if st =="T": return "03"
505
if st =="G": return "04"
506
raise Exception()
Loggerhead is a web-based interface for Breezy
Version: 2.0.1