<s> import ngsutils . bam <EOL> import ngsutils . bam . export <EOL> import ngsutils . bam . basecall <EOL> import ngsutils . bed <EOL> import StringIO <EOL> import collections <EOL> PileupRecords = collections . namedtuple ( '<STR_LIT>' , '<STR_LIT>' ) <EOL> PileupRead = collections . namedtuple ( '<STR_LIT>' , '<STR_LIT>' ) <EOL> def assertIn ( foo , bar ) : <EOL> for el in bar : <EOL> if foo == el : <EOL> return True <EOL> return False <EOL> def assertNotIn ( foo , bar ) : <EOL> for el in bar : <EOL> if foo == el : <EOL> return False <EOL> return True <EOL> def _matches ( valid , queries ) : <EOL> extra = False <EOL> check = [ False , ] * len ( valid ) <EOL> for q in queries : <EOL> found = False <EOL> for i , v in enumerate ( valid ) : <EOL> if v == q : <EOL> check [ i ] = True <EOL> found = True <EOL> break <EOL> if not found : <EOL> extra = True <EOL> if extra : <EOL> return False <EOL> for c in check : <EOL> if not c : <EOL> return False <EOL> return True <EOL> class MockBam ( object ) : <EOL> def __init__ ( self , refs , lengths = None , insert_order = False ) : <EOL> self . _refs = refs [ : ] <EOL> self . _ref_lengths = [ ] <EOL> if lengths : <EOL> for length in lengths : <EOL> self . _ref_lengths . append ( length ) <EOL> else : <EOL> for ref in refs : <EOL> self . _ref_lengths . append ( <NUM_LIT:0> ) <EOL> self . _reads = { } <EOL> self . _read_keys = [ ] <EOL> self . filename = '<STR_LIT>' <EOL> self . __iter = None <EOL> self . __pos = <NUM_LIT:0> <EOL> self . _insert_order = insert_order <EOL> self . _read_insert_order = [ ] <EOL> def seek ( self , pos ) : <EOL> self . __pos = pos <EOL> def tell ( self ) : <EOL> return self . __pos <EOL> def write ( self , read ) : <EOL> if not self . _insert_order : <EOL> k = ( read . tid , read . pos , len ( self . _reads ) ) <EOL> self . _reads [ k ] = read <EOL> self . _read_keys . append ( k ) <EOL> self . _read_keys . sort ( ) <EOL> else : <EOL> self . _reads [ read . qname ] = read <EOL> self . _read_keys . append ( read . qname ) <EOL> if self . _ref_lengths [ read . tid ] < read . aend : <EOL> self . _ref_lengths [ read . tid ] = read . aend <EOL> def add_read ( self , read , * args , ** kwargs ) : <EOL> if type ( read ) != MockRead : <EOL> read = MockRead ( read , * args , ** kwargs ) <EOL> self . write ( read ) <EOL> self . _read_insert_order . append ( read . qname ) <EOL> return self <EOL> def __iter__ ( self ) : <EOL> self . __iter = self . fetch ( ) <EOL> return self <EOL> def next ( self ) : <EOL> try : <EOL> return self . __iter . next ( ) <EOL> except : <EOL> self . __iter = None <EOL> raise StopIteration <EOL> def pileup ( self , ref , start , end ) : <EOL> '''<STR_LIT>''' <EOL> basecaller = ngsutils . bam . basecall . BamBaseCaller ( self , regions = ngsutils . bed . BedFile ( region = '<STR_LIT>' % ( ref , start + <NUM_LIT:1> , end ) ) ) <EOL> for basepos in basecaller . fetch ( ) : <EOL> pileups = [ ] <EOL> for record in basepos . reads : <EOL> indel = <NUM_LIT:0> <EOL> level = <NUM_LIT:0> <EOL> is_del = False <EOL> is_head = False <EOL> is_tail = False <EOL> if record . qpos == <NUM_LIT:0> : <EOL> if record . read . is_reverse : <EOL> is_tail = True <EOL> else : <EOL> is_head = True <EOL> if basepos . insertions : <EOL> indel = max ( [ len ( ins ) for ins in basepos . insertions ] ) <EOL> if basepos . deletions : <EOL> indel = - max ( basepos . deletions ) <EOL> is_del = True <EOL> pr = PileupRead ( record . read , indel , is_del , is_head , is_tail , level , record . qpos ) <EOL> pileups . append ( pr ) <EOL> yield PileupRecords ( basepos . tid , basepos . pos , basepos . total , pileups ) <EOL> def fetch ( self , ref = None , start = - <NUM_LIT:1> , end = - <NUM_LIT:1> , region = None ) : <EOL> if region : <EOL> print "<STR_LIT>" <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> self . __pos = <NUM_LIT:0> <EOL> tid = - <NUM_LIT:1> <EOL> if ref and start > - <NUM_LIT:1> and end > - <NUM_LIT:1> : <EOL> for i , val in enumerate ( self . _refs ) : <EOL> if val == ref : <EOL> tid = i <EOL> break <EOL> if tid == - <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" % ref ) <EOL> for k in self . _read_keys : <EOL> self . __pos += <NUM_LIT:1> <EOL> read = self . _reads [ k ] <EOL> if tid == - <NUM_LIT:1> : <EOL> yield read <EOL> else : <EOL> if self . _reads [ k ] . tid == tid : <EOL> if start <= read . pos <= end or start <= read . aend <= end or ( read . pos < start and read . aend > end ) : <EOL> yield read <EOL> @ property <EOL> def references ( self ) : <EOL> return self . _refs [ : ] <EOL> @ property <EOL> def lengths ( self ) : <EOL> return self . _ref_lengths [ : ] <EOL> def getrname ( self , tid ) : <EOL> return self . _refs [ tid ] <EOL> def close ( self ) : <EOL> pass <EOL> class MockRead ( object ) : <EOL> def __init__ ( self , qname , seq = None , qual = None , tid = - <NUM_LIT:1> , pos = - <NUM_LIT:1> , aend = None , cigar = None , tags = None , mapq = <NUM_LIT:0> , rnext = - <NUM_LIT:1> , pnext = - <NUM_LIT:1> , isize = <NUM_LIT:0> , tlen = <NUM_LIT:0> , is_reverse = False , is_paired = False , is_secondary = False , is_qcfail = False , mate_is_unmapped = False , mate_is_reverse = False , is_read1 = False , is_read2 = False , flag = <NUM_LIT:0> , bam = None ) : <EOL> self . qname = qname <EOL> self . seq = seq . upper ( ) if seq else '<STR_LIT>' <EOL> self . qual = qual <EOL> self . tid = tid <EOL> self . pos = pos <EOL> self . aend = aend <EOL> self . mapq = mapq <EOL> self . rnext = rnext <EOL> self . pnext = pnext <EOL> self . isize = isize <EOL> self . tlen = tlen <EOL> self . bam = bam <EOL> self . is_paired = is_paired <EOL> self . is_reverse = is_reverse <EOL> self . is_secondary = is_secondary <EOL> self . is_qcfail = is_qcfail <EOL> self . mate_is_reverse = mate_is_reverse <EOL> self . is_read1 = is_read1 <EOL> self . is_read2 = is_read2 <EOL> if flag : <EOL> self . _flag = flag <EOL> self . is_paired = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_unmapped = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . mate_is_unmapped = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_reverse = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . mate_is_reverse = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_read1 = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_read2 = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_secondary = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> self . is_qcfail = <NUM_LIT> & self . flag > <NUM_LIT:0> <EOL> else : <EOL> self . _flag = <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_paired else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_unmapped else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . mate_is_unmapped else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_reverse else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . mate_is_reverse else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_read1 else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_read2 else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_secondary else <NUM_LIT:0> <EOL> self . _flag |= <NUM_LIT> if self . is_qcfail else <NUM_LIT:0> <EOL> if self . tid == - <NUM_LIT:1> or self . pos == - <NUM_LIT:1> : <EOL> self . is_unmapped = True <EOL> if self . is_paired : <EOL> if self . rnext == - <NUM_LIT:1> or self . pnext == - <NUM_LIT:1> : <EOL> self . mate_is_unmapped = True <EOL> if cigar : <EOL> self . cigar = ngsutils . bam . cigar_fromstr ( cigar ) <EOL> if not aend : <EOL> self . aend = self . pos + ngsutils . bam . cigar_read_len ( self . cigar ) <EOL> else : <EOL> self . cigar = None <EOL> if tags : <EOL> self . tags = tags <EOL> else : <EOL> self . tags = [ ( '<STR_LIT>' , <NUM_LIT:1> ) , ( '<STR_LIT>' , <NUM_LIT:1> ) ] <EOL> def opt ( self , k ) : <EOL> for k1 , val in self . tags : <EOL> if k1 == k : <EOL> return val <EOL> return None <EOL> @ property <EOL> def flag ( self ) : <EOL> return self . _flag <EOL> @ property <EOL> def is_unmapped ( self ) : <EOL> if self . tid > - <NUM_LIT:1> and self . pos > - <NUM_LIT:1> : <EOL> return False <EOL> return True <EOL> @ is_unmapped . setter <EOL> def is_unmapped ( self , val ) : <EOL> if val : <EOL> self . tid = - <NUM_LIT:1> <EOL> self . pos = - <NUM_LIT:1> <EOL> self . _flag |= <NUM_LIT> <EOL> else : <EOL> self . tid = <NUM_LIT:0> <EOL> self . pos = <NUM_LIT:0> <EOL> self . _flag = self . _flag & <NUM_LIT> <EOL> @ property <EOL> def mate_is_unmapped ( self ) : <EOL> if not self . is_paired : <EOL> return False <EOL> if self . rnext > - <NUM_LIT:1> and self . pnext > - <NUM_LIT:1> : <EOL> return False <EOL> return True <EOL> @ mate_is_unmapped . setter <EOL> def mate_is_unmapped ( self , val ) : <EOL> if val : <EOL> self . rnext = - <NUM_LIT:1> <EOL> self . pnext = - <NUM_LIT:1> <EOL> self . _flag |= <NUM_LIT> <EOL> else : <EOL> self . rnext = <NUM_LIT:0> <EOL> self . pnext = <NUM_LIT:0> <EOL> self . _flag = self . _flag & <NUM_LIT> <EOL> def __repr__ ( self ) : <EOL> out = StringIO . StringIO ( '<STR_LIT>' ) <EOL> ngsutils . bam . export . export_read ( self . bam , self , [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] , out = out ) <EOL> return out . getvalue ( ) . rstrip ( ) </s>
