Mercurial > repos > iuc > ngsutils_bam_filter

comparison ngsutils/bam/__init__.py @ 2:7a68005de299 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ngsutils commit 9a243c616a4a3156347e38fdb5f35863ae5133f9
author iuc
date Sun, 27 Nov 2016 15:01:21 -0500
parents 4e4e4093d65d
children
comparison
equal deleted inserted replaced
1:8187a729d9f4 2:7a68005de299
1 import sys
2 import os 1 import os
3 import re 2 import re
3 import sys
4
5 import ngsutils.support
4 import pysam 6 import pysam
5 try: 7 try:
6 from eta import ETA 8 from eta import ETA
7 except: 9 except:
8 pass 10 pass
9 import ngsutils.support
10 11
11 12
12 def bam_open(fname, mode='r', *args, **kwargs): 13 def bam_open(fname, mode='r', *args, **kwargs):
13 if fname.lower()[-4:] == '.bam': 14 if fname.lower()[-4:] == '.bam':
14 return pysam.Samfile(fname, '%sb' % mode, *args, **kwargs) 15 return pysam.Samfile(fname, '%sb' % mode, *args, **kwargs)
157 >>> cigar_tostr(((0, 10), (1, 5), (0, 20))) 158 >>> cigar_tostr(((0, 10), (1, 5), (0, 20)))
158 '10M5I20M' 159 '10M5I20M'
159 >>> cigar_tostr(((0, 1), (1, 1), (2, 1), (3, 1), (4, 1), (5, 1), (6, 1))) 160 >>> cigar_tostr(((0, 1), (1, 1), (2, 1), (3, 1), (4, 1), (5, 1), (6, 1)))
160 '1M1I1D1N1S1H1P' 161 '1M1I1D1N1S1H1P'
161 ''' 162 '''
162 163
163 s = '' 164 s = ''
164 165
165 for op, size in cigar: 166 for op, size in cigar:
166 s += '%s%s' % (size, bam_cigar[op]) 167 s += '%s%s' % (size, bam_cigar[op])
167 168
228 229
229 def _extract_md_matches(md, maxlength): 230 def _extract_md_matches(md, maxlength):
230 md_pos = 0 231 md_pos = 0
231 232
232 while md and md_pos < maxlength: 233 while md and md_pos < maxlength:
233 tmp = '0' # preload a zero so that immediate mismatches will be caught 234 # preload a zero so that immediate mismatches will be caught
234 # the zero will have no affect otherwise... 235 # the zero will have no affect otherwise...
236 tmp = '0'
235 237
236 # look for matches 238 # look for matches
237 while md and md[0] in '0123456789': 239 while md and md[0] in '0123456789':
238 tmp += md[0] 240 tmp += md[0]
239 md = md[1:] 241 md = md[1:]
623 chr_cigar.append((op, frag_end - cur_pos)) 625 chr_cigar.append((op, frag_end - cur_pos))
624 length -= (frag_end - cur_pos) 626 length -= (frag_end - cur_pos)
625 cur_pos = frag_end 627 cur_pos = frag_end
626 frag_idx += 1 628 frag_idx += 1
627 if len(fragments) <= frag_idx: 629 if len(fragments) <= frag_idx:
628 print 'ERROR converting: ', name, fragments 630 print 'ERROR converting: ', name, fragments
629 return (chrom, 0, chr_cigar) 631 return (chrom, 0, chr_cigar)
630 frag_start, frag_end = fragments[frag_idx] 632 frag_start, frag_end = fragments[frag_idx]
631 chr_cigar.append((3, frag_start - cur_pos)) 633 chr_cigar.append((3, frag_start - cur_pos))
632 cur_pos = frag_start 634 cur_pos = frag_start
633 635
862 newread.is_paired = True 864 newread.is_paired = True
863 865
864 if not read.is_unmapped and read.is_reverse: 866 if not read.is_unmapped and read.is_reverse:
865 newread.seq = ngsutils.support.revcomp(read.seq) 867 newread.seq = ngsutils.support.revcomp(read.seq)
866 newread.qual = read.qual[::-1] 868 newread.qual = read.qual[::-1]
867 else: 869 else:
868 newread.seq = read.seq 870 newread.seq = read.seq
869 newread.qual = read.qual 871 newread.qual = read.qual
870 872
871 newread.tags = tags 873 newread.tags = tags
872 874
873 return newread 875 return newread
874
875 876
876 877
877 if __name__ == '__main__': 878 if __name__ == '__main__':
878 import doctest 879 import doctest
879 doctest.testmod() 880 doctest.testmod()