blob: 405b3dca1d274c095ab91e7a07b226138716e251 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
# Usage: python make_parameter_dapseq3.py
# Because dap-seq dose not include all TFs, so include other TFs. The idea is that TFs within a same family are very conservative in binding.
#
import sys, glob, os, operator
from geneid2name import make_gene_name_AGI_map_dict
DAPSEQ_DIR = '/home/hui/network/dapseq_merged'
MAP_FILE = '/home/hui/network/dapseq_merged/tffamily.simple.txt'
def get_name(s):
lst = s.split('_')
result = []
for x in lst:
if x != 'tnt':
result.append(x)
return '_'.join(result)
def make_dapseq_dictionary(dirname):
d = {}
files = glob.glob(os.path.join(dirname, '*.narrowPeak'))
for f in files:
lst = f.split('/')
tf_name = lst[-1].split('.')[0]
tf_name = get_name(tf_name)
if not tf_name in d:
d[tf_name] = f
else:
print('ERROR: transcription factor name not unique.')
sys.exit()
return d
d = make_dapseq_dictionary(DAPSEQ_DIR)
f = open(MAP_FILE)
lines = f.readlines()
f.close()
# since MAP_FILE contain duplicate lines
d_family = {}
for line in lines:
line = line.strip()
lst = line.split()
tf = lst[0].upper()
tf_name = lst[1]
family = lst[2]
if not tf in d_family:
d_family[tf] = (tf_name, family)
else:
if family != d_family[tf][1]:
print('WARNING: %s conflict [%s %s]!' % (tf, family, d_family[tf][1]))
count = 1
for k in sorted(d_family.keys()):
g = k
gname = d_family[k][0]
key = d_family[k][1]
if key in d:
cid = 'C0003%09d' % (count)
count += 1
print('@%s' % (cid))
print('PROTEIN_ID:%s' % (g))
print('PROTEIN_NAME:%s' % (gname))
print('DATA_NAME:%s' % (gname))
print('DATA_FORMAT:%s' % ('narrowPeak'))
print('DESCRIPTION:inferred from dapseq')
#print('LOCATION:%s' % (os.path.join(DAPSEQ_DIR, d[key])))
print('LOCATION:%s' % (d[key]))
print('NOTE:')
print('')
|
