Creates a DNA duplex from a nucleotide sequence.

        Parameters
        ----------
        sequence: str
            Nucleotide sequence.
        phos_3_prime: bool, optional
            If false the 5' and the 3' phosphor will be omitted.
487847
def from_sequence(cls, sequence, phos_3_prime=False):
        
        strand1 = NucleicAcidStrand(sequence, phos_3_prime=phos_3_prime)
        duplex = cls(strand1)
        return duplex
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 = []
217486
def from_sequence(cls, seq):
        
        args = tuple(islice(seq, 4))
        if len(args) > 3:
            raise ValueError('When creating a Point from sequence, it '
                             'must not have more than 3 items.')
        return cls(*args)
1208257
def from_variant_and_transcript_and_sequence_key(
            cls, variant, transcript, sequence_key):
        
        # get the interbase range of offsets which capture all reference
        # bases modified by the variant
        variant_start_offset, variant_end_offset = \
            interbase_range_affected_by_variant_on_transcript(
                variant=variant,
                transcript=transcript)

        start_codon_idx = min(transcript.start_codon_spliced_offsets)

        # skip any variants which occur in the 5' UTR or overlap the start codon
        # since
        #   (1) UTR variants have unclear coding effects and
        #   (2) start-loss variants may result in a new start codon / reading
        #       frame but we're not sure which one!
        if variant_start_offset < start_codon_idx + 3:
            logger.info(
                "Skipping transcript %s for variant %s, must be after start codon",
                transcript.name,
                variant)
            return None

        stop_codon_idx = min(transcript.stop_codon_spliced_offsets)

        # skip variants which affect the 3' UTR of the transcript since
        # they don't have obvious coding effects on the protein sequence
        if variant_start_offset >= stop_codon_idx + 3:
            logger.info(
                "Skipping transcript %s for variant %s, occurs in 3' UTR",
                transcript,
                variant)
            return None

        n_prefix = len(sequence_key.sequence_before_variant_locus)
        prefix_start_idx = variant_start_offset - n_prefix
        n_bases_between_start_and_variant = variant_start_offset - start_codon_idx
        n_full_codons_before_variant = n_bases_between_start_and_variant // 3

        # if the sequence before the variant contains more bases than the
        # distance to the start codon, then by definition it must contain
        # some untranslated bases
        contains_five_prime_utr = (n_prefix > n_bases_between_start_and_variant)
        # allows for the possibility that the first base in the sequence might
        # be the first nucleotide of the start codon
        contains_start_codon = (n_prefix >= n_bases_between_start_and_variant)
        # the sequence context might only include the 2nd or 3rd bases of
        # the start codon
        overlaps_start_codon = (n_prefix > n_bases_between_start_and_variant - 3)

        if contains_start_codon:
            offset_to_first_complete_codon = start_codon_idx - prefix_start_idx
            amino_acids_before_variant = transcript.protein_sequence[:n_full_codons_before_variant]
        else:
            reading_frame = (prefix_start_idx - start_codon_idx) % 3
            offset_to_first_complete_codon = reading_frame_to_offset(reading_frame)
            n_codons_in_prefix = (n_prefix - offset_to_first_complete_codon) // 3
            amino_acids_before_variant = transcript.protein_sequence[
                n_full_codons_before_variant - n_codons_in_prefix:
                n_full_codons_before_variant]

        return cls(
            strand=sequence_key.strand,
            sequence_before_variant_locus=sequence_key.sequence_before_variant_locus,
            sequence_at_variant_locus=sequence_key.sequence_at_variant_locus,
            sequence_after_variant_locus=sequence_key.sequence_after_variant_locus,
            offset_to_first_complete_codon=offset_to_first_complete_codon,
            contains_start_codon=contains_start_codon,
            overlaps_start_codon=overlaps_start_codon,
            contains_five_prime_utr=contains_five_prime_utr,
            amino_acids_before_variant=amino_acids_before_variant)
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)))
396978
def protein_sequence(self, protein_id):
        
        if self.protein_sequences is None:
            raise ValueError(
                "No protein FASTA supplied to this Genome: %s" % self)
        return self.protein_sequences.get(protein_id)
923381
def find_first_stop_codon(nucleotide_sequence):
    
    n_mutant_codons = len(nucleotide_sequence) // 3
    for i in range(n_mutant_codons):
        codon = nucleotide_sequence[3 * i:3 * i + 3]
        if codon in STOP_CODONS:
            return i
    return -1
505726
def reverse_complement(self):
        '''

        '''
        copy = self.copy()
        # Note: if sequence is double-stranded, swapping strand is basically
        # (but not entirely) the same thing - gaps affect accuracy.
        copy.top = self.top.reverse_complement()
        copy.bottom = self.bottom.reverse_complement()

        # Remove all features - the reverse complement isn't flip!
        copy.features = []

        return copy
229669
def InputDNAQuantitySplines(seq_length, n_bases=10, name="DNASmoothPosition", **kwargs):
    
    return Input((seq_length, n_bases), name=name, **kwargs)
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))
1129037
def _create_variant(self,
                        start,
                        end,
                        ref,
                        alt,
                        fsext_len=None,
                        is_dup=False,
                        acc=None,
                        is_ambiguous=False,
                        is_sub=False,
                        is_ext=False,
                        is_no_protein=False,
                        is_init_met=False):
        """"""

        if is_init_met:
            posedit = AARefAlt(ref=ref, alt=alt, init_met=True)
        elif is_ambiguous:
            posedit = None
        else:
            interval = Interval(start=start, end=end)
            # Note - order matters
            if is_no_protein:
                edit = '0'
            elif is_sub:
                edit = AASub(ref=ref, alt=alt)
            elif is_ext:
                edit = AAExt(ref=ref, alt=alt, aaterm='*', length=fsext_len)
            elif self._is_frameshift:
                edit = AAFs(ref=ref, alt=alt, length=fsext_len)
            elif is_dup:
                edit = Dup()
            elif ref == alt == '':
                edit = AARefAlt(ref='', alt='')
            else:
                edit = AARefAlt(ref=ref, alt=alt)
            posedit = PosEdit(
                pos=interval,
                edit=edit,
                uncertain=hgvs.global_config.mapping.inferred_p_is_uncertain)
        var_p = hgvs.sequencevariant.SequenceVariant(acc, 'p', posedit)
        return var_p
1174249
def six_frame(genome, table, minimum = 10):
    
    for seq in parse_fasta(genome):
        dna = Seq(seq[1].upper().replace('U', 'T'), IUPAC.ambiguous_dna)
        counter = 0
        for sequence in ['f', dna], ['rc', dna.reverse_complement()]:
            direction, sequence = sequence
            for frame in range(0, 3):
                for prot in \
                            sequence[frame:].\
                            translate(table = table, to_stop = False).split('*'):
                    if len(prot) < minimum:
                        continue
                    counter += 1
                    header = '%s_%s table=%s frame=%s-%s %s' % \
                                (seq[0].split()[0], counter, table, frame+1, \
                                direction, ' '.join(seq[0].split()[1:]))
                    yield [header, prot]
1074166
def make_degenerate_regex( motif_seq, molecule='dna' ):
	
	
	if not isinstance( motif_seq, six.string_types ):
		raise TypeError( "motif_seq must be a string!" )
		
	if molecule == 'dna':
		degenerate_code = { "A":"A", "B":"[CGT]", "C":"C", "D":"[AGT]", "G":"G", "H":"[ACT]", "K":"[GT]", "M":"[AC]", "N":"[ACGT]", "R":"[AG]", "S":"[GC]", "T":"T", "V":"[ACG]", "W":"[AT]", "Y":"[CT]" }
	elif molecule == 'rna':
		degenerate_code = { "A":"A", "B":"[CGU]", "C":"C", "D":"[AGU]", "G":"G", "H":"[ACU]", "K":"[GU]", "M":"[AC]", "N":"[ACGU]", "R":"[AG]", "S":"[GC]", "U":"U", "V":"[ACG]", "W":"[AU]", "Y":"[CU]" }
	else:
		raise ValueError( "make_degenerate_regex requires molecule to be dna or rna" )
	
	regex_string = ''
	
	idx = 0
	
	while idx < len( motif_seq ):
		curr = motif_seq[idx]
		
		count = 1
		
		for next_idx in range( idx+1, len( motif_seq ) ):
			next = motif_seq[next_idx]
			
			if next == curr:
				count += 1
			else:
				break
		
		regex_string += degenerate_code[curr]
		
		if count > 1:
			idx = idx + count - 1
			regex_string += "{%s}" % ( count )
		
		idx += 1
	
	return regex_string
1235108
public void setParentDNASequence(AbstractSequence<NucleotideCompound> parentDNASequence, Integer begin, Integer end) {
		this.setParentSequence(parentDNASequence);
		setBioBegin(begin);
		setBioEnd(end);
	}
620417
def add_flanking_seqs(self, ref_seq, new_start, new_end):
        ''''''
        if new_start > self.POS or new_end < self.ref_end_pos():
            raise Error('new start and end positions must not try to shrink VCF record. new_start=' + str(new_start) + ', new_end=' + str(new_end) + '. VCF=' + str(self))

        new_start_nucleotides = ref_seq[new_start:self.POS]
        new_end_nucleotodes = ref_seq[self.ref_end_pos() + 1:new_end + 1]
        self.POS = new_start
        self.REF = new_start_nucleotides + self.REF + new_end_nucleotodes
        self.ALT = [new_start_nucleotides + x + new_end_nucleotodes for x in self.ALT]
1229892
public static int countGC(Sequence<NucleotideCompound> sequence) {
		CompoundSet<NucleotideCompound> cs = sequence.getCompoundSet();
		NucleotideCompound G = cs.getCompoundForString("G");
		NucleotideCompound C = cs.getCompoundForString("C");
		NucleotideCompound g = cs.getCompoundForString("g");
		NucleotideCompound c = cs.getCompoundForString("c");
		return countCompounds(sequence, G, C, g, c);
	}
652768
def random_dna(n):
    '''

    '''
    return coral.DNA(''.join([random.choice('ATGC') for i in range(n)]))