Takes a dna sequence with polymorphismes and returns all the possible sequences that it can yield
162581
def getSequenceCombinaisons(polymorphipolymorphicDnaSeqSeq, pos = 0) :
	""""""

	if type(polymorphipolymorphicDnaSeqSeq) is not types.ListType :
		seq = list(polymorphipolymorphicDnaSeqSeq)
	else :
		seq = polymorphipolymorphicDnaSeqSeq

	if pos >= len(seq) :
		return [''.join(seq)]

	variants = []
	if seq[pos] in polymorphicNucleotides :
		chars = decodePolymorphicNucleotide(seq[pos])
	else :
		chars = seq[pos]#.split('/')

	for c in chars :
		rseq = copy.copy(seq)
		rseq[pos] = c
		variants.extend(getSequenceCombinaisons(rseq, pos + 1))

	return variants
677102
def scan(self, seqs, nreport=100, scan_rc=True, normalize=False):
        

        if not self.threshold:
            sys.stderr.write(
                "Using default threshold of 0.95. "
                "This is likely not optimal!\n"
                )
            self.set_threshold(threshold=0.95)

        seqs = as_fasta(seqs, genome=self.genome)
           
        it = self._scan_sequences(seqs.seqs, 
                    nreport, scan_rc)
       
        if normalize:
            if len(self.meanstd) == 0:
                self.set_meanstd()
            mean_std = [self.meanstd.get(m_id) for m_id in self.motif_ids]
            means = [x[0] for x in  mean_std]
            stds = [x[1] for x in  mean_std]


        for result in it:
            if normalize:
                zresult = [] 
                for i,mrow in enumerate(result):
                    mrow = [((x[0] - means[i]) / stds[i], x[1], x[2]) for x in mrow]
                    zresult.append(mrow)
                yield zresult
            else:
                yield result
230638
def gt_bases(self):
        
        result = []
        for a in self.gt_alleles:
            if a is None:
                result.append(None)
            elif a == 0:
                result.append(self.site.REF)
            else:
                result.append(self.site.ALT[a - 1].value)
        return tuple(result)
394493
def sequence_type(seq):
    '''

    '''
    if isinstance(seq, coral.DNA):
        material = 'dna'
    elif isinstance(seq, coral.RNA):
        material = 'rna'
    elif isinstance(seq, coral.Peptide):
        material = 'peptide'
    else:
        raise ValueError('Input was not a recognized coral.sequence object.')
    return material
650083
def get_single_allele_from_reads(allele_reads):
    
    allele_reads = list(allele_reads)

    if len(allele_reads) == 0:
        raise ValueError("Expected non-empty list of AlleleRead objects")

    seq = allele_reads[0].allele
    if any(read.allele != seq for read in allele_reads):
        raise ValueError("Expected all AlleleRead objects to have same allele '%s', got %s" % (
            seq, allele_reads))
    return seq
277183
def rs_generator_poly_all(max_nsym, fcr=0, generator=2):
    ''''''
    g_all = {}
    g_all[0] = g_all[1] = [1]
    for nsym in xrange(max_nsym):
        g_all[nsym] = rs_generator_poly(nsym, fcr, generator)
    return g_all
52409
def variants(self, case_id, skip=0, count=1000, filters=None):
        """"""
        filters = filters or {}
        logger.debug("Fetching case with case_id: {0}".format(case_id))
        case_obj = self.case(case_id)
        plugin, case_id = self.select_plugin(case_obj)
        self.filters = plugin.filters

        gene_lists = (self.gene_list(list_id) for list_id
                      in filters.get('gene_lists', []))
        nested_geneids = (gene_list.gene_ids for gene_list in gene_lists)
        gene_ids = set(itertools.chain.from_iterable(nested_geneids))

        if filters.get('gene_ids'):
            filters['gene_ids'].extend(gene_ids)
        else:
            filters['gene_ids'] = gene_ids
        variants = plugin.variants(case_id, skip, count, filters)
        return variants
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())
1176321
def generate_antisense_sequence(sequence):
    """"""
    dna_antisense = {
        'A': 'T',
        'T': 'A',
        'C': 'G',
        'G': 'C'
    }
    antisense = [dna_antisense[x] for x in sequence[::-1]]
    return ''.join(antisense)
587628
def encodePolymorphicNucleotide(polySeq) :
	"""'"""
	
	if type(polySeq) is types.StringType :
		if polySeq.find("/") < 0 :
			sseq = list(polySeq)
		else :
			sseq = polySeq.split('/')
			
	else :
		sseq = polySeq
	
	seq = []
	for n in sseq :
		try :
			for n2 in polymorphicNucleotides[n] :
				seq.append(n2)
		except KeyError :
			seq.append(n)
	
	seq = set(seq)
	
	if len(seq) == 4:
		return 'N'
	elif len(seq) == 3 :
		if 'T' not in seq :
			return 'V'
		elif 'G' not in seq :
			return 'H'
		elif 'C' not in seq :
			return 'D'
		elif 'A' not in seq :
			return 'B'
	elif len(seq) == 2 :
		if 'A' in seq and 'G' in seq :
			return 'R'
		elif 'C' in seq and 'T' in seq :
			return 'Y'
		elif 'A' in seq and 'C' in seq :
			return 'M'
		elif 'T' in seq and 'G' in seq :
			return 'K'
		elif 'A' in seq and 'T' in seq :
			return 'W'
		elif 'C' in seq and 'G' in seq :
			return 'S'
	elif polySeq[0] in nucleotides :
		return polySeq[0]
	else :
		raise UnknownNucleotide(polySeq)
165069
def from_sequence(cls, seq, order=None, **kwargs):
        """"""
        for arg in seq:
            yield cls(arg, order=order, **kwargs)
1109120
def protein_sequences_generator_to_dataframe(variant_and_protein_sequences_generator):
    
    return dataframe_from_generator(
        element_class=ProteinSequence,
        variant_and_elements_generator=variant_and_protein_sequences_generator,
        converters=dict(
            gene=lambda x: ";".join(x)))
805578
def from23andme(args):
    
    p = OptionParser(from23andme.__doc__)
    p.set_ref()
    opts, args = p.parse_args(args)

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

    txtfile, seqid = args
    ref_dir = opts.ref
    fastafile = op.join(ref_dir, "hs37d5.fa")
    fasta = Fasta(fastafile)

    pf = txtfile.rsplit(".", 1)[0]
    px = CM[seqid]
    chrvcf = pf + ".{0}.vcf".format(px)
    legend = op.join(ref_dir, "1000GP_Phase3/{0}.rsids".format(px))
    register = read_rsid(seqid, legend)

    fw = open(chrvcf, "w")
    print(get_vcfstanza(fastafile, fasta, txtfile), file=fw)

    fp = open(txtfile)
    seen = set()
    duplicates = skipped = missing = 0
    for row in fp:
        if row[0] == '#':
            continue
        rsid, chr, pos, genotype = row.split()
        if chr != seqid:
            continue
        pos = int(pos)
        if (chr, pos) in seen:
            duplicates += 1
            continue
        seen.add((chr, pos))
        genotype = list(genotype)
        if "-" in genotype:  # missing daa
            missing += 1
            continue

        # Y or MT
        if not register:
            assert len(genotype) == 1
            ref = fasta[chr][pos - 1].seq.upper()
            if "D" in genotype or "I" in genotype:
                skipped += 1
                continue
            genotype = genotype[0]
            code = "0/0" if ref == genotype else "1/1"
            alt = "." if ref == genotype else genotype
            print("\t".join(str(x) for x in \
                    (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)
            continue

        # If rsid is seen in the db, use that
        if rsid in register:
            pos, ref, alt = register[rsid]
        elif pos in register:
            pos, ref, alt = register[pos]
        else:
            skipped += 1  # Not in reference panel
            continue

        assert fasta[chr][pos - 1:pos + len(ref) - 1].seq.upper() == ref
        # Keep it bi-allelic
        not_seen = [x for x in alt if x not in genotype]
        while len(alt) > 1 and not_seen:
            alt.remove(not_seen.pop())
        if len(alt) > 1:
            alt = [alt[0]]
        alleles = [ref] + alt

        if len(genotype) == 1:
            genotype = [genotype[0]] * 2

        alt = ",".join(alt) or "."
        if "D" in genotype or "I" in genotype:
            max_allele = max((len(x), x) for x in alleles)[1]
            alleles = [("I" if x == max_allele else "D") for x in alleles]
            assert "I" in alleles and "D" in alleles
        a, b = genotype
        try:
            ia, ib = alleles.index(a), alleles.index(b)
        except ValueError:  # alleles not seen
            logging.error("{0}: alleles={1}, genotype={2}".\
                            format(rsid, alleles, genotype))
            skipped += 1
            continue
        code = "/".join(str(x) for x in sorted((ia, ib)))

        print("\t".join(str(x) for x in \
                (chr, pos, rsid, ref, alt, ".", ".", "PR", "GT", code)), file=fw)

    logging.debug("duplicates={0} skipped={1} missing={2}".\
                    format(duplicates, skipped, missing))
811157
def parse_codons(ref, start, end, strand):
    
    codon = []
    c = cycle([1, 2, 3])
    ref = ref[start - 1:end]
    if strand == -1:
        ref = rc_stats(ref)
    for pos in ref:
        n = next(c)
        codon.append(pos)
        if n == 3:
            yield codon
            codon = []
1270789
def group_unique_sequences(
        allele_reads,
        max_prefix_size=None,
        max_suffix_size=None):
    
    groups = defaultdict(set)
    for r in allele_reads:
        prefix = r.prefix
        allele = r.allele
        suffix = r.suffix
        if max_prefix_size and len(prefix) > max_prefix_size:
            prefix = prefix[-max_prefix_size:]
        if max_suffix_size and len(suffix) > max_suffix_size:
            suffix = suffix[:max_suffix_size]
        key = (prefix, allele, suffix)
        groups[key].add(r)
    return groups
487847
def from_sequence(cls, sequence, phos_3_prime=False):
        
        strand1 = NucleicAcidStrand(sequence, phos_3_prime=phos_3_prime)
        duplex = cls(strand1)
        return duplex
352769
def expand_alleles(p, tolerance=0):
    
    _p = set()
    for x in p:
        _p |= set(range(x - tolerance, x + tolerance + 1))
    return _p