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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
|
#!/usr/bin/python3
# -*- coding: utf-8 -*-
import argparse
import sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['svg.fonttype'] = 'none'
from deeptools import cm # noqa: F401
import matplotlib.pyplot as plt
import plotly.offline as py
import plotly.graph_objs as go
# own tools
from deeptools.parserCommon import writableFile
from deeptools.getFragmentAndReadSize import get_read_and_fragment_length
from deeptools._version import __version__
def parse_arguments():
parser = argparse.ArgumentParser(
description='This tool calculates the fragment sizes for read pairs given a BAM file from paired-end sequencing.'
'Several regions are sampled depending on the '
'size of the genome and number of processors to estimate the'
'summary statistics on the fragment lengths. '
'Properly paired reads are preferred for computation, i.e., '
'it will only use discordant pairs if no concordant alignments '
'overlap with a given region. '
'The default setting simply prints the summary statistics to the screen.')
parser.add_argument('--bamfiles', '-b',
help='List of BAM files to process',
nargs='+',
metavar='bam files')
parser.add_argument('--histogram', '-hist', '-o',
help='Save a .png file with a histogram '
'of the fragment length distribution.',
metavar='FILE')
parser.add_argument('--plotFileFormat',
metavar='FILETYPE',
help='Image format type. If given, this option '
'overrides the image format based on the plotFile '
'ending. The available options are: png, '
'eps, pdf, svg and plotly.',
default=None,
choices=['png', 'pdf', 'svg', 'eps', 'plotly'])
parser.add_argument('--numberOfProcessors', '-p',
help='Number of processors to use. The default is '
'to use 1. (Default: %(default)s)',
metavar="INT",
type=int,
default=1,
required=False)
parser.add_argument('--samplesLabel',
help='Labels for the samples plotted. The '
'default is to use the file name of the '
'sample. The sample labels should be separated '
'by spaces and quoted if a label itself'
'contains a space E.g. --samplesLabel label-1 "label 2" ',
nargs='+')
parser.add_argument('--plotTitle', '-T',
help='Title of the plot, to be printed on top of '
'the generated image. Leave blank for no title. (Default: %(default)s)',
default='')
parser.add_argument('--maxFragmentLength',
help='The maximum fragment length in the histogram. A value of 0 (the default) indicates to use twice the mean fragment length. (Default: %(default)s)',
default=0,
type=int)
parser.add_argument('--logScale',
help='Plot on the log scale',
action='store_true')
parser.add_argument('--binSize', '-bs',
metavar='INT',
help='Length in bases of the window used to sample the genome. (Default: %(default)s)',
default=1000,
type=int)
parser.add_argument('--distanceBetweenBins', '-n',
metavar='INT',
help='To reduce the computation time, not every possible genomic '
'bin is sampled. This option allows you to set the distance '
'between bins actually sampled from. Larger numbers are sufficient '
'for high coverage samples, while smaller values are useful for '
'lower coverage samples. Note that if you specify a value that '
'results in too few (<1000) reads sampled, the value will be '
'decreased. (Default: %(default)s)',
default=1000000,
type=int)
parser.add_argument('--blackListFileName', '-bl',
help="A BED file containing regions that should be excluded from all analyses. Currently this works by rejecting genomic chunks that happen to overlap an entry. Consequently, for BAM files, if a read partially overlaps a blacklisted region or a fragment spans over it, then the read/fragment might still be considered.",
metavar="BED file",
required=False)
parser.add_argument('--table',
metavar='FILE',
help='In addition to printing read and fragment length metrics to the screen, write them to the given file in tabular format.',
required=False)
parser.add_argument('--outRawFragmentLengths',
metavar='FILE',
required=False,
type=writableFile,
help='Save the fragment (or read if the input is single-end) length and their associated number of occurrences to a tab-separated file. Columns are length, number of occurrences, and the sample label.')
parser.add_argument('--verbose',
help='Set if processing data messages are wanted.',
action='store_true',
required=False)
parser.add_argument('--version', action='version',
version='%(prog)s {}'.format(__version__))
return parser
def getDensity(lengths, minVal, maxVal):
"""
This is essentially computing what hist() in matplotlib is doing and returning the results.
This then allows us to free up the memory consumed by each sample rather than returning it all back to main() for plotting.
"""
n, bins, patches = plt.hist(lengths, bins=100, range=(minVal, maxVal), density=True)
plt.clf()
return (n, bins)
def getFragSize(bam, args, idx, outRawFrags):
fragment_len_dict, read_len_dict = get_read_and_fragment_length(bam, return_lengths=True,
blackListFileName=args.blackListFileName,
numberOfProcessors=args.numberOfProcessors,
verbose=args.verbose,
binSize=args.binSize,
distanceBetweenBins=args.distanceBetweenBins)
if outRawFrags:
label = bam
if args.samplesLabel and idx < len(args.samplesLabel):
label = args.samplesLabel[idx]
if fragment_len_dict:
fragment_len_dict['lengths'] = [int(x) for x in fragment_len_dict['lengths']]
cnts = np.bincount(fragment_len_dict['lengths'], minlength=int(fragment_len_dict['max']) + 1)
else:
read_len_dict['lengths'] = [int(x) for x in read_len_dict['lengths']]
cnts = np.bincount(read_len_dict['lengths'], minlength=int(read_len_dict['max']) + 1)
for idx, v in enumerate(cnts):
if v > 0:
outRawFrags.write("{}\t{}\t{}\n".format(idx, v, label))
if args.samplesLabel and idx < len(args.samplesLabel):
print("\n\nSample label: {}".format(args.samplesLabel[idx]))
else:
print("\n\nBAM file : {}".format(bam))
if fragment_len_dict:
if fragment_len_dict['mean'] == 0:
print("No pairs were found. Is the data from a paired-end sequencing experiment?")
print("Sample size: {}\n".format(fragment_len_dict['sample_size']))
print("Fragment lengths:")
print("Min.: {}\n1st Qu.: {}\nMean: {}\nMedian: {}\n"
"3rd Qu.: {}\nMax.: {}\nStd: {}".format(fragment_len_dict['min'],
fragment_len_dict['qtile25'],
fragment_len_dict['mean'],
fragment_len_dict['median'],
fragment_len_dict['qtile75'],
fragment_len_dict['max'],
fragment_len_dict['std']))
print("MAD: {}\nLen. 10%: {}\nLen. 20%: {}\nLen. 30%: {}\nLen. 40%: {}\nLen. 60%: {}\nLen. 70%: {}\nLen. 80%: {}\nLen. 90%: {}\nLen. 99%: {}\n".format(fragment_len_dict['mad'],
fragment_len_dict['qtile10'],
fragment_len_dict['qtile20'],
fragment_len_dict['qtile30'],
fragment_len_dict['qtile40'],
fragment_len_dict['qtile60'],
fragment_len_dict['qtile70'],
fragment_len_dict['qtile80'],
fragment_len_dict['qtile90'],
fragment_len_dict['qtile99']))
else:
print("No pairs were found. Is the data from a paired-end sequencing experiment?")
print("\nRead lengths:")
print("Sample size: {}\n".format(read_len_dict['sample_size']))
print("Min.: {}\n1st Qu.: {}\nMean: {}\nMedian: {}\n"
"3rd Qu.: {}\nMax.: {}\nStd: {}".format(read_len_dict['min'],
read_len_dict['qtile25'],
read_len_dict['mean'],
read_len_dict['median'],
read_len_dict['qtile75'],
read_len_dict['max'],
read_len_dict['std']))
print("MAD: {}\nLen. 10%: {}\nLen. 20%: {}\nLen. 30%: {}\nLen. 40%: {}\nLen. 60%: {}\nLen. 70%: {}\nLen. 80%: {}\nLen. 90%: {}\nLen. 99%: {}\n".format(read_len_dict['mad'],
read_len_dict['qtile10'],
read_len_dict['qtile20'],
read_len_dict['qtile30'],
read_len_dict['qtile40'],
read_len_dict['qtile60'],
read_len_dict['qtile70'],
read_len_dict['qtile80'],
read_len_dict['qtile90'],
read_len_dict['qtile99']))
# The read and fragment lists will just eat up memory if not removed!
if args.histogram:
if fragment_len_dict:
maxVal = fragment_len_dict['mean'] * 2
minVal = fragment_len_dict['min']
else:
maxVal = read_len_dict['mean'] * 2
minVal = read_len_dict['min']
if args.maxFragmentLength > 0:
maxVal = args.maxFragmentLength
if fragment_len_dict:
fragment_len_dict['lengths'] = getDensity(fragment_len_dict['lengths'], minVal, maxVal)
if read_len_dict:
read_len_dict['lengths'] = getDensity(read_len_dict['lengths'], minVal, maxVal)
else:
if fragment_len_dict:
del fragment_len_dict['lengths']
if read_len_dict:
del read_len_dict['lengths']
return (fragment_len_dict, read_len_dict)
def printTable(args, fragDict, readDict):
"""
Print the read and fragment dictionary in more easily parsable tabular format to a file.
"""
of = open(args.table, "w")
of.write("\tFrag. Sampled")
of.write("\tFrag. Len. Min.\tFrag. Len. 1st. Qu.\tFrag. Len. Mean\tFrag. Len. Median\tFrag. Len. 3rd Qu.\tFrag. Len. Max\tFrag. Len. Std.")
of.write("\tFrag. Med. Abs. Dev.\tFrag. Len. 10%\tFrag. Len. 20%\tFrag. Len. 30%\tFrag. Len. 40%\tFrag. Len. 60%\tFrag. Len. 70%\tFrag. Len. 80%\tFrag. Len. 90%\tFrag. Len. 99%")
of.write("\tReads Sampled")
of.write("\tRead Len. Min.\tRead Len. 1st. Qu.\tRead Len. Mean\tRead Len. Median\tRead Len. 3rd Qu.\tRead Len. Max\tRead Len. Std.")
of.write("\tRead Med. Abs. Dev.\tRead Len. 10%\tRead Len. 20%\tRead Len. 30%\tRead Len. 40%\tRead Len. 60%\tRead Len. 70%\tRead Len. 80%\tRead Len. 90%\tRead Len. 99%\n")
for idx, bam in enumerate(args.bamfiles):
if args.samplesLabel and idx < len(args.samplesLabel):
of.write(args.samplesLabel[idx])
else:
of.write(bam)
if fragDict is not None and fragDict[bam] is not None:
d = fragDict[bam]
of.write("\t{}".format(d['sample_size']))
of.write("\t{}\t{}\t{}\t{}\t{}\t{}\t{}".format(d['min'],
d['qtile25'],
d['mean'],
d['median'],
d['qtile75'],
d['max'],
d['std']))
of.write("\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}".format(d['mad'],
d['qtile10'],
d['qtile20'],
d['qtile30'],
d['qtile40'],
d['qtile60'],
d['qtile70'],
d['qtile80'],
d['qtile90'],
d['qtile99']))
else:
of.write("\t0")
of.write("\t0\t0\t0\t0\t0\t0\t0")
of.write("\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0")
d = readDict[bam]
of.write("\t{}".format(d['sample_size']))
of.write("\t{}\t{}\t{}\t{}\t{}\t{}\t{}".format(d['min'],
d['qtile25'],
d['mean'],
d['median'],
d['qtile75'],
d['max'],
d['std']))
of.write("\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\n".format(d['mad'],
d['qtile10'],
d['qtile20'],
d['qtile30'],
d['qtile40'],
d['qtile60'],
d['qtile70'],
d['qtile80'],
d['qtile90'],
d['qtile99']))
of.close()
def main(args=None):
args = parse_arguments().parse_args(args)
fraglengths = {}
readlengths = {}
of = None
if args.outRawFragmentLengths is not None:
of = open(args.outRawFragmentLengths, "w")
of.write("#bamPEFragmentSize\nSize\tOccurrences\tSample\n")
for idx, bam in enumerate(args.bamfiles):
f, r = getFragSize(bam, args, idx, of)
fraglengths[bam] = f
readlengths[bam] = r
if args.table is not None:
printTable(args, fraglengths, readlengths)
if args.histogram:
if args.samplesLabel:
if len(args.bamfiles) != len(args.samplesLabel):
sys.exit("The number of labels does not match the number of BAM files.")
else:
labels = args.samplesLabel
else:
labels = list(fraglengths.keys())
i = 0
data = []
for bam in fraglengths.keys():
d = fraglengths[bam]
if d is None:
d = readlengths[bam]
if args.maxFragmentLength > 0:
maxVal = args.maxFragmentLength
else:
maxVal = d['mean'] * 2
if args.plotFileFormat == 'plotly':
trace = go.Histogram(x=d['lengths'],
histnorm='probability',
opacity=0.5,
name=labels[i],
nbinsx=100,
xbins=dict(start=d['min'], end=maxVal))
data.append(trace)
else:
plt.bar(d['lengths'][1][:-1], height=d['lengths'][0],
width=d['lengths'][1][1:] - d['lengths'][1][:-1],
align='edge', log=args.logScale,
alpha=0.5, label=labels[i])
i += 1
if args.plotFileFormat == 'plotly':
fig = go.Figure()
fig['data'] = data
fig['layout']['yaxis1'].update(title='Frequency')
fig['layout']['xaxis1'].update(title='Fragment Length')
fig['layout'].update(title=args.plotTitle)
fig['layout'].update(showlegend=True)
if args.logScale:
fig['layout']['yaxis1'].update(type='log')
py.plot(fig, filename=args.histogram, auto_open=False)
else:
plt.xlabel('Fragment Length')
plt.ylabel('Frequency')
plt.legend(loc='upper right')
plt.title(args.plotTitle)
plt.savefig(args.histogram, bbox_inches=0, format=args.plotFileFormat)
plt.close()
if __name__ == "__main__":
main()
|
