github_scripts/NS_domain_SNP_hist.py at master · nick297/github_scripts · GitHub

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
#!/usr/bin/env python
import sys
import matplotlib.pyplot as plt
import numpy as np

bed_file = sys.argv[1]


class bedDetail:
    def __init__(self, GStart, GStop, EStart, EStop, GeneName, ExonName):
        self.GS = GStart
        self.GP = GStop
        self.ES = EStart
        self.EP = EStop
        self.GN = GeneName
        self.EN = ExonName


beds = []
for line in open(bed_file, 'r'):
    l = line.split("\t")
    gstart = int(l[1])
    gstop = int(l[2])
    estart = int(l[-2])
    estop = int(l[-1])
    feature = str(l[3])
    f = feature.split(" ")
    gene = str(f[0])
    d = f[1:]
    domain = str(' '.join(d))
    x = bedDetail(gstart, gstop, estart, estop, gene, domain)
    beds.append(x)


nDomains = len(beds)
#print nDomains
results_file = sys.argv[2]

names = []
results = []
for line in open(results_file, "rU"):
    line_results = line.split("\t")
    resultsL = []
    if line.startswith('H'):
        for l in line_results:
            names.append(l.replace("\n",""))
    else:
        for l in line_results:
            resultsL.append(int(l))
        results.append(resultsL)

#print results
#print names
animals = len(results[0]) # calculate number of animals
ind = np.arange(nDomains)    # the x locations for the groups
width = 0.90       # the width of the bars: can also be len(x) sequence
ax = plt.subplot2grid((24,12), (0,0), rowspan=18, colspan=11)
score50 = []

for f in results:
    score50.append(0)


#colors = ["#FF0000", "#009999", "#9FEE00", "#BF3030", "#1D7373", "#86B32D", "#A60000",\
#"#006363", "#679B00", "#FF4040", "#33CCCC", "#B9F73E", "#FF7373", "#5CCCCC", "#C9F76F"]


colors = ["k", \
"#a6cee3", \
"#1f78b4", \
"0.8", \
"#b2df8a", \
"#33a02c", \
"#fb9a99", \
"0.6", \
"#e31a1c", \
"#fdbf6f", \
"#ff7f00", \
"#cab2d6", \
"#6a3d9a", \
"#ffff99", \
"#a6cee3", \
"#1f78b4", \
"0.8", \
"#b2df8a", \
"#33a02c", \
"#fb9a99", \
"0.6", \
"#e31a1c", \
"#fdbf6f", \
"#ff7f00", \
"#cab2d6", \
"#6a3d9a", \
"#ffff99", \
"#b15928"]


for i in range(animals):
    pi = "p" + str(i)
    score75 = []
    for f in results:
        score75.append(f[i])
    #print pi, score75
    pi = ax.bar(ind, score75, width, bottom=score50,color=colors[i], label=names[i])
    lists = []
    lists.append(score75)
    lists.append(score50)
    score50 = map(sum, zip(*lists))

plt.ylabel('Number of SNPs per domain')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[::-1], labels[::-1], bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., prop={'size':11})
ax.tick_params(axis='x', which='major', labelsize=8)
ax.tick_params(axis='x', which='minor', labelsize=6)
plt.xlim([0,len(score75)])
domains = []
for b in beds:
    domains.append(str(b.EN))
plt.xticks(ind+width/2., (domains), rotation=90)
plt.show()
#plt.savefig(savename, format='pdf')