For a given CIGAR string, return the start positions of
        each aligned segment in ref and tgt, and a list of CIGAR operators.
321263
def _parse_cigar(self, cigar):
        
        ces = [m.groupdict() for m in cigar_re.finditer(cigar)]
        ref_pos = [None] * len(ces)
        tgt_pos = [None] * len(ces)
        cigar_op = [None] * len(ces)
        ref_cur = tgt_cur = 0
        for i, ce in enumerate(ces):
            ref_pos[i] = ref_cur
            tgt_pos[i] = tgt_cur
            cigar_op[i] = ce["op"]
            step = int(ce["len"])
            if ce["op"] in "=MINX":
                ref_cur += step
            if ce["op"] in "=MDX":
                tgt_cur += step
        ref_pos.append(ref_cur)
        tgt_pos.append(tgt_cur)
        return ref_pos, tgt_pos, cigar_op
1167468
def find_splice(cigar):
    ''''''

    offset = 0
    # a soft clip at the end of the read is taken as splicing
    # where as a soft clip at the start is not.
    if cigar[0][0] == 4:
        offset = cigar[0][1]
        cigar = cigar[1:]

    for op, bases in cigar:
        if op in (3, 4):
            # N or S: found the splice
            return offset
        elif op in (0, 2, 7, 8):
            # M, D, = or X: reference consuming
            offset += bases
        elif op in (1, 5, 6):
            # I, H, P: non-reference consuming
            continue
        else:
            raise ValueError("Bad Cigar operation: %i" % op)

    return False
157699
def _segment_with_tokens(text, tokens):
        """"""
        list_form = []
        text_ptr = 0
        for token in tokens:
            inter_token_string = []
            while not text[text_ptr:].startswith(token):
                inter_token_string.append(text[text_ptr])
                text_ptr += 1
                if text_ptr >= len(text):
                    raise ValueError("Tokenization produced tokens that do not belong in string!")
            text_ptr += len(token)
            if inter_token_string:
                list_form.append(''.join(inter_token_string))
            list_form.append(token)
        if text_ptr < len(text):
            list_form.append(text[text_ptr:])
        return list_form
810448
def beginning_of_code_block(node, line_number, full_contents, is_code_node):
    """"""
    if SPHINX_INSTALLED and not is_code_node:
        delta = len(node.non_default_attributes())
        current_line_contents = full_contents.splitlines()[line_number:]
        blank_lines = next(
            (i for (i, x) in enumerate(current_line_contents) if x),
            0)
        return (
            line_number +
            delta - 1 +
            blank_lines - 1 +
            SPHINX_CODE_BLOCK_DELTA)
    else:
        lines = full_contents.splitlines()
        code_block_length = len(node.rawsource.splitlines())

        try:
            # Case where there are no extra spaces.
            if lines[line_number - 1].strip():
                return line_number - code_block_length + 1
        except IndexError:
            pass

        # The offsets are wrong if the RST text has multiple blank lines after
        # the code block. This is a workaround.
        for line_number in range(line_number, 1, -1):
            if lines[line_number - 2].strip():
                break

        return line_number - code_block_length
210950
def extract_finditer(pos_seq, regex=SimpleNP):
	""""""
	ss = coarse_tag_str(pos_seq)

	def gen():
		for m in re.finditer(regex, ss):
			yield (m.start(), m.end())

	return list(gen())
1018727
def getStartNodes(fdefs,calls):
    ''''''
    s=[]
    for source in fdefs:
        for fn in fdefs[source]:
            inboundEdges=False
            for call in calls:
                if call.target==fn:
                    inboundEdges=True
            if not inboundEdges:
                s.append(fn)
    return s
340103
public static List<ReadBaseWithReference> extractReadBases(Read read) {

    // Make sure this read has a valid alignment with Cigar Units
    if (!read.hasAlignment() || (read.getAlignment().getCigarCount() == 0)) {
      return null;
    }

    ImmutableList.Builder<ReadBaseWithReference> bases = ImmutableList.builder();

    String readSeq = read.getAlignedSequence();
    List<Integer> readQual = read.getAlignedQualityList();
    String refSeq = UNINITIALIZED_REFERENCE_SEQUENCE;

    int refPosAbsoluteOffset = 0;
    int readOffset = 0;
    for (CigarUnit unit : read.getAlignment().getCigarList()) {
      switch(unit.getOperation()) {
        case ALIGNMENT_MATCH:
        case SEQUENCE_MISMATCH:
        case SEQUENCE_MATCH:
          for (int i = 0; i < unit.getOperationLength(); i++) {
            String refBase = "";
            if (unit.getOperation().equals(CigarUnit.Operation.SEQUENCE_MATCH)) {
              refBase = readSeq.substring(readOffset, readOffset + 1);
            } else if (!unit.getReferenceSequence().isEmpty()) {
              // try to get the ref sequence from the Cigar unit
              refBase = unit.getReferenceSequence().substring(i, i + 1);
            } else {
              // try to get the ref sequence by fully parsing the MD tag if not already cached
              if (refSeq != null && refSeq.equals(UNINITIALIZED_REFERENCE_SEQUENCE)) {
                refSeq = com.google.cloud.genomics.utils.grpc.ReadUtils
                    .inferReferenceSequenceByParsingMdFlag(read);
              }
              if (refSeq != null) {
                refBase = refSeq.substring(readOffset, readOffset + 1);
              }
            }
            String name = read.getAlignment().getPosition().getReferenceName();
            Matcher m = Pattern.compile("^(chr)?(X|Y|([12]?\\d))$").matcher(name);
            if (m.matches()) {
              name = m.group(m.groupCount() - 1);
            }
            Position refPosition = Position.newBuilder()
                .setReferenceName(name)
                .setPosition(read.getAlignment().getPosition().getPosition()
                  + refPosAbsoluteOffset)
                  .build();
            bases.add(new ReadBaseWithReference(new ReadBaseQuality(
                readSeq.substring(readOffset, readOffset + 1),
                readQual.get(readOffset)),
                refBase,
                refPosition));
            refPosAbsoluteOffset++;
            readOffset++;
          }
          break;

        case PAD:        // padding (silent deletion from padded reference)
        case CLIP_HARD:  // hard clipping (clipped sequences NOT present in seq)
          break;

        case CLIP_SOFT:  // soft clipping (clipped sequences present in SEQ)
        case INSERT:     // insertion to the reference
          readOffset += unit.getOperationLength();
          break;

        case DELETE:  // deletion from the reference
        case SKIP:    // intron (mRNA-to-genome alignment only)
          refPosAbsoluteOffset += unit.getOperationLength();
          break;

        default:
          throw new IllegalArgumentException("Illegal cigar code: " + unit.getOperation());
      }
    }

    return bases.build();
  }
297588
def _get_alignment_lines(self):
        ''' '''

        # Strip the boilerplate lines
        lines = self.alignment_output.split("\n")
        assert(lines[0].startswith('CLUSTAL'))
        lines = '\n'.join(lines[1:]).lstrip().split('\n')

        # The sequence IDs should be unique. Reassert this here
        assert(len(self.sequence_ids.values()) == len(set(self.sequence_ids.values())))

        # Create the list of sequence IDs
        id_list = [v for k, v in sorted(self.sequence_ids.iteritems())]

        # Determine the indentation level
        first_id = id_list[0]
        header_regex = re.compile("(.*?\s+)(.*)")
        alignment_regex = re.compile("^([A-Z\-]+)\s*$")
        mtchs = header_regex.match(lines[0])
        assert(mtchs.group(1).strip() == first_id)
        indentation = len(mtchs.group(1))
        sequence = mtchs.group(2)
        assert(sequence)
        assert(alignment_regex.match(sequence))

        # Create empty lists for the sequences
        sequences = {}
        for id in id_list:
            sequences[id] = []
        sequences[-1] = []

        # Get the lists of sequences
        num_ids = len(id_list)
        for x in range(0, len(lines), num_ids + 2):
            for y in range(num_ids):
                id = id_list[y]
                assert(lines[x + y][:indentation].strip() == id)
                assert(lines[x + y][indentation - 1] == ' ')
                sequence = lines[x + y][indentation:].strip()
                assert(alignment_regex.match(sequence))
                sequences[id].append(sequence)

            # Get the length of the sequence lines
            length_of_sequences = list(set(map(len, [v[-1] for k, v in sequences.iteritems() if k != -1])))
            assert(len(length_of_sequences) == 1)
            length_of_sequences = length_of_sequences[0]

            # Parse the Clustal match line
            assert(lines[x + num_ids][:indentation].strip() == '')
            match_sequence = lines[x + num_ids][indentation:indentation + length_of_sequences]
            assert(match_sequence.strip() == lines[x + num_ids].strip())
            assert(lines[x + y][indentation - 1] == ' ')
            sequences[-1].append(match_sequence)

            # Check for the empty line
            assert(lines[x + num_ids + 1].strip() == '')

        # Create the sequences, making sure that all sequences are the same length
        lengths = set()
        for k, v in sequences.iteritems():
            sequences[k] = "".join(v)
            lengths.add(len(sequences[k]))
        assert(len(lengths) == 1)

        return sequences
1291735
private static List<Segmentation> getSegmentations(String str, int begin, int end, int nbMaxComponents, int minComponentSize, List<CuttingPoint> cuttingPoints) {
		List<Segmentation> s = new ArrayList<Segmentation>(nbMaxComponents);
		for(CuttingPoint cp:getPossibleCuttingPoints(str, begin, end, nbMaxComponents, minComponentSize)) {
			s.add(new Segmentation(str, addToDuplicate(cuttingPoints, cp, nbMaxComponents - 1)));
			if(nbMaxComponents > 2) {
				if(cp.isHypen())
					// recursivity on the left part only if the cp was an hypen
					if((cp.getIndex() - begin) >= 2*minComponentSize)
						s.addAll(getSegmentations(str, begin, cp.getIndex(), nbMaxComponents -1, minComponentSize, addToDuplicate(cuttingPoints, cp, nbMaxComponents - 1)));
				
				if((end - (cp.getIndex() + cp.getOffset())) >= 2*minComponentSize)
					s.addAll(getSegmentations(str, cp.getIndex() + cp.getOffset(), end, nbMaxComponents -1, minComponentSize, addToDuplicate(cuttingPoints, cp, nbMaxComponents - 1)));
			}
		}
		return s;
	}
958031
def alignments(self):
        
        referenceIds = self.referenceIds
        dropUnmapped = self.dropUnmapped
        dropSecondary = self.dropSecondary
        dropSupplementary = self.dropSupplementary
        dropDuplicates = self.dropDuplicates
        keepQCFailures = self.keepQCFailures
        storeQueryIds = self.storeQueryIds
        filterRead = self.filterRead
        minScore = self.minScore
        maxScore = self.maxScore
        scoreTag = self.scoreTag

        if storeQueryIds:
            self.queryIds = queryIds = set()

        lastAlignment = None
        count = 0
        with samfile(self.filename) as samAlignment:
            for count, alignment in enumerate(samAlignment.fetch(), start=1):
                if storeQueryIds:
                    queryIds.add(alignment.query_name)

                if minScore is not None or maxScore is not None:
                    try:
                        score = alignment.get_tag(scoreTag)
                    except KeyError:
                        continue
                    else:
                        if ((minScore is not None and score < minScore) or
                                (maxScore is not None and score > maxScore)):
                            continue

                # Secondary and supplementary alignments may have a '*'
                # (pysam returns this as None) SEQ field, indicating that
                # the previous sequence should be used. This is best
                # practice according to section 2.5.2 of
                # https://samtools.github.io/hts-specs/SAMv1.pdf So we use
                # the last alignment query and quality strings if we get
                # None as a query sequence.
                if alignment.query_sequence is None:
                    if lastAlignment is None:
                        raise InvalidSAM(
                            'pysam produced an alignment (number %d) with no '
                            'query sequence without previously giving an '
                            'alignment with a sequence.' % count)
                    # Use the previous query sequence and quality. I'm not
                    # making the call to _hardClip dependent on
                    # alignment.cigartuples (as in the else clause below)
                    # because I don't think it's possible for
                    # alignment.cigartuples to be None in this case. If we
                    # have a second match on a query, then it must be
                    # aligned to something (i.e., it cannot be unmapped
                    # with no CIGAR string). The assertion will tell us if
                    # this is ever not the case.
                    assert alignment.cigartuples
                    (alignment.query_sequence,
                     alignment.query_qualities, _) = _hardClip(
                         lastAlignment.query_sequence,
                         lastAlignment.query_qualities,
                         alignment.cigartuples)
                else:
                    lastAlignment = alignment
                    if alignment.cigartuples:
                        (alignment.query_sequence,
                         alignment.query_qualities, _) = _hardClip(
                             alignment.query_sequence,
                             alignment.query_qualities,
                             alignment.cigartuples)

                if ((filterRead is None or
                     filterRead(Read(alignment.query_name,
                                     alignment.query_sequence,
                                     alignment.qual))) and
                    not (
                        (referenceIds and
                         alignment.reference_name not in referenceIds) or
                        (alignment.is_unmapped and dropUnmapped) or
                        (alignment.is_secondary and dropSecondary) or
                        (alignment.is_supplementary and dropSupplementary) or
                        (alignment.is_duplicate and dropDuplicates) or
                        (alignment.is_qcfail and not keepQCFailures))):
                    yield alignment

        self.alignmentCount = count
901253
def get_read_position(read, soft_clip_threshold):
    '''  '''
    is_spliced = False

    if read.is_reverse:
        pos = read.aend
        if read.cigar[-1][0] == 4:
            pos = pos + read.cigar[-1][1]
        start = read.pos

        if ('N' in read.cigarstring or
            (read.cigar[0][0] == 4 and
             read.cigar[0][1] > soft_clip_threshold)):

            cigar = read.cigar[::-1]
            is_spliced = find_splice(cigar)
    else:
        pos = read.pos
        if read.cigar[0][0] == 4:
            pos = pos - read.cigar[0][1]
        start = pos

        if ('N' in read.cigarstring or
            (read.cigar[-1][0] == 4 and
             read.cigar[-1][1] > soft_clip_threshold)):
            is_spliced = find_splice(read.cigar)

    return start, pos, is_spliced
525497
private int findSpanStart(TextCursor cursor, int limit) {
        for (int i = cursor.currentOffset; i < limit; i++) {
            char c = cursor.text.charAt(i);
            if (c == '*' || c == '_') {
                // Check prev and next symbols
                if (isGoodAnchor(cursor.text, i - 1) && isNotSymbol(cursor.text, i + 1, c)) {
                    return i;
                }
            }
        }
        return -1;
    }
1253522
def parse_segments(vcffile):
    
    from cStringIO import StringIO
    from cyvcf2 import VCF

    output = StringIO()
    for v in VCF(vcffile):
        chrom = v.CHROM
        start = v.start
        end = v.INFO.get('END') - 1
        cn, = v.format('CN')[0]
        print("\t".join(str(x) for x in (chrom, start, end, cn)), file=output)

    beds = BedTool(output.getvalue(), from_string=True)
    return beds
865761
def getReadAlignments(self, reference, start=None, end=None):
        
        return self._getReadAlignments(reference, start, end, self, None)
1042608
def findlinestarts(code, dup_lines=False):
    
    if not isinstance(code.co_lnotab, str):
        byte_increments = list(code.co_lnotab[0::2])
        line_increments = list(code.co_lnotab[1::2])
    else:
        byte_increments = [ord(c) for c in code.co_lnotab[0::2]]
        line_increments = [ord(c) for c in code.co_lnotab[1::2]]

    lastlineno = None
    lineno = code.co_firstlineno
    offset = 0
    for byte_incr, line_incr in zip(byte_increments, line_increments):
        if byte_incr:
            if (lineno != lastlineno or
                (dup_lines and 0 < byte_incr < 255)):
                yield (offset, lineno)
                lastlineno = lineno
                pass
            offset += byte_incr
            pass
        lineno += line_incr
    if (lineno != lastlineno or
        (dup_lines and 0 < byte_incr < 255)):
        yield (offset, lineno)
623020
def queriessam2dalignbed(cfg):
    
    datatmpd=cfg['datatmpd']
    alignmentbedp=cfg['alignmentbedp']
    dalignbedp=cfg['dalignbedp']
    logging.info(basename(dalignbedp))
    if not exists(alignmentbedp) or cfg['force']:
        #input/s
        queriessamps=glob(f'{datatmpd}/01_queries_queryl*.sam')
        for queriessamp in queriessamps:
            if stat(queriessamp).st_size != 0:
                samfile=pysam.AlignmentFile(queriessamp, "rb")
                dalignbed=pd.DataFrame(columns=bed_colns)
                for read in samfile.fetch():
                    algnids=[]
                    
#                     read_position=read.positions[0]
#                     tag_NM=read.get_tag('NM')                    
                    if len(read.positions)!=0:
                        read_position=read.positions[0]
                    else:
                        logging.error('no alignments found')
                        print(read)
                        continue
#                         read_position=[None]
                    try: 
                        tag_NM=read.get_tag('NM')
                    except: 
                        logging.error('no NM tag found')
                        print(read)
                        continue
#                         tag_NM=None
                    algnids.append(f"{read.reference_name}|{'-' if read.is_reverse else '+'}{read_position}|{read.cigarstring}|{tag_NM}")
                    if read.has_tag('XA'):
                        algnids+=['|'.join(s.split(',')) for s in read.get_tag('XA').split(';') if len(s.split(','))>1]
#                     print(len(algnids))
                    chroms=[]
                    starts=[]
                    ends=[]
                    algnids=algnids[:20]
                    NMs=[]
                    strands=[]
                    for a in algnids:
                        strand=a.split('|')[1][0]
                        chroms.append(a.split('|')[0])
                        if strand=='+':
                            offset=0
                        elif strand=='-':
                            offset=0                    
                        starts.append(int(a.split('|')[1][1:])+offset)
                        ends.append(int(a.split('|')[1][1:])+str2num(a.split('|')[2])+offset)
                        NMs.append(a.split('|')[3])
                        strands.append(strand)
                        del strand,offset
                    col2dalignbed={'chromosome':chroms,
                                   'start':starts,
                                   'end':ends,
                                   'id':algnids,
                                   'NM':NMs,
                                   'strand':strands}
                #     col2dalignbed=dict(zip(cols,[a.split('|')[0],a.split('|')[1],a.split('|')[2],a,a.split('|')[3],a.split('|')[4] for a in algnids]))
                    dalignbed_=pd.DataFrame(col2dalignbed)
                    dalignbed_['query id']=read.qname.replace('_',' ')
                    dalignbed = dalignbed.append(dalignbed_,ignore_index=True,sort=True)
                #     break
                samfile.close()
            else:
                logging.warning(f"file is empty: {queriessamp}")
        dalignbed.to_csv(dalignbedp,sep='\t')
        from rohan.dandage.io_nums import str2numorstr
        dalignbed['chromosome']=dalignbed.apply(lambda x : str2numorstr(x['chromosome']),axis=1)
        dalignbed=dalignbed.sort_values(['chromosome','start','end'], ascending=[True, True, True])
        dalignbed.loc[:,bed_colns].to_csv(alignmentbedp,sep='\t',
                        header=False,index=False,
                        chunksize=5000)
    return cfg
666038
def match(tgt, opts=None):
    '''
    
    '''
    if not opts:
        opts = __opts__
    nodegroups = opts.get('nodegroups', {})
    matchers = salt.loader.matchers(opts)

    if not isinstance(tgt, six.string_types) and not isinstance(tgt, (list, tuple)):
        log.error('Compound target received that is neither string, list nor tuple')
        return False
    log.debug('compound_match: %s ? %s', opts['id'], tgt)
    ref = {'G': 'grain',
           'P': 'grain_pcre',
           'I': 'pillar',
           'J': 'pillar_pcre',
           'L': 'list',
           'N': None,      # Nodegroups should already be expanded
           'S': 'ipcidr',
           'E': 'pcre'}
    if HAS_RANGE:
        ref['R'] = 'range'

    results = []
    opers = ['and', 'or', 'not', '(', ')']

    if isinstance(tgt, six.string_types):
        words = tgt.split()
    else:
        # we make a shallow copy in order to not affect the passed in arg
        words = tgt[:]

    while words:
        word = words.pop(0)
        target_info = salt.utils.minions.parse_target(word)

        # Easy check first
        if word in opers:
            if results:
                if results[-1] == '(' and word in ('and', 'or'):
                    log.error('Invalid beginning operator after "(": %s', word)
                    return False
                if word == 'not':
                    if not results[-1] in ('and', 'or', '('):
                        results.append('and')
                results.append(word)
            else:
                # seq start with binary oper, fail
                if word not in ['(', 'not']:
                    log.error('Invalid beginning operator: %s', word)
                    return False
                results.append(word)

        elif target_info and target_info['engine']:
            if 'N' == target_info['engine']:
                # if we encounter a node group, just evaluate it in-place
                decomposed = salt.utils.minions.nodegroup_comp(target_info['pattern'], nodegroups)
                if decomposed:
                    words = decomposed + words
                continue

            engine = ref.get(target_info['engine'])
            if not engine:
                # If an unknown engine is called at any time, fail out
                log.error(
                    'Unrecognized target engine "%s" for target '
                    'expression "%s"', target_info['engine'], word
                )
                return False

            engine_args = [target_info['pattern']]
            engine_kwargs = {}
            if target_info['delimiter']:
                engine_kwargs['delimiter'] = target_info['delimiter']

            results.append(
                six.text_type(matchers['{0}_match.match'.format(engine)](*engine_args, **engine_kwargs))
            )

        else:
            # The match is not explicitly defined, evaluate it as a glob
            results.append(six.text_type(matchers['glob_match.match'](word)))

    results = ' '.join(results)
    log.debug('compound_match %s ? "%s" => "%s"', opts['id'], tgt, results)
    try:
        return eval(results)  # pylint: disable=W0123
    except Exception:
        log.error(
            'Invalid compound target: %s for results: %s', tgt, results)
        return False
    return False