forked from chiulab/surpi
-
Notifications
You must be signed in to change notification settings - Fork 0
/
coveragePlot.py
executable file
·125 lines (108 loc) · 3.5 KB
/
coveragePlot.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
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
#!/usr/bin/python
#
# coveragePlot.py
#
# This program generates genomic coverage plots
# Chiu Laboratory
# University of California, San Francisco
# January, 2014
#
# Copyright (C) 2014 Charles Y Chiu - All Rights Reserved
# SURPI has been released under a modified BSD license.
# Please see license file for details.
# Last revised 1/26/2014
#
# Added option 3/22/2014 for log y-axes
import matplotlib
#import matplotlib.numerix as nx
matplotlib.use('Agg')
from pylab import *
from pylab import figure, show, legend
from matplotlib import pyplot as plt
import numpy as np
import sys, os
import re
def smart_truncate1(text, max_length=100, suffix='...'):
"""Returns a string of at most `max_length` characters, cutting
only at word-boundaries. If the string was truncated, `suffix`
will be appended.
"""
if len(text) > max_length:
pattern = r'^(.{0,%d}\S)\s.*' % (max_length-len(suffix)-1)
return re.sub(pattern, r'\1' + suffix, text)
else:
return text
if len(sys.argv) < 3:
print "usage: coveragePlot.py <data file .map/.report> <title of plot> <log y-axes Y/N/B=both>"
sys.exit(-1)
dataFile = sys.argv[1]
data = mlab.load(dataFile)
outputFile = os.path.splitext(dataFile)[0]+".ps"
reportFile = os.path.splitext(dataFile)[0]+".report"
with open(reportFile) as f:
reportContent = f.readlines()
reportText = ""
logPlot = sys.argv[3]
for line in reportContent:
stripped_line = line.rstrip('\r\n\t ')
reportText = reportText + smart_truncate1(stripped_line, max_length=100, suffix='...') + "\n"
print "Loaded " + dataFile
hold(True)
if logPlot=='N':
fig=plt.figure(figsize=[8.5,4.5])
ax = fig.add_subplot(111)
fig.text(0.1,0.0,reportText, fontsize=9)
color ='k-'
plot(data[:,0],data[:,1],color)
xlabel("base position",fontsize=8)
ylabel("fold coverage",fontsize=8)
title_text = sys.argv[2]
suptitle(title_text,fontsize=9)
xMin, xMax, yMin, yMax = min(data[:,0]),max(data[:,0]),min(data[:,1]),max(data[:,1])
# add a 10% buffer to yMax
yMax *= 1.1
axis([xMin,xMax,yMin,yMax])
gcf().subplots_adjust(bottom=0.60)
plt.show()
if logPlot=='B':
fig=plt.figure(figsize=[8.5,4.5])
ax1 = fig.add_subplot(211)
color ='k-'
plot(data[:,0],data[:,1],color)
xlabel("base position",fontsize=8)
ylabel("fold coverage",fontsize=8)
xMin, xMax, yMin, yMax = min(data[:,0]),max(data[:,0]),min(data[:,1]),max(data[:,1])
yMax *= 1.1
axis([xMin,xMax,yMin,yMax])
plt.show()
ax2 = fig.add_subplot(212)
ax2.set_yscale('symlog')
fig.text(0.1,0.0,reportText, fontsize=9)
color ='k-'
plot(data[:,0],data[:,1],color)
xlabel("base position",fontsize=8)
ylabel("fold coverage",fontsize=8)
title_text = sys.argv[2]
suptitle(title_text,fontsize=9)
xMin, xMax, yMin, yMax = min(data[:,0]),max(data[:,0]),min(data[:,1]),max(data[:,1])
yMax *= 1.1
axis([xMin,xMax,yMin,yMax])
gcf().subplots_adjust(bottom=0.40)
plt.show()
if logPlot=='Y':
fig=plt.figure(figsize=[8.5,4.5])
ax = fig.add_subplot(111)
ax.set_yscale('symlog')
fig.text(0.1,0.0,reportText, fontsize=9)
color ='k-'
plot(data[:,0],data[:,1],color)
xlabel("base position",fontsize=8)
ylabel("fold coverage",fontsize=8)
title_text = sys.argv[2]
suptitle(title_text,fontsize=9)
xMin, xMax, yMin, yMax = min(data[:,0]),max(data[:,0]),min(data[:,1]),max(data[:,1])
yMax *= 1.1
axis([xMin,xMax,yMin,yMax])
gcf().subplots_adjust(bottom=0.60)
plt.show()
savefig(outputFile)