m.bestscanseq(seq) -- Return score,sequence of the best match to the motif in the supplied sequence
194720
def bestscanseq(self,seq):
        
        matches,endpoints,scores = self.scan(seq,-100)
        t = zip(scores,matches)
        t.sort()
        bestseq   = t[-1][1]
        bestscore = t[-1][0]
        return bestscore, bestseq
318212
def findBest(self, pattern):
        
        findFailedRetry = True
        while findFailedRetry:
            best_match = None
            all_matches = self.findAll(pattern)
            for match in all_matches:
                if best_match is None or best_match.getScore() < match.getScore():
                    best_match = match
            self._lastMatch = best_match
            if best_match is not None:
                break
            path = pattern.path if isinstance(pattern, Pattern) else pattern
            findFailedRetry = self._raiseFindFailed("Could not find pattern '{}'".format(path))
            if findFailedRetry:
                time.sleep(self._repeatWaitTime)
        return best_match
706195
def seq_info(names, id2names, insertions, sequences):
    
    seqs = {} # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns]], ...]]
    for name in names:
        id = id2names[name]
        gene = name.split('fromHMM::', 1)[0].rsplit(' ', 1)[1]
        model = name.split('fromHMM::', 1)[1].split('=', 1)[1].split()[0]
        i_gene_pos = insertions[id] # coordinates of each insertion wrt gene
        i_model_pos = name.split('fromHMM::', 1)[1].split('model-pos(ins-len)=')[1].split()[0].split(';') # model overlap
        i_info = []
        for i, ins in enumerate(i_gene_pos):
            model_pos = i_model_pos[i].split('-')[1].split('(')[0]
            length = i_model_pos[i].split('(')[1].split(')')[0]
            iheader = '>%s_%s insertion::seq=%s type=insertion strand=n/a gene-pos=%s-%s model-pos=%s'\
                    % (id, (i + 1), (i + 1), ins[0], ins[1], model_pos)
            iseq = sequences[id][1][ins[0]:(ins[1] + 1)]
            iseq = [iheader, iseq]
            info = [ins, model_pos, length, iseq, [], []]
            i_info.append(info)
        seqs[id] = [gene, model, i_info]
    return seqs
569359
def best_model(seq2hmm):
    
    for seq in seq2hmm:
        best = []
        for model in seq2hmm[seq]:
            best.append([model, sorted([i[-1] for i in seq2hmm[seq][model]], reverse = True)[0]])
        best_model = sorted(best, key = itemgetter(1), reverse = True)[0][0]
        seq2hmm[seq] = [best_model] + [seq2hmm[seq][best_model]]
    return seq2hmm
1319415
def short_sequence_repeat_extractor(string, min_length=1):
    
    length = len(string)

    k_max = length // 2 + 1
    if k_max > THRESHOLD:
        k_max = THRESHOLD // 2

    repeats = []

    i = 0
    last_repeat = i
    while i < length:
        max_count = 0
        max_k = 1
        for k in range(min_length, k_max):
            count = 0
            for j in range(i + k, length - k + 1, k):
                if string[i:i + k] != string[j:j + k]:
                    break
                count += 1

            if count > 0 and count >= max_count:
                max_count = count
                max_k = k

        if max_count > 0:
            if last_repeat < i:
                repeats.append(Repeat(last_repeat, i))
            repeats.append(Repeat(i, i + max_k, max_count))
            last_repeat = i + max_k * (max_count + 1)

        i += max_k * (max_count + 1)

    if last_repeat < i:
        repeats.append(Repeat(last_repeat, i))

    return repeats
732795
def pick(seq, func, maxobj=None):
    """"""
    maxscore = None
    for obj in seq:
        score = func(obj)
        if maxscore is None or maxscore < score:
            (maxscore, maxobj) = (score, obj)
    return maxobj
599145
def get_sequences_from_cluster(c1, c2, data):
    """"""
    seqs1 = data[c1]['seqs']
    seqs2 = data[c2]['seqs']
    seqs = list(set(seqs1 + seqs2))
    names = []
    for s in seqs:
        if s in seqs1 and s in seqs2:
            names.append("both")
        elif s in seqs1:
            names.append(c1)
        else:
            names.append(c2)
    return seqs, names
949192
def find(find_seq, elem_list):
        
        seq_pos = 0
        for index, elem in enumerate(elem_list):
            if Sequence.match(elem, find_seq[seq_pos]):
                seq_pos += 1
                if seq_pos == len(find_seq):  # found matching sequence
                    return index - seq_pos + 1
            else:  # exited sequence
                seq_pos = 0
        raise LookupError('Failed to find sequence in elem_list')
1000013
def seqannotation(self, seqrecord, allele, loc):
        
        #seqrecord = self.seqrecord(allele, loc)
        complete_annotation = get_features(seqrecord)
        annotation = Annotation(annotation=complete_annotation,
                                method='match',
                                complete_annotation=True)

        if self.alignments:
            alignment = {f: self.annoated_alignments[loc][allele][f]['Seq']
                         for f in self.annoated_alignments[loc][allele].keys()}
            annotation.aligned = alignment

        return annotation
604848
def locate(self, pattern):
        '''

        '''
        if self.circular:
            if len(pattern) >= 2 * len(self):
                raise ValueError('Search pattern longer than searchable ' +
                                 'sequence.')
            seq = self + self[:len(pattern) - 1]
            return super(NucleicAcid, seq).locate(pattern)
        else:
            return super(NucleicAcid, self).locate(pattern)
502434
def _py_expand_long(subsequence, sequence, max_l_dist):
    """"""
    # The additional optimization in this version is to limit the part of
    # the sub-sequence inspected for each sequence character.  The start and
    # end of the iteration are limited to the range where the scores are
    # smaller than the maximum allowed distance.  Additionally, once a good
    # expansion has been found, the range is further reduced to where the
    # scores are smaller than the score of the best expansion found so far.

    subseq_len = len(subsequence)
    if subseq_len == 0:
        return (0, 0)

    # Initialize the scores array with values for just skipping sub-sequence
    # chars.
    scores = list(range(1, subseq_len + 1))

    min_score = subseq_len
    min_score_idx = -1
    max_good_score = max_l_dist
    new_needle_idx_range_start = 0
    new_needle_idx_range_end = subseq_len - 1

    for seq_index, char in enumerate(sequence):
        # calculate scores, one for each character in the sub-sequence
        needle_idx_range_start = new_needle_idx_range_start
        needle_idx_range_end = min(subseq_len, new_needle_idx_range_end + 1)

        a = seq_index
        c = a + 1

        if c <= max_good_score:
            new_needle_idx_range_start = 0
            new_needle_idx_range_end = 0
        else:
            new_needle_idx_range_start = None
            new_needle_idx_range_end = -1

        for subseq_index in range(needle_idx_range_start, needle_idx_range_end):
            b = scores[subseq_index]
            c = scores[subseq_index] = min(
                a + (char != subsequence[subseq_index]),
                b + 1,
                c + 1,
            )
            a = b

            if c <= max_good_score:
                if new_needle_idx_range_start is None:
                    new_needle_idx_range_start = subseq_index
                new_needle_idx_range_end = max(
                    new_needle_idx_range_end,
                    subseq_index + 1 + (max_good_score - c),
                )

        # bail early when it is impossible to find a better expansion
        if new_needle_idx_range_start is None:
            break

        # keep the minimum score found for matches of the entire sub-sequence
        if needle_idx_range_end == subseq_len and c <= min_score:
            min_score = c
            min_score_idx = seq_index
            if min_score < max_good_score:
                max_good_score = min_score

    return (min_score, min_score_idx + 1) if min_score <= max_l_dist else (None, None)
939057
def mask_seq(self,seq):
        
        masked = ''
        matches, endpoints, scores = self.scan(seq)
        cursor = 0
        for start, stop in endpoints:
            masked = masked + seq[cursor:start] + 'N'*self.width
            cursor = stop+1
        masked = masked + seq[cursor:]
        return masked
49950
def bestscore(self,seq, fwd=''):
        
        matches, endpoints, scores = self._scan(seq,threshold=-100000,forw_only=fwd)
        if scores: return max(scores)
        else:      return -1000
618869
public static <S extends Sequence<S>> int calculateScore(S seq1, Range seq1Range, Mutations<S> mutations,
                                                             AlignmentScoring<S> scoring) {
        if (scoring instanceof LinearGapAlignmentScoring)
            return calculateScore(seq1, seq1Range, mutations, (LinearGapAlignmentScoring<S>) scoring);
        else if (scoring instanceof AffineGapAlignmentScoring)
            return calculateScore(seq1, seq1Range, mutations, (AffineGapAlignmentScoring<S>) scoring);
        else
            throw new IllegalArgumentException("Unknown scoring type");
    }
144158
def longest_increasing_subsequence(sequence):
    
    length = len(sequence)
    counts = [1 for _ in range(length)]
    for i in range(1, length):
        for j in range(0, i):
            if sequence[i] > sequence[j]:
                counts[i] = max(counts[i], counts[j] + 1)
                print(counts)
    return max(counts)
638711
def read_next_maf( file, species_to_lengths=None, parse_e_rows=False ):
    
    alignment = Alignment(species_to_lengths=species_to_lengths)
    # Attributes line
    line = readline( file, skip_blank=True )
    if not line: return None
    fields = line.split() 
    if fields[0] != 'a': raise Exception("Expected 'a ...' line")
    alignment.attributes = parse_attributes( fields[1:] )
    if 'score' in alignment.attributes:
        alignment.score = alignment.attributes['score']
        del alignment.attributes['score']
    else:
        alignment.score = 0
    # Sequence lines
    last_component = None
    while 1:
        line = readline( file )
        # EOF or Blank line terminates alignment components
        if not line or line.isspace(): break
        if line.isspace(): break 
        # Parse row
        fields = line.split()
        if fields[0] == 's':
            # An 's' row contains sequence for a component
            component = Component()
            component.src = fields[1]
            component.start = int( fields[2] )
            component.size = int( fields[3] )
            component.strand = fields[4]
            component.src_size = int( fields[5] )
            if len(fields) > 6: component.text = fields[6].strip()
            # Add to set
            alignment.add_component( component )
            last_component = component
        elif fields[0] == 'e':
            # An 'e' row, when no bases align for a given species this tells
            # us something about the synteny 
            if parse_e_rows:
                component = Component()
                component.empty = True
                component.src = fields[1]
                component.start = int( fields[2] )
                component.size = int( fields[3] )
                component.strand = fields[4]
                component.src_size = int( fields[5] )
                component.text = None
                synteny = fields[6].strip()
                assert len( synteny ) == 1, \
                    "Synteny status in 'e' rows should be denoted with a single character code"
                component.synteny_empty = synteny
                alignment.add_component( component )
                last_component = component
        elif fields[0] == 'i':
            # An 'i' row, indicates left and right synteny status for the 
            # previous component, we hope ;)
            assert fields[1] == last_component.src, "'i' row does not follow matching 's' row"
            last_component.synteny_left = ( fields[2], int( fields[3] ) )
            last_component.synteny_right = ( fields[4], int( fields[5] ) )
        elif fields[0] == 'q':
            assert fields[1] == last_component.src, "'q' row does not follow matching 's' row"
            # TODO: Should convert this to an integer array?
            last_component.quality = fields[2]
            
    return alignment
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)