-
Notifications
You must be signed in to change notification settings - Fork 2
/
visualize_fisher_shifts.py
executable file
·59 lines (52 loc) · 2.35 KB
/
visualize_fisher_shifts.py
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
#!/usr/bin/env python
# encoding: utf-8
"""
From the shifts and merged sites, create a bed12 for each comparison with a
trackline to visualize shift in UCSC browser. Unsorted output to stdout.
"""
from toolshed import reader
from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
class Bed(object):
def __init__(self, toks):
self.chrom = toks[0]
self.start, self.stop = map(int, toks[1:3])
self.name, self.score, self.strand = toks[3:6]
def bed12line(chrom, start, stop, strand, shift):
score = 0
color = "255,0,0"
diff = stop - start - 1
if strand == "+":
arrow = "-" if shift == "proximal" else "+"
if strand == "-":
arrow = "+" if shift == "proximal" else "-"
return [chrom, start, stop, shift, score, arrow, start, stop, color, 2, "1,1", "0,{diff}".format(diff=diff)]
def sites_to_dict(bed):
d = {}
for l in reader(bed, header=Bed):
d[l.name] = l
return d
def main(shifts, sites, cutoff=0.05, min_length=5):
refsites = sites_to_dict(sites)
for l in reader(shifts):
if not l['shift'] == "proximal" and not l['shift'] == "distal": continue
if float(l['q']) >= cutoff: continue
sites = [l['SiteA'], l['SiteB']]
# sort by ascending site IDs
sites.sort(key=lambda x: int(x.split(".")[-1]))
downstream, upstream = sites
a = refsites[downstream]
b = refsites[upstream]
if abs(a.start - b.start) < min_length: continue
# if the strand is negative, sites are classified in descending order
if a.strand == "-":
print "\t".join(map(str, bed12line(a.chrom, b.start, a.stop, a.strand, l['shift'])))
else:
print "\t".join(map(str, bed12line(a.chrom, a.start, b.stop, a.strand, l['shift'])))
if __name__ == '__main__':
p = ArgumentParser(description=__doc__, formatter_class=ArgumentDefaultsHelpFormatter)
p.add_argument("shifts", help="output of `fisher_test.py`")
p.add_argument("sites", help="output of `merge_sites.py`")
p.add_argument("-q", dest="cutoff", default=0.05, type=float, help="suppress peaks with a q-value less than cutoff")
p.add_argument("-l", dest="min_length", default=5, type=int, help="suppress peaks with a distance between them less than min_length")
args = p.parse_args()
main(args.shifts, args.sites, args.cutoff, args.min_length)