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
123
124
125
126
127
128
129
130
131
132
133
134
135
|
#!/usr/bin/env python
import h5py, sys, numpy, random
from pbsuite.honey.HSpots import *
import matplotlib.pyplot as plt
#transposon
#start = 1976517 - buffer
#end = 1977304 + buffer
#tandem duplication
#start = 1096192 - buffer
#end = 1096825 + buffer
#tandem deletion
#start = 4294209 - buffer
#end = 4294314 + buffer
#p-element inversion
#appx 8192387
#start, end = 304417-buffer, 304460+buffer #tails
#start, end = 1456604 - buffer, 1456650 + buffer
#start, end = 3255896 - buffer, 3255932 + buffer
avgWindow = numpy.ones(50, dtype=numpy.float16) / 50
def makeTransformPlots(key, data, start, end, buffer, binsize, fignum, normalize=True):
orig = numpy.convolve(data, avgWindow, "same") #smooth
if normalize:
norm = orig / cov
else:
norm = orig
slopeT = numpy.convolve(norm, slopWindow, "same") #slope transform
abs = numpy.abs(slopeT) # absolute value
smo = numpy.convolve(abs, avgWindow, "same") #smooth
space = numpy.convolve(smo * orig, avgWindow, "same") # Take back to absolute space and smooth again
win = range(start-buffer, end+buffer)
plt.figure(fignum)
plt.subplot(711)
plt.plot(win, data, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(712)
plt.plot(win, orig, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(713)
plt.plot(win, norm, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(714)
plt.plot(win, slopeT, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(715)
plt.plot(win, abs, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(716)
plt.plot(win, smo, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.subplot(717)
plt.plot(win, space, 'b-')
plt.axhline(1, color='k'); plt.axvline(start, color='k'); plt.axvline(end, color='k')
plt.show()
plt.savefig("transform_%s.png" % key)
def makeLinePlots(dataOrig, start, end, offset):
plt.figure()
print start, end, offset
print len(dataOrig)
cov = numpy.convolve(data[COV], avgWindow, "same")
print "cov", numpy.max(numpy.abs(data[COV][10:-10])), numpy.mean(cov)
insR, imu, isd = preprocessSignal(data[INS], data[COV])
insS = signalTransform(insR)
delR, dmu, dsd = preprocessSignal(data[DEL], data[COV])
delS = signalTransform(delR)
win = range(start, end)
covRg = plt.plot(win, cov, 'k-')
misRg = plt.plot(win, misR * cov, 'r-')
insRg = plt.plot(win, insR * cov, 'g-')
delRg = plt.plot(win, delR * cov, 'b-')
ticks = range(start, end, (end-start)/6)#[:-1]
labels = range(offset, offset + (end-start)+1, (end-start)/6)
print ticks, labels
plt.xticks(ticks, ticks, horizontalalignment="left", rotation=17)
plt.xlabel("position")
plt.ylabel("rate")
plt.legend([covRg, insRg, delRg], ["COV", "INS", "DEL"],)
plt.suptitle("%d bp sv (%d - %d)" % (end - start, start, end))
plt.show()
plt.savefig("rates.png")
plt.figure()
insSg = plt.plot(win, insS, 'g-')
delSg = plt.plot(win, delS, 'b-')
ticks = range(start, end+1, (end-start)/6)#[:-1]
#labels = range(offset, offset + (end-start)+1, (end-start)/6)
plt.xticks(ticks, ticks, horizontalalignment="center", rotation=17)
plt.xlabel("position")
plt.ylabel("signal")
plt.legend([insSg, delSg], ["INS", "DEL"])
plt.suptitle("%d bp sv (%d - %d)" % (end - start, start, end))
plt.show()
plt.savefig("signals.png")
if __name__ == '__main__':
h5 = h5py.File(sys.argv[1], 'r')
cols = h5.attrs["columns"]
key = sys.argv[2]
offset = int(sys.argv[3])
start = int(sys.argv[3]) - h5[key].attrs["start"]
end = int(sys.argv[4]) - h5[key].attrs["start"]
#buffer = int((start-end) *.1)
buffer = 100
start = max(0, start-buffer)
end = min(h5[key].attrs["end"], end+buffer)
data = h5[key]["data"][: , start:end]
h5.close()
makeKernals()
makeLinePlots(data, start, end, offset)
|
