returns the nb of variants of sequences between x1 and x2
14146
def getNbVariants(self, x1, x2 = -1) :
		""""""
		if x2 == -1 :
			xx2 = len(self.defaultSequence)
		else :
			xx2 = x2
		
		nbP = 1
		for p in self.polymorphisms:
			if x1 <= p[0] and p[0] <= xx2 :
				nbP *= len(p[1])
		
		return nbP
942566
def filter_variant_sequences(
        variant_sequences,
        preferred_sequence_length,
        min_variant_sequence_coverage=MIN_VARIANT_SEQUENCE_COVERAGE,):
    
    variant_sequences = trim_variant_sequences(
        variant_sequences, min_variant_sequence_coverage)

    return filter_variant_sequences_by_length(
        variant_sequences=variant_sequences,
        preferred_sequence_length=preferred_sequence_length)
527085
def get_possible_splits(n):
    
    get_bin = lambda x, n: x >= 0 and str(bin(x))[2:].zfill(n) or "-" + str(bin(x))[3:].zfill(n)
    possible_splits = []
    for i in range(2**(n - 1)):
        possible_splits.append(get_bin(i, n - 1))
    return possible_splits
935885
def _bisect(seq):
        
        if len(seq) == 0:
            return [], None, []
        mid_idx = int(len(seq)/2)
        mid = seq[mid_idx]
        greater = seq[mid_idx+1:]
        lesser = seq[:mid_idx]
        return mid, greater, lesser
652116
def idx_sequence(self):
        
        return [x[1] for x in sorted(zip(self._record, list(range(len(self._record)))))]
1271521
def sequence_matcher_similarity(state_a, state_b):
        
        addrs_a = tuple(state_a.history.bbl_addrs)
        addrs_b = tuple(state_b.history.bbl_addrs)
        return SequenceMatcher(a=addrs_a, b=addrs_b).ratio()
112048
def compare_seqs_leven(seqs):
    
    A, B, ignore_gaps = seqs
    a, b = remove_gaps(A[1], B[1]) # actual sequences
    if len(a) != len(b):
        print('# reads are not the same length', file=sys.stderr)
        exit()
    pident = lr(a, b) * 100
    return A[0], B[0], pident
1195194
def sorensen(seq1, seq2):
	
	set1, set2 = set(seq1), set(seq2)
	return 1 - (2 * len(set1 & set2) / float(len(set1) + len(set2)))
664423
def split(self, sequence):
    

    major_idx = sequence.idx
    idx2 = 0
    for start, end in zip(major_idx[:-1], major_idx[1:]):
      idx1 = self.idx.index(start, idx2)
      idx2 = self.idx.index(end, idx2)
      seq = Sequence(self.text[start:end])
      seq.idx = [x-start for x in self.idx[idx1:idx2]]
      yield seq
175910
def all_pairs(seq1, seq2=None):
    
    if seq2 is None: seq2 = seq1
    for item1 in seq1:
        for item2 in seq2:
            yield (item1, item2)
725817
def chop_sequence(sequence, limit_length):
        
        return [sequence[i:i + limit_length] for i in range(0, len(sequence), limit_length)]
1015833
def hyp_2F1(a, b, c, z):
    
    return mp.hyp2f1(a, b, c, z)
974996
def binned_bitsets_proximity( f, chrom_col=0, start_col=1, end_col=2, strand_col=5, upstream=0, downstream=0 ):
    """"""
    last_chrom = None
    last_bitset = None
    bitsets = dict()
    for line in f:
        if line.startswith("#"): continue
#        print "input=%s" % ( line ),
        fields = line.split()
        strand = "+"
        if len(fields) >= strand_col + 1:
            if fields[strand_col] == "-": strand = "-"
        chrom = fields[chrom_col]
        if chrom != last_chrom:
            if chrom not in bitsets:
                bitsets[chrom] = BinnedBitSet( MAX )
            last_chrom = chrom
            last_bitset = bitsets[chrom]
        start, end = int( fields[start_col] ), int( fields[end_col] )
        if strand == "+":
            if upstream: start = max( 0, start - upstream )
            if downstream: end = min( MAX, end + downstream )
        if strand == "-":
            if upstream: end = min( MAX, end + upstream )
            if downstream: start = max( 0, start - downstream )
#        print "set: start=%d\tend=%d" % ( start, end )
        if end-start > 0:
            last_bitset.set_range( start, end-start )
    return bitsets
437650
def inner_join_impl(other, sequence):
    
    seq_dict = {}
    for element in sequence:
        seq_dict[element[0]] = element[1]
    seq_kv = seq_dict
    other_kv = dict(other)
    keys = seq_kv.keys() if len(seq_kv) < len(other_kv) else other_kv.keys()
    result = {}
    for k in keys:
        if k in seq_kv and k in other_kv:
            result[k] = (seq_kv[k], other_kv[k])
    return six.viewitems(result)
1266863
def iter_nearest_neighbours(umis, substr_idx):
    '''
    
    '''
    for i, u in enumerate(umis, 1):
        neighbours = set()
        for idx, substr_map in substr_idx.items():
            u_sub = u[slice(*idx)]
            neighbours = neighbours.union(substr_map[u_sub])
        neighbours.difference_update(umis[:i])
        for nbr in neighbours:
            yield u, nbr
896392
def reduce_sequences(object_a, object_b):
    
    def is_seq(obj):
        """Returns true if the object passed is a sequence."""
        return hasattr(obj, "__getitem__") or hasattr(obj, "__iter__")
    if object_a is None or object_b is None:
        return None
    elif is_seq(object_a) and is_seq(object_b):
        reduced = []
        for element_a, element_b in zip(object_a, object_b):
            if element_a is not None and element_b is not None:
                reduced.append(element_a + element_b)
            else:
                reduced.append(None)
        return reduced
    else:
        return object_a + object_b
916906
def diff(args):
    
    from jcvi.utils.table import banner

    p = OptionParser(diff.__doc__)
    p.add_option("--ignore_case", default=False, action="store_true",
            help="ignore case when comparing sequences [default: %default]")
    p.add_option("--ignore_N", default=False, action="store_true",
            help="ignore N and X's when comparing sequences [default: %default]")
    p.add_option("--ignore_stop", default=False, action="store_true",
            help="ignore stop codon when comparing sequences [default: %default]")
    p.add_option("--rc", default=False, action="store_true",
            help="also consider reverse complement [default: %default]")
    p.add_option("--quiet", default=False, action="store_true",
            help="don't output comparison details [default: %default]")

    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    afasta, bfasta = args

    afastan = len(Fasta(afasta))
    bfastan = len(Fasta(bfasta))

    if afastan == bfastan:
        print(green("Two sets contain the same number of sequences ({0}, {1})".\
                format(afastan, bfastan)))
    else:
        print(red("Two sets contain different number of sequences ({0}, {1})".\
                format(afastan, bfastan)))

    ah = SeqIO.parse(afasta, "fasta")
    bh = SeqIO.parse(bfasta, "fasta")

    problem_ids = []
    for arec, brec in zip(ah, bh):

        if opts.ignore_stop:
            arec.seq = arec.seq.rstrip("*")
            brec.seq = brec.seq.rstrip("*")

        asize, bsize = len(arec), len(brec)

        if not opts.quiet:
            print(banner(str(arec), [str(brec)]))
            if asize == bsize:
                print(green("Two sequence size match (%d)" % asize))
            else:
                print(red("Two sequence size do not match (%d, %d)" % (asize, bsize)))

        # print out the first place the two sequences diff
        fd = print_first_difference(arec, brec, ignore_case=opts.ignore_case,
                ignore_N=opts.ignore_N, rc=opts.rc, report_match=not opts.quiet)
        if not fd:
            logging.error("Two sets of sequences differ at `{0}`".format(arec.id))
            problem_ids.append("\t".join(str(x) for x in (arec.id, asize, bsize,
                    abs(asize - bsize))))

    if problem_ids:
        print(red("A total of {0} records mismatch.".format(len(problem_ids))))
        fw = must_open("Problems.ids", "w")
        print("\n".join(problem_ids), file=fw)