longest common sequence ratio

    >>> lcsr('ab', 'abcd')
    0.5
175760
def lcsr(s1, s2):
    '''
    '''
    if s1 == s2:
        return 1.0
    return llcs(s1, s2) / max(1, len(s1), len(s2))
914983
def bin2str(b):
    
    ret = []
    for pos in range(0, len(b), 8):
        ret.append(chr(int(b[pos:pos + 8], 2)))
    return ''.join(ret)
512801
def concat_align(fastas):
    
    # read in sequences
    fa2len = {}
    seqs = {}
    IDs = []
    for fasta in fastas:
        seqs[fasta] = {}
        for seq in parse_fasta(fasta):
            ID = seq[0].split('>')[1].split()[0]
            IDs.append(ID)
            seqs[fasta][ID] = seq[1]
        fa2len[fasta] = len(seq[1])
    # concat sequences
    IDs = set(IDs)
    concat = {}
    for fasta in fastas:
        for ID in IDs:
            if ID not in concat:
                concat[ID] = []
            if ID not in seqs[fasta]:
                concat[ID].append('-'*fa2len[fasta])
            else:
                concat[ID].append(seqs[fasta][ID])
    return concat
661743
def fasta(self):
        """"""
        max_line_length = 79
        fasta_str = '>{0}:{1}|PDBID|CHAIN|SEQUENCE\n'.format(
            self.parent.id.upper(), self.id)
        seq = self.sequence
        split_seq = [seq[i: i + max_line_length]
                     for i in range(0, len(seq), max_line_length)]
        for seq_part in split_seq:
            fasta_str += '{0}\n'.format(seq_part)
        return fasta_str
650900
def long2str(l):
    """"""
    if type(l) not in (types.IntType, types.LongType):
        raise ValueError('the input must be an integer')

    if l < 0:
        raise ValueError('the input must be greater than 0')
    s = ''
    while l:
        s = s + chr(l & 255)
        l >>= 8

    return s
101155
public static int cusparseXcsr2bsrNnz(
        cusparseHandle handle, 
        int dirA, 
        int m, 
        int n, 
        cusparseMatDescr descrA, 
        Pointer csrSortedRowPtrA, 
        Pointer csrSortedColIndA, 
        int blockDim, 
        cusparseMatDescr descrC, 
        Pointer bsrSortedRowPtrC, 
        Pointer nnzTotalDevHostPtr)
    {
        return checkResult(cusparseXcsr2bsrNnzNative(handle, dirA, m, n, descrA, csrSortedRowPtrA, csrSortedColIndA, blockDim, descrC, bsrSortedRowPtrC, nnzTotalDevHostPtr));
    }
1055703
def utr5(self):
        
        if not self.is_coding or len(self.exons) < 2: return (None, None)
        if self.strand == "+":
            s, e = (self.txStart, self.cdsStart)
        else:
            s, e = (self.cdsEnd, self.txEnd)
        if s == e: return (None, None)
        return s, e
230181
def sequenceToWord(sequence):
    
    reberString = ''
    for i in xrange(len(sequence)):
        index = np.where(sequence[i]==1.)[0][0]
        reberString += chars[index]
    return reberString
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)
738643
private function longestCommonPrefixLength($s, $t)
    {
        $length = min(strlen($s), strlen($t));
        $i = 0;
        while ($i < $length) {
            if ($s [$i] !== $t [$i]) {
                break;
            }
            $i++;
        }
        return $i;
    }
401831
def kmer_lca(seqs_path: 'path to (optionally gzipped) fasta/fastq input',
             fastq: 'input is fastq; disable autodetection' = False,
             progress: 'show progress bar (sent to stderr)' = False):
    '''
    
    '''
    conf = tictax.config()
    records = tictax.parse_seqs(seqs_path, fastq)
    print('Classifying sequences…', file=sys.stderr)
    asyncio.get_event_loop().run_until_complete(tictax.oc_classify(records,
                                                                   conf['one_codex_api_key'],
                                                                   progress,
                                                                   True))
    print('✓📌 ✓📌 ✓📌 ✓📌 ✓📌 ✓📌 ✓📌 ✓📌 ✓📌 ✓📌', file=sys.stderr)
150800
def _get_SEQRES_sequences(self):
        ''''''

        pdb_id = self.get_pdb_id()
        SEQRES_lines = self.parsed_lines["SEQRES"]

        modified_residue_mapping_3 = self.modified_residue_mapping_3
        # I commented this out since we do not need it for my current test cases
        #for k, v in self.modified_residues.iteritems():
        #    assert(v['modified_residue'] not in modified_residues)
        #    modified_residues[v['modified_residue']] = v['original_residue_3']

        for x in range(0, len(SEQRES_lines)):
            assert(SEQRES_lines[x][7:10].strip().isdigit())

        if not SEQRES_lines:
            #colortext.warning("WARNING: No SEQRES records were found. Kyle is trying to handle this gracefully, but Shane may need to fix it")
            return

        seqres_chain_order = []
        SEQRES_lines = [line[11:].rstrip() for line in SEQRES_lines] # we cannot strip the left characters as some cases e.g. 2MBP are missing chain identifiers

        # Collect all residues for all chains, remembering the chain order
        chain_tokens = {}
        for line in SEQRES_lines:
            chainID = line[0]
            if missing_chain_ids.get(self.pdb_id):
                chainID = missing_chain_ids[self.pdb_id]
            if chainID not in seqres_chain_order:
                seqres_chain_order.append(chainID)
            chain_tokens[chainID] = chain_tokens.get(chainID, [])
            chain_tokens[chainID].extend(line[6:].strip().split())

        sequences = {}
        self.chain_types = {}

        canonical_acid_types = set(residue_type_3to1_map.keys())
        canonical_acid_types.remove('UNK')
        assert(len(canonical_acid_types) == 20)

        for chain_id, tokens in chain_tokens.iteritems():

            # Determine whether the chain is DNA, RNA, or a protein chain
            # 1H38 is a good test for this - it contains DNA (chains E and G and repeated by H, K, N, J, M, P), RNA (chain F, repeated by I, L, O) and protein (chain D, repeated by A,B,C) sequences
            # 1ZC8 is similar but also has examples of DU
            # 4IHY has examples of DI (I is inosine)
            # 2GRB has RNA examples of I and U
            # 1LRP has protein chains with only CA atoms
            # This will throw an exception when a non-canonical is found which is not listed in basics.py. In that case, the list in basics.py should be updated.

            chain_type = None
            set_of_tokens = set(tokens)
            if (set_of_tokens.union(all_recognized_dna) == all_recognized_dna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'DNA'
            elif (set_of_tokens.union(all_recognized_rna) == all_recognized_rna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'RNA'
            elif len(set_of_tokens) == 1 and 'UNK' in set_of_tokens:
                chain_type = 'Unknown'
            elif not(set_of_tokens.intersection(canonical_acid_types)):
                # Zero canonical residues may imply a ligand or a heterogen chain
                chain_type = PDB._determine_heterogen_chain_type(set_of_tokens)
            else:
                assert(len(set(tokens).intersection(dna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                chain_type = 'Protein'
                if not self.chain_atoms.get(chain_id):
                    # possible for biological unit files
                    continue
                if self.chain_atoms[chain_id] == set(['CA']):
                    chain_type = 'Protein skeleton'

            # Get the sequence, mapping non-canonicals to the appropriate letter
            self.chain_types[chain_id] = chain_type

            sequence = []
            if chain_type == 'DNA':
                for r in tokens:
                    if dna_nucleotides_2to1_map.get(r):
                        sequence.append(dna_nucleotides_2to1_map[r])
                    else:
                        if non_canonical_dna.get(r):
                            sequence.append(non_canonical_dna[r])
                        else:
                            raise Exception("Unknown DNA residue %s." % r)
            elif chain_type == 'RNA':
                for r in tokens:
                    if r in rna_nucleotides:
                        sequence.append(r)
                    else:
                        if non_canonical_rna.get(r):
                            sequence.append(non_canonical_rna[r])
                        else:
                            raise Exception("Unknown RNA residue %s." % r)
            else:
                token_counter = 0
                for r in tokens:
                    token_counter += 1
                    if residue_type_3to1_map.get(r):
                        sequence.append(residue_type_3to1_map[r])
                    else:

                        if self.modified_residue_mapping_3.get(r):
                            sequence.append(residue_type_3to1_map[self.modified_residue_mapping_3.get(r)])
                        elif non_canonical_amino_acids.get(r):
                            #print('Mapping non-canonical residue %s to %s.' % (r, non_canonical_amino_acids[r]))
                            #print(SEQRES_lines)
                            #print(line)
                            sequence.append(non_canonical_amino_acids[r])
                        elif r == 'UNK':
                            continue
                        # Skip these residues
                        elif r == 'ACE' and token_counter == 1:
                            # Always allow ACE as the first residue of a chain
                            sequence.append('X')
                        elif r == 'ACE' and pdb_id in cases_with_ACE_residues_we_can_ignore:
                            sequence.append('X')
                            #continue
                        # End of skipped residues
                        else:
                            #print(modified_residue_mapping_3)
                            if modified_residue_mapping_3.get(r):
                                if modified_residue_mapping_3[r] == 'UNK':
                                    sequence.append('X')
                                else:
                                    assert(modified_residue_mapping_3[r] in residue_types_3)
                                    sequence.append(residue_type_3to1_map[modified_residue_mapping_3[r]])
                            else:
                                raise Exception("Unknown protein residue %s in chain %s." % (r, chain_id))
            sequences[chain_id] = "".join(sequence)

        self.seqres_chain_order = seqres_chain_order

        # Create Sequence objects for the SEQRES sequences
        for chain_id, sequence in sequences.iteritems():
            self.seqres_sequences[chain_id] = Sequence.from_sequence(chain_id, sequence, self.chain_types[chain_id])
830697
public static List<int[]> hunks(String s0, String s1) {
        List<int[]> lcs = longestCommonSubsequence(s0 ,s1);
        List<int[]> hunks = new ArrayList<int[]>();
        int inf0 = -1;
        int inf1 = -1;
        int last0 = -1;
        int last1 = -1;
        for (int i = 0; i < lcs.size(); i++) {
            int[] match = lcs.get(i);
            if (inf0 == -1 || inf1 == -1) {
                inf0 = match[0];
                inf1 = match[1];
            } else if (last0 + 1 != match[0] || last1 + 1 != match[1]) {
                hunks.add(new int[] {inf0, last0 + 1, inf1, last1 + 1});
                inf0 = match[0];
                inf1 = match[1];
            } else if (i == lcs.size() - 1) {
                hunks.add(new int[] {inf0, match[0] + 1, inf1, match[1] + 1});
                break;
            }
            last0 = match[0];
            last1 = match[1];
        }
        return hunks;
    }
1365326
def sift4(s1, s2, max_offset=5):
    
    t1, t2 = list(s1), list(s2)
    l1, l2 = len(t1), len(t2)

    if not s1:
        return l2

    if not s2:
        return l1

    # Cursors for each string
    c1, c2 = 0, 0

    # Largest common subsequence
    lcss = 0

    # Local common substring
    local_cs = 0

    # Number of transpositions ('ab' vs 'ba')
    trans = 0

    # Offset pair array, for computing the transpositions
    offsets = []

    while c1 < l1 and c2 < l2:
        if t1[c1] == t2[c2]:
            local_cs += 1

            # Check if current match is a transposition
            is_trans = False
            i = 0
            while i < len(offsets):
                ofs = offsets[i]
                if c1 <= ofs['c1'] or c2 <= ofs['c2']:
                    is_trans = abs(c2-c1) >= abs(ofs['c2'] - ofs['c1'])
                    if is_trans:
                        trans += 1
                    elif not ofs['trans']:
                        ofs['trans'] = True
                        trans += 1
                    break
                elif c1 > ofs['c2'] and c2 > ofs['c1']:
                    del offsets[i]
                else:
                    i += 1
            offsets.append({
                'c1': c1,
                'c2': c2,
                'trans': is_trans
            })

        else:
            lcss += local_cs
            local_cs = 0
            if c1 != c2:
                c1 = c2 = min(c1, c2)

            for i in range(max_offset):
                if c1 + i >= l1 and c2 + i >= l2:
                    break
                elif c1 + i < l1 and s1[c1+i] == s2[c2]:
                    c1 += i - 1
                    c2 -= 1
                    break

                elif c2 + i < l2 and s1[c1] == s2[c2 + i]:
                    c2 += i - 1
                    c1 -= 1
                    break

        c1 += 1
        c2 += 1

        if c1 >= l1 or c2 >= l2:
            lcss += local_cs
            local_cs = 0
            c1 = c2 = min(c1, c2)

    lcss += local_cs
    return round(max(l1, l2) - lcss + trans)
41952
def guess_phred_format(self):
        """"""
        # Possibilities #
        self.intervals = {
            'Sanger':       (33,  74),         # <- This one is the gold standard
            'Solexa':       (59, 105),
            'Illumina-1.3': (64, 105),         # <- These were abandoned after they wised up.
            'Illumina-1.5': (67, 105),         # <- These were abandoned after they wised up.
        }
        # Initialize variables #
        glob_min  = 9999
        glob_max  = 0
        take_next = False
        # Stop after #
        max_sequences_to_consider = 2000000
        # Error message #
        message = "Multiple PHRED encodings possible for '%s'.\n"
        # Loop #
        self.open()
        for i, line in enumerate(self.handle):
            # Are we on the PHRED line? #
            if not take_next:
                take_next = line.startswith('+')
                continue
            # Now we are on the PHRED line! #
            take_next = False
            current_min, current_max = self.get_qual_range(line.rstrip('\n'))
            if current_min < glob_min or current_max > glob_max:
                glob_min = min(current_min, glob_min)
                glob_max = max(current_max, glob_max)
                valid_encodings = [e for e, r in self.intervals.items() if glob_min >= r[0] and glob_max <= r[1]]
                if len(valid_encodings) == 0:
                    message = "Illegal PHRED encoding for '%s'.\n"
                    break
                if len(valid_encodings) == 1:
                    return valid_encodings[0]
            # Stop condition #
            if i >= max_sequences_to_consider: break
        # It didn't work #
        message += "Maximum detected value is: %s.\n"
        message += "Minimum detected value is: %s.\n"
        message += "Possibilities are: %s."
        message  = message % (self, glob_max, glob_min, ' or '.join(valid_encodings))
        raise Exception(message)
1198343
def normalize_cjk_fullwidth_ascii(seq: str) -> str:
    

    def convert(char: str) -> str:
        code_point = ord(char)
        if not 0xFF01 <= code_point <= 0xFF5E:
            return char
        return chr(code_point - 0xFEE0)

    return ''.join(map(convert, seq))
123334
function parse_RgceLoc(blob, length, opts) {
	if(opts && opts.biff >= 2 && opts.biff <= 5) return parse_RgceLoc_BIFF2(blob, length, opts);
	var r = blob.read_shift(opts && opts.biff == 12 ? 4 : 2);
	var c = parse_ColRelU(blob, 2);
	return {r:r, c:c[0], cRel:c[1], rRel:c[2]};
}