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
|
#!/usr/bin/python3
# -*- coding: utf-8 -*-
import os
import sys
import argparse
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
from importlib.metadata import version
import deeptools.countReadsPerBin as countR
from deeptools import parserCommon
from deeptools.utilities import smartLabels
old_settings = np.seterr(all='ignore')
def parse_arguments(args=None):
parent_parser = parserCommon.getParentArgParse(binSize=False)
read_options_parser = parserCommon.read_options()
parser = \
argparse.ArgumentParser(
parents=[required_args(), parent_parser, read_options_parser],
formatter_class=argparse.RawDescriptionHelpFormatter,
add_help=False,
description="""
This tool is useful to assess the sequencing depth of a given sample.
It samples 1 million bp, counts the number of overlapping reads and can report
a histogram that tells you how many bases are covered how many times.
Multiple BAM files are accepted, but they all should correspond to the same genome assembly.
detailed usage help:
$ plotCoverage -h
""",
epilog='example usages:\nplotCoverage '
'--bamfiles file1.bam file2.bam -o results.png\n\n'
' \n\n',
conflict_handler='resolve',
usage='plotCoverage -b sample1.bam sample2.bam -o coverage.png \n'
'help: plotCoverage -h / plotCoverage --help\n')
parser.add_argument('--version', action='version',
version='%(prog)s {}'.format(version('deeptools')))
return parser
def process_args(args=None):
args = parse_arguments().parse_args(args)
if not args.labels:
if args.smartLabels:
args.labels = smartLabels(args.bamfiles)
else:
args.labels = [os.path.basename(x) for x in args.bamfiles]
if args.labels and len(args.bamfiles) != len(args.labels):
sys.exit("The number of labels does not match the number of BAM files.")
return args
def required_args():
parser = argparse.ArgumentParser(add_help=False)
required = parser.add_argument_group('Required arguments')
required.add_argument('--bamfiles', '-b',
metavar='FILE1 FILE2',
help='List of indexed BAM files separated by spaces.',
nargs='+',
required=True)
optional = parser.add_argument_group('Optional arguments')
optional.add_argument("--help", "-h", action="help",
help="show this help message and exit")
optional.add_argument('--plotFile', '-o',
type=parserCommon.writableFile,
help='File name to save the plot to.')
optional.add_argument('--labels', '-l',
metavar='sample1 sample2',
help='User defined labels instead of default labels from '
'file names. '
'Multiple labels have to be separated by spaces, e.g. '
'--labels sample1 sample2 sample3',
nargs='+')
optional.add_argument('--smartLabels',
action='store_true',
help='Instead of manually specifying labels for the input '
'BAM files, this causes deepTools to use the file name '
'after removing the path and extension.')
optional.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='')
optional.add_argument('--skipZeros',
help='By setting this option, genomic regions '
'that have zero or nan values in _all_ samples '
'are excluded.',
action='store_true',
required=False)
optional.add_argument('--numberOfSamples', '-n',
help='Number of 1 bp regions to sample. (Default: %(default)s)',
required=False,
type=int,
default=1000000)
optional.add_argument('--BED',
help='Limits the coverage analysis to '
'the regions specified in these files. This overrides --numberOfSamples. '
'Due to memory requirements, it is inadvised to combine this with '
'--outRawCounts or many tens of thousands of regions, as per-base '
'coverage is used!',
metavar='FILE1.bed FILE2.bed',
nargs='+')
optional.add_argument('--outRawCounts',
help='Save raw counts (coverages) to file.',
type=parserCommon.writableFile,
metavar='FILE')
optional.add_argument('--outCoverageMetrics',
help='Save percentage of bins/regions above the specified thresholds to '
'the specified file. The coverage thresholds are specified by '
'--coverageThresholds. If no coverage thresholds are specified, the file '
'will be empty.',
type=parserCommon.writableFile,
metavar='FILE')
optional.add_argument('--coverageThresholds', '-ct',
type=int,
action="append",
help='The percentage of reported bins/regions with signal at least as '
'high as the given threshold. This can be specified multiple times.')
optional.add_argument('--plotHeight',
help='Plot height in cm. (Default: %(default)s)',
type=float,
default=5.0)
optional.add_argument('--plotWidth',
help='Plot width in cm. The minimum value is 1 cm. (Default: %(default)s)',
type=float,
default=15.0)
optional.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'])
return parser
def main(args=None):
args = process_args(args)
if not args.outRawCounts and not args.plotFile and not args.outCoverageMetrics:
sys.exit("At least one of --plotFile, --outRawCounts and --outCoverageMetrics are required.\n")
if 'BED' in args:
bed_regions = args.BED
else:
bed_regions = None
cr = countR.CountReadsPerBin(args.bamfiles,
binLength=1,
bedFile=bed_regions,
numberOfSamples=args.numberOfSamples,
numberOfProcessors=args.numberOfProcessors,
verbose=args.verbose,
region=args.region,
blackListFileName=args.blackListFileName,
extendReads=args.extendReads,
minMappingQuality=args.minMappingQuality,
ignoreDuplicates=args.ignoreDuplicates,
center_read=args.centerReads,
samFlag_include=args.samFlagInclude,
samFlag_exclude=args.samFlagExclude,
minFragmentLength=args.minFragmentLength,
maxFragmentLength=args.maxFragmentLength,
bed_and_bin=True,
out_file_for_raw_data=args.outRawCounts)
num_reads_per_bin = cr.run()
if args.outCoverageMetrics and args.coverageThresholds:
args.coverageThresholds.sort() # Galaxy in particular tends to give things in a weird order
of = open(args.outCoverageMetrics, "w")
of.write("Sample\tThreshold\tPercent\n")
nbins = float(num_reads_per_bin.shape[0])
for thresh in args.coverageThresholds:
vals = np.sum(num_reads_per_bin >= thresh, axis=0)
for lab, val in zip(args.labels, vals):
of.write("{}\t{}\t{:6.3f}\n".format(lab, thresh, 100. * val / nbins))
of.close()
if args.outRawCounts:
# append to the generated file the
# labels
header = "#plotCoverage --outRawCounts\n#'chr'\t'start'\t'end'\t"
header += "'" + "'\t'".join(args.labels) + "'\n"
f = open(args.outRawCounts, 'r+')
content = f.read()
f.seek(0, 0)
f.write(header + content)
f.close()
if num_reads_per_bin.shape[0] < 2:
exit("ERROR: too few non-zero bins found.\n"
"If using --region please check that this "
"region is covered by reads.\n")
if args.skipZeros:
num_reads_per_bin = countR.remove_row_of_zeros(num_reads_per_bin)
if args.plotFile:
if args.plotFileFormat == 'plotly':
fig = go.Figure()
fig['layout']['xaxis1'] = {'domain': [0.0, 0.48], 'anchor': 'x1', 'title': 'coverage (#reads per base)'}
fig['layout']['xaxis2'] = {'domain': [0.52, 1.0], 'anchor': 'x2', 'title': 'coverage (#reads per base)'}
fig['layout']['yaxis1'] = {'domain': [0.0, 1.0], 'anchor': 'x1', 'title': 'fraction of bases sampled'}
fig['layout']['yaxis2'] = {'domain': [0.0, 1.0], 'anchor': 'x2', 'title': 'fraction of bases sampled >= coverage'}
fig['layout'].update(title=args.plotTitle)
else:
fig, axs = plt.subplots(1, 2, figsize=(args.plotWidth, args.plotHeight))
plt.suptitle(args.plotTitle)
# plot up to two std from mean
num_reads_per_bin = num_reads_per_bin.astype(int)
sample_mean = num_reads_per_bin.mean(axis=0)
sample_std = num_reads_per_bin.std(axis=0)
sample_max = num_reads_per_bin.max(axis=0)
sample_min = num_reads_per_bin.min(axis=0)
sample_25 = np.percentile(num_reads_per_bin, 25, axis=0)
sample_50 = np.percentile(num_reads_per_bin, 50, axis=0)
sample_75 = np.percentile(num_reads_per_bin, 75, axis=0)
# use the largest 99th percentile from all samples to set the x_max value
x_max = np.max(np.percentile(num_reads_per_bin, 99, axis=0))
# plot coverage
# print headers for text output
print("sample\tmean\tstd\tmin\t25%\t50%\t75%\tmax")
# the determination of a sensible value for y_max of the first plot (fraction of bases sampled vs.
# coverage) is important because, depending on the data,
# it becomes very difficult to see the lines in the plot. For example, if the coverage of a sample
# is a nice gaussian curve with a large mean of 50. Then a sensible range for the y axis (fraction of
# reads having coverage=x) is (0, 0.02) which nicely shows the coverage curve. If instead the coverage is
# very por and centers close to 1 then a good y axis range is (0,1).
# the current implementation aims to find the y_value for which 50% of the reads >= x (coverage) and
# sets that as the x_axis range.
y_max = []
data = []
# We need to manually set the line colors so they're shared between the two plots.
plotly_colors = ["#d73027", "#fc8d59", "#f33090", "#e0f3f8", "#91bfdb", "#4575b4"]
plotly_styles = sum([6 * ["solid"], 6 * ["dot"], 6 * ["dash"], 6 * ["longdash"], 6 * ["dashdot"], 6 * ["longdashdot"]], [])
for idx, col in enumerate(num_reads_per_bin.T):
if args.plotFile:
frac_reads_per_coverage = np.bincount(col.astype(int)).astype(float) / num_reads_per_bin.shape[0]
csum = np.bincount(col.astype(int))[::-1].cumsum()
csum_frac = csum.astype(float)[::-1] / csum.max()
if args.plotFileFormat == 'plotly':
color = plotly_colors[idx % len(plotly_colors)]
dash = plotly_styles[idx % len(plotly_styles)]
trace = go.Scatter(x=np.arange(0, int(x_max) - 1),
y=frac_reads_per_coverage[:int(x_max)],
mode='lines',
xaxis='x1',
yaxis='y1',
line=dict(color=color, dash=dash),
name="{}, mean={:.1f}".format(args.labels[idx], sample_mean[idx]),
legendgroup="{}".format(idx))
data.append(trace)
trace = go.Scatter(x=np.arange(0, int(x_max) - 1),
y=csum_frac[:int(x_max)],
mode='lines',
xaxis='x2',
yaxis='y2',
line=dict(color=color, dash=dash),
name=args.labels[idx],
showlegend=False,
legendgroup="{}".format(idx))
data.append(trace)
else:
axs[0].plot(frac_reads_per_coverage, label="{}, mean={:.1f}".format(args.labels[idx], sample_mean[idx]))
axs[1].plot(csum_frac, label=args.labels[idx])
# find the indexes (i.e. the x values) for which the cumulative distribution 'fraction of bases
# sampled >= coverage' where fraction of bases sampled = 50%: `np.flatnonzero(csum_frac>0.5)`
# then find the fraction of bases sampled that that have the largest x
y_max.append(frac_reads_per_coverage[max(np.flatnonzero(csum_frac > 0.5))])
print("{}\t{:0.2f}\t{:0.2f}\t{}\t{}\t{}\t{}\t{}\t".format(args.labels[idx],
sample_mean[idx],
sample_std[idx],
sample_min[idx],
sample_25[idx],
sample_50[idx],
sample_75[idx],
sample_max[idx],
))
if args.plotFile:
# Don't clip plots
y_max = max(y_max)
if args.plotFileFormat == "plotly":
fig.add_traces(data)
fig['layout']['yaxis1'].update(range=[0.0, min(1, y_max + (y_max * 0.10))])
fig['layout']['yaxis2'].update(range=[0.0, 1.0])
py.plot(fig, filename=args.plotFile, auto_open=False)
else:
axs[0].set_ylim(0, min(1, y_max + (y_max * 0.10)))
axs[0].set_xlim(0, x_max)
axs[0].set_xlabel('coverage (#reads per bp)')
axs[0].legend(fancybox=True, framealpha=0.5)
axs[0].set_ylabel('fraction of bases sampled')
# plot cumulative coverage
axs[1].set_xlim(0, x_max)
axs[1].set_xlabel('coverage (#reads per bp)')
axs[1].set_ylabel('fraction of bases sampled >= coverage')
axs[1].legend(fancybox=True, framealpha=0.5)
plt.savefig(args.plotFile, format=args.plotFileFormat)
plt.close()
if __name__ == "__main__":
main()
|