-
Notifications
You must be signed in to change notification settings - Fork 11
/
Copy pathbedtools.py
executable file
·93 lines (78 loc) · 3.4 KB
/
bedtools.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
#!/usr/bin/env python
from optparse import OptionParser
import os, subprocess, tempfile
import pysam
################################################################################
# bedtools.py
#
#
################################################################################
################################################################################
# main
################################################################################
def main():
usage = 'usage: %prog [options] arg'
parser = OptionParser(usage)
#parser.add_option()
(options,args) = parser.parse_args()
################################################################################
# abam_f1
#
# Intersect the BAM file with the BED file using the "-f 1" option, but correct
# for the loss of spliced reads.
################################################################################
def abam_f1(bam_file, bed_file, out_file):
############################################
# divide BAM by splicing
############################################
spliced_bam_fd, spliced_bam_file = tempfile.mkstemp(dir='%s/research/scratch/temp' % os.environ['HOME'])
unspliced_bam_fd, unspliced_bam_file = tempfile.mkstemp(dir='%s/research/scratch/temp' % os.environ['HOME'])
# open BAMs
bam_in = pysam.Samfile(bam_file, 'rb')
spliced_bam_out = pysam.Samfile(spliced_bam_file, 'wb', template=bam_in)
unspliced_bam_out = pysam.Samfile(unspliced_bam_file, 'wb', template=bam_in)
# divide
for aligned_read in bam_in:
if spliced(aligned_read):
spliced_bam_out.write(aligned_read)
else:
unspliced_bam_out.write(aligned_read)
# close BAMs
bam_in.close()
spliced_bam_out.close()
unspliced_bam_out.close()
############################################
# intersect and merge
############################################
spliced_is_bam_fd, spliced_is_bam_file = tempfile.mkstemp(dir='%s/research/scratch/temp' % os.environ['HOME'])
unspliced_is_bam_fd, unspliced_is_bam_file = tempfile.mkstemp(dir='%s/research/scratch/temp' % os.environ['HOME'])
subprocess.call('intersectBed -f 1 -abam %s -b %s > %s' % (unspliced_bam_file, bed_file, unspliced_is_bam_file), shell=True)
subprocess.call('intersectBed -abam %s -b %s > %s' % (spliced_bam_file, bed_file, spliced_is_bam_file), shell=True)
subprocess.call('samtools merge -f %s %s %s' % (out_file, unspliced_is_bam_file, spliced_is_bam_file), shell=True)
############################################
# clean
############################################
os.close(spliced_bam_fd)
os.remove(spliced_bam_file)
os.close(unspliced_bam_fd)
os.remove(unspliced_bam_file)
os.close(spliced_is_bam_fd)
os.remove(spliced_is_bam_file)
os.close(unspliced_is_bam_fd)
os.remove(unspliced_is_bam_file)
################################################################################
# spliced
#
# Return true if the read is spliced.
################################################################################
def spliced(aligned_read):
spliced = False
for code,size in aligned_read.cigar:
if code == 3:
spliced = True
return spliced
################################################################################
# __main__
################################################################################
if __name__ == '__main__':
main()