-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathboiler.py
executable file
·146 lines (123 loc) · 7.19 KB
/
boiler.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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#! /usr/bin/env python
import sys
import argparse
import time
import logging
import os
VERSION = '1.0.1'
def go(args):
if sys.version_info < (3,0):
print('Boiler requires Python version 3 or better to run')
exit()
if args.command == 'compress':
import compress
if args.preprocess:
modes = ['hisat']
if args.preprocess.lower() == 'hisat':
print('Preprocessing HISAT alignments')
import enumeratePairs
prefix = args.alignments[:args.alignments.index('.')] + '.processed'
enumeratePairs.processHISAT(args.alignments, prefix + '.sam')
os.system('samtools view -bS ' + prefix + '.sam | samtools sort - ' + prefix)
os.system('samtools view -h -o ' + prefix + '.sam ' + prefix + '.bam')
else:
print('Preprocessing mode not recognized: %s' % args.preprocess)
print('Supported preprocessing modes include: ' + ', '.join(modes))
exit()
args.alignments = prefix + '.sam'
if args.verbose:
print('Compressing')
start = time.time()
compressor = compress.Compressor(args.frag_len_cutoff)
compressor.compress(args.alignments, args.compressed, args.gtf, None, args.frag_len_z_cutoff, args.split_diff_strands, args.split_discordant)
if args.verbose:
end = time.time()
print('Compression took %0.3f s' % (end-start))
elif args.command == 'query':
import expand
expander = expand.Expander()
if args.counts:
if not args.gtf:
print('Missing required argument for counts query: --gtf path/to/reference.gtf')
exit()
if not args.output:
print('Counts query cannot write to standard output -- please include a file prefix for output')
exit()
exon_counts, junc_counts = expander.getCounts(args.compressed, args.gtf)
expander.write_counts(exon_counts, junc_counts, args.output)
else:
if not args.chrom:
print('Missing required argument for query: --chrom chromosome')
exit()
if args.output:
#logging.info('Opening %s to write results' % args.output)
try:
f = open(args.output, 'w')
except IOError:
#logging.info('Couldn\'t open file %s for writing. Using standard out instead.' % args.output)
f = sys.stdout
else:
#logging.warning('Writing results to standard out')
f = sys.stdout
if args.bundles:
bundles = expander.getGeneBounds(args.compressed, args.chrom, args.start, args.end)
for b in bundles:
f.write(str(b[0])+'\t'+str(b[1])+'\n')
if args.coverage:
cov = expander.getCoverage(args.compressed, args.chrom, args.start, args.end)
f.write(','.join([str(c) for c in cov]) + '\n')
if args.reads:
aligned, unpaired, paired = expander.getReads(args.compressed, args.chrom, args.start, args.end)
aligned.writeSAM(f, unpaired, paired, False, False, 0)
elif args.command == 'decompress':
import expand
if args.verbose:
print('Decompressing')
start = time.time()
expander = expand.Expander(args.force_xs)
expander.expand(args.compressed, args.expanded)
if args.verbose:
end = time.time()
print('Decompression took %0.3f s' % (end-start))
if __name__ == '__main__':
if '--version' in sys.argv:
print('Boiler v' + VERSION)
sys.exit(0)
# Check for Python version
version = sys.version_info[0]
if version < 3:
print('Python 3 is required to run Boiler')
exit()
# Print file's docstring if -h is invoked
parser = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter, prog='Boiler')
subparsers = parser.add_subparsers(help='Commands', dest='command')
parser_compress = subparsers.add_parser('compress', help="Compress a SAM file")
parser_compress.add_argument("-c", "--frag-len-cutoff", type=int, help='Store any fragments longer than this in a bundle-spanning bucket')
parser_compress.add_argument("-z", "--frag-len-z-cutoff", type=float, help='Store any fragments above this z-score in a bundle-spanning bucket')
parser_compress.add_argument("-s", "--split-diff-strands", action="store_true", help='Split any pairs with different XS values')
parser_compress.add_argument("-d", "--split-discordant", action="store_true", help='Treat discordant pairs as unpaired reads')
parser_compress.add_argument("-p", "--preprocess", type=str, help="Set to 'tophat' to preprocess TopHat alignments, 'hisat' to preprocess HISAT alignments")
parser_compress.add_argument("-g", "--gtf", type=str, help="Path to reference GTF to improve compression accuracy (this will result in larger file size)")
parser_compress.add_argument("-v", "--verbose", help="Print timing information", action="store_true")
parser_compress.add_argument("alignments", type=str, help='Full path of SAM file containing aligned reads')
parser_compress.add_argument("compressed", type=str, nargs='?', default='compressed.bin', help="Compressed filename. Default: compressed.bin")
parser_query = subparsers.add_parser('query', help="Query compressed file")
group = parser_query.add_mutually_exclusive_group()
group.add_argument('-b', '--bundles', help="Query bundles", action="store_true")
group.add_argument('-c', '--coverage', help="Query coverage", action="store_true")
group.add_argument('-r', '--reads', help="Query reads", action="store_true")
group.add_argument('-f', '--counts', help='Query read counts over exons and splice junctions in a GTF', action="store_true")
parser_query.add_argument('--chrom', help="Chromosome to query", type=str)
parser_query.add_argument('--start', help="Beginning of range to query", type=int)
parser_query.add_argument('--end', help="End of range to query", type=int)
parser_query.add_argument('--gtf', help="Path to reference GTF over which to query counts")
parser_query.add_argument('compressed', help="Path to compressed file created by Boiler", type=str)
parser_query.add_argument('output', nargs='?', default=None, help="File to write result to. Default: Standard out")
parser_decompress = subparsers.add_parser('decompress', help="Decompress to a SAM file")
parser_decompress.add_argument("-f", "--force-xs", help="If we decompress a spliced read with no XS value, assign it a random one (so Cufflinks can run)", action="store_true")
parser_decompress.add_argument("-v", "--verbose", help="Print timing information", action="store_true")
parser_decompress.add_argument("compressed", type=str, help="Compressed filename")
parser_decompress.add_argument("expanded", type=str, nargs='?', default='expanded.sam', help="Write decompressed SAM to this filename. Default: expanded.sam")
args = parser.parse_args(sys.argv[1:])
go(args)