Mercurial > repos > devteam > sam2interval
comparison sam2interval.py @ 0:8c737b8ddc45 draft
Uploaded tool tarball.
| author | devteam |
|---|---|
| date | Mon, 26 Aug 2013 15:12:38 -0400 |
| parents | |
| children | 75557c0908a9 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:8c737b8ddc45 |
|---|---|
| 1 #!/usr/bin/env python | |
| 2 | |
| 3 import sys | |
| 4 import optparse | |
| 5 import re | |
| 6 | |
| 7 def stop_err( msg ): | |
| 8 sys.stderr.write( msg ) | |
| 9 sys.exit() | |
| 10 | |
| 11 def main(): | |
| 12 usage = """%prog [options] | |
| 13 | |
| 14 options (listed below) default to 'None' if omitted | |
| 15 """ | |
| 16 parser = optparse.OptionParser(usage=usage) | |
| 17 | |
| 18 parser.add_option( | |
| 19 '-f','--input_sam_file', | |
| 20 metavar="INPUT_SAM_FILE", | |
| 21 dest='input_sam', | |
| 22 default = False, | |
| 23 help='Name of the SAM file to be filtered. STDIN is default') | |
| 24 | |
| 25 parser.add_option( | |
| 26 '-c','--flag_column', | |
| 27 dest='flag_col', | |
| 28 default = '2', | |
| 29 help='Column containing SAM bitwise flag. 1-based') | |
| 30 | |
| 31 parser.add_option( | |
| 32 '-s','--start_column', | |
| 33 dest='start_col', | |
| 34 default = '4', | |
| 35 help='Column containing position. 1-based') | |
| 36 | |
| 37 parser.add_option( | |
| 38 '-g','--cigar_column', | |
| 39 dest='cigar_col', | |
| 40 default = '6', | |
| 41 help='Column containing CIGAR or extended CIGAR string') | |
| 42 | |
| 43 parser.add_option( | |
| 44 '-r','--ref_column', | |
| 45 dest='ref_col', | |
| 46 default = '3', | |
| 47 help='Column containing name of the reference sequence coordinate. 1-based') | |
| 48 | |
| 49 parser.add_option( | |
| 50 '-e','--read_column', | |
| 51 dest='read_col', | |
| 52 default = '1', | |
| 53 help='Column containing read name. 1-based') | |
| 54 | |
| 55 parser.add_option( | |
| 56 '-p','--print_all', | |
| 57 dest='prt_all', | |
| 58 action='store_true', | |
| 59 default = False, | |
| 60 help='Print coordinates and original SAM?') | |
| 61 | |
| 62 options, args = parser.parse_args() | |
| 63 | |
| 64 if options.input_sam: | |
| 65 infile = open ( options.input_sam, 'r') | |
| 66 else: | |
| 67 infile = sys.stdin | |
| 68 | |
| 69 cigar = re.compile( '\d+M|\d+N|\d+D|\d+P' ) | |
| 70 | |
| 71 print '#chrom\tstart\tend\tstrand\tread_name' # provide a (partial) header so that strand is automatically set in metadata | |
| 72 | |
| 73 for line in infile: | |
| 74 line = line.rstrip( '\r\n' ) | |
| 75 if line and not line.startswith( '#' ) and not line.startswith( '@' ) : | |
| 76 fields = line.split( '\t' ) | |
| 77 start = int( fields[ int( options.start_col ) - 1 ] ) - 1 | |
| 78 end = 0 | |
| 79 for op in cigar.findall( fields[ int( options.cigar_col) - 1 ] ): | |
| 80 end += int( op[ 0:len( op ) - 1 ] ) | |
| 81 | |
| 82 strand = '+' | |
| 83 if bool( int( fields[ int( options.flag_col ) - 1 ] ) & 0x0010 ): | |
| 84 strand = '-' | |
| 85 read_name = fields[ int( options.read_col ) - 1 ] | |
| 86 ref_name = fields[ int( options.ref_col ) - 1 ] | |
| 87 | |
| 88 if ref_name != '*': | |
| 89 # Do not print lines with unmapped reads that contain '*' instead of chromosome name | |
| 90 if options.prt_all: | |
| 91 print '%s\t%s\t%s\t%s\t%s' % (ref_name, str(start), str(end+start), strand, line) | |
| 92 else: | |
| 93 print '%s\t%s\t%s\t%s\t%s' % (ref_name, str(start), str(end+start), strand, read_name) | |
| 94 | |
| 95 if __name__ == "__main__": main() | |
| 96 |