/
bigwig.py
113 lines (98 loc) · 3.28 KB
/
bigwig.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
"""
Module to help create scaled bigWig files from BAM
"""
import pybedtools
import os
import subprocess
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 stderr:
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]
os.system(' '.join(cmds))
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)
p = subprocess.Popen(cmds, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
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]
os.system(' '.join(cmds))
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]
os.system(' '.join(cmds))