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
|
"""
Module to help create scaled bigWig files from BAM
"""
import pybedtools
import os
import subprocess
try:
FileNotFoundError
except NameError:
# python2
FileNotFoundError = OSError
def mapped_read_count(bam, force=False):
"""
Scale is cached in a bam.scale file containing the number of mapped reads.
Use force=True to override caching.
"""
scale_fn = bam + ".scale"
if os.path.exists(scale_fn) and not force:
for line in open(scale_fn):
if line.startswith("#"):
continue
readcount = float(line.strip())
return readcount
cmds = ["samtools", "view", "-c", "-F", "0x4", bam]
p = subprocess.Popen(cmds, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
if p.returncode:
raise ValueError("samtools says: %s" % stderr)
readcount = float(stdout)
# write to file so the next time you need the lib size you can access
# it quickly
if not os.path.exists(scale_fn):
fout = open(scale_fn, "w")
fout.write(str(readcount) + "\n")
fout.close()
return readcount
def bedgraph_to_bigwig(bedgraph, genome, output):
genome_file = pybedtools.chromsizes_to_file(pybedtools.chromsizes(genome))
cmds = ["bedGraphToBigWig", bedgraph.fn, genome_file, output]
try:
p = subprocess.Popen(cmds, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
except FileNotFoundError:
raise FileNotFoundError(
"bedGraphToBigWig was not found on the path. This is an external "
"tool from UCSC which can be downloaded from "
"http://hgdownload.soe.ucsc.edu/admin/exe/. Alternatatively, use "
"`conda install ucsc-bedgraphtobigwig`"
)
if p.returncode:
raise ValueError(
"cmds: %s\nstderr:%s\nstdout:%s" % (" ".join(cmds), stderr, stdout)
)
return output
def bigwig_to_bedgraph(fn, chrom=None, start=None, end=None, udcDir=None):
cmds = ["bigWigToBedGraph", fn]
if chrom is not None:
cmds.extend(["-chrom", chrom])
if start is not None:
cmds.extend(["-start", start])
if end is not None:
cmds.extend(["-end", end])
if udcDir is not None:
cmds.extend(["-udcDir", udcDir])
outfn = pybedtools.BedTool._tmp()
cmds.append(outfn)
try:
p = subprocess.Popen(cmds, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
except FileNotFoundError:
raise FileNotFoundError(
"bigWigToBedGraph was not found on the path. This is an external "
"tool from UCSC which can be downloaded from "
"http://hgdownload.soe.ucsc.edu/admin/exe/. Alternatatively, use "
"`conda install ucsc-bedgraphtobigwig`"
)
if p.returncode:
raise ValueError(
"cmds: %s\nstderr:%s\nstdout:%s" % (" ".join(cmds), stderr, stdout)
)
return pybedtools.BedTool(outfn)
def wig_to_bigwig(wig, genome, output):
genome_file = pybedtools.chromsizes_to_file(pybedtools.chromsizes(genome))
cmds = ["wigToBigWig", wig.fn, genome_file, output]
try:
p = subprocess.Popen(cmds, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
except FileNotFoundError:
raise FileNotFoundError(
"bigWigToBedGraph was not found on the path. This is an external "
"tool from UCSC which can be downloaded from "
"http://hgdownload.soe.ucsc.edu/admin/exe/. Alternatatively, use "
"`conda install ucsc-bedgraphtobigwig`"
)
if p.returncode:
raise ValueError(
"cmds: %s\nstderr:%s\nstdout:%s" % (" ".join(cmds), stderr, stdout)
)
return output
def bam_to_bigwig(bam, genome, output, scale=False):
"""
Given a BAM file `bam` and assembly `genome`, create a bigWig file scaled
such that the values represent scaled reads -- that is, reads per million
mapped reads.
(Disable this scaling step with scale=False; in this case values will
indicate number of reads)
Assumes that `bedGraphToBigWig` from UCSC tools is installed; see
http://genome.ucsc.edu/goldenPath/help/bigWig.html for more details on the
format.
"""
genome_file = pybedtools.chromsizes_to_file(pybedtools.chromsizes(genome))
kwargs = dict(bg=True, split=True, g=genome_file)
if scale:
readcount = mapped_read_count(bam)
_scale = 1 / (readcount / 1e6)
kwargs["scale"] = _scale
x = pybedtools.BedTool(bam).genome_coverage(**kwargs)
cmds = ["bedGraphToBigWig", x.fn, genome_file, output]
try:
p = subprocess.Popen(
cmds,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True,
)
stdout, stderr = p.communicate()
except FileNotFoundError:
raise FileNotFoundError(
"bedGraphToBigWig was not found on the path. This is an external "
"tool from UCSC which can be downloaded from "
"http://hgdownload.soe.ucsc.edu/admin/exe/. Alternatatively, use "
"`conda install ucsc-bedgraphtobigwig`"
)
if p.returncode and "bedSort" in stderr:
print("BAM header was not sorted; sorting bedGraph")
y = x.sort()
cmds[1] = y.fn
try:
p = subprocess.Popen(
cmds,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True,
)
stdout, stderr = p.communicate()
except FileNotFoundError:
raise FileNotFoundError(
"bedSort was not found on the path. This is an external "
"tool from UCSC which can be downloaded from "
"http://hgdownload.soe.ucsc.edu/admin/exe/. Alternatatively, use "
"`conda install ucsc-bedgraphtobigwig`"
)
if p.returncode:
raise ValueError(
"cmds: %s\nstderr: %s\nstdout: %s" % (" ".join(cmds), stderr, stdout)
)
|
