-
Notifications
You must be signed in to change notification settings - Fork 2
/
extract_junc.py
executable file
·61 lines (53 loc) · 2 KB
/
extract_junc.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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
Usage: extract_junc.py [options] <bam>
Options:
-h --help Show help message.
--version Show version.
-o output_dir Output file directory [default: ./].
--url Extract from remote url.
--bb Convert to BigWig.
--uniq Only fetch unique mapped reads.
--stranded Retain strand information.
--min-reads=min Minimum junction reads [default: 0].
'''
import sys
import os
import os.path
import tempfile
import pysam
from docopt import docopt
from seqlib.ngs import fetch_juncfile
from seqlib.path import which, check_dir
from seqlib.version import __version__
__author__ = 'Xiao-Ou Zhang <[email protected]>'
def main():
# parse options
options = docopt(__doc__, version=__version__)
# check output_dir
if options['-o'] == './':
dir = os.getcwd()
else:
dir = check_dir(options['-o'])
# fetch junction bed file
junc_f = fetch_juncfile(options['<bam>'], url=options['--url'], dir=dir,
uniq=options['--uniq'],
stranded=options['--stranded'],
min=int(options['--min-reads']))
# create junction bigbed file in case
if options['--bb'] and which('bedToBigBed') is not None:
prefix = os.path.splitext(os.path.split(options['<bam>'])[-1])[0]
bamf = pysam.AlignmentFile(options['<bam>'], 'rb')
with tempfile.NamedTemporaryFile() as chrom_size:
for seq in bamf.header['SQ']:
chrom_size.write('%s\t%s\n' % (seq['SN'], seq['LN']))
chrom_size.seek(0)
bb_path = os.path.join(dir, prefix + '_junc.bb')
return_code = os.system('bedToBigBed -type=bed12 %s %s %s' %
(junc_f, chrom_size.name, bb_path)) >> 8
if return_code:
sys.exit('Error: cannot convert bed to BigBed!')
bamf.close()
if __name__ == '__main__':
main()