Python Bio.Seq类代码示例

def cutOligos(GeneName, cutsite, DNA):
    
# This function asks user for a target cut site, and generates oligos to clone that cut site into 
# a CRISPR plasmid, and to screen for positive clones containing that cut site. It may also generate
# universal homology regions to integrate or otherwise alter that target region. 
#
# The current cut site architecture is:  >>>>> YtRNAp-HDV ribozyme- >20nt< -gRNA <<<<<

    GeneName=input("Name, using quotes: ")
    cutsite=input("20-mer cut sequence, using quotes: ")
    DNA=input("Locus sequence +/- a few kb, using quotes: ")

    if DNA.find(cutsite)==-1:                   # If cutiste sequence found in ANTISENSE
        DNA=Seq(DNA).reverse_complement()       # then reverse DNA, and turn it into a string
    
    index=DNA.find(cutsite)+16                  # index gives the start position of the string, e.g., 0. 
                                                # we add 16 since index+0=start of 20-mer, so index+16=cut site, 
                                                # 3 nt before last of 20mer
        
    Lup=DNA[index-520:index-490]                # This primer binds 500bp upstream of cut site

    
    cutSequence=Seq("cgggtggcgaatgggacttt")+cutsite+Seq("gttttagagctagaaatagc")
    seqprimer=Seq("gacttt")+cutsite
    

    print("cut" + GeneName + "  " + cutSequence)
    print("Lcolony" + GeneName + "  " + seqprimer)
    print("Lup" + GeneName + "  " + Lup)

def ReadingFrameFinder(DNASTRING):
    CleanDNA = DNASTRING.rstrip("\n")
    OpenLocations = []
    CloseLocations = []
    stringlen = len(CleanDNA)
    TtoU = CleanDNA.replace("T", 'U')
    readingframeRange = xrange(0, stringlen)
    PossibleGenes = []
    for item in readingframeRange:
        if TtoU[item:item+3] == "AUG":
            Newthing = xrange(item, stringlen, 3)
            storage = item
            for number in Newthing:
                if TtoU[number:number+3] == "UAA" or TtoU[number:number+3] == "UAG" or TtoU[number:number+3] == "UGA":
                    PossibleGenes.append(TtoU[storage:number+3])
                    break
    for Seqeu in PossibleGenes:
        if len(Seqeu) % 3 == 0:
            LETGO = Seq(Seqeu, generic_rna)
            FinalizedProt.append(str(LETGO.translate()))
        else:
            Removal_Len = len(Seqeu) % 3
            UpdatedSequence = Seqeu[:-Removal_Len]
            ETGO2 = Seq(UpdatedSequence, generic_rna)
            FinalizedProt.append(str(ETGO2.translate()))

def assign_fitness(nodes):
	'''
	loops over all viruses, translates their sequences and calculates the virus fitness
	'''
	aa, sites, wt_aa, aa_prob = load_mutational_tolerance()
	aln = AlignIO.read('source-data/H1_H3.fasta', 'fasta')
	# returns true whenever either of the sequences have a gap
	aligned = (np.array(aln)!='-').min(axis=0)
	# map alignment positions to sequence positions, subset to aligned amino acids
	indices = {}
	for seq in aln:
		indices[seq.name] = (np.cumsum(np.fromstring(str(seq.seq), dtype='S1')!='-')-1)[aligned]

	# make a reduced set of amino-acid probabilities that only contains aligned positions
	aa_prob=aa_prob[indices['H1'],:]
	# attach another column for non-canonical amino acids
	aa_prob = np.hstack((aa_prob, 1e-5*np.ones((aa_prob.shape[0],1))))
	if isinstance(nodes, list):
		for node in nodes:
			node['tol'] = calc_fitness_tolerance(Seq.translate(node['seq']), 
															aa_prob, aa, indices['H3'])
	elif isinstance(nodes, dendropy.Tree):
		for node in nodes.postorder_node_iter():
			node.tol = calc_fitness_tolerance(Seq.translate(node.seq), 
															aa_prob, aa, indices['H3'])

    def test_reverse_complement_on_proteins(self):
        """Test reverse complement shouldn't work on a protein!"""
        for s in protein_seqs:
            with self.assertRaises(ValueError):
                Seq.reverse_complement(s)

            with self.assertRaises(ValueError):
                s.reverse_complement()

    def test_translation_to_stop(self):
        for nucleotide_seq in self.test_seqs:
            nucleotide_seq = nucleotide_seq[:3 * (len(nucleotide_seq) // 3)]
            if isinstance(nucleotide_seq, Seq.Seq) and 'X' not in str(nucleotide_seq):
                short = Seq.translate(nucleotide_seq, to_stop=True)
                self.assertEqual(str(short), str(Seq.translate(nucleotide_seq).split('*')[0]))

        seq = "GTGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"
        self.assertEqual("VAIVMGRWKGAR", Seq.translate(seq, table=2, to_stop=True))

    def test_back_transcription_of_proteins(self):
        """Test back-transcription shouldn't work on a protein!"""
        for s in protein_seqs:
            with self.assertRaises(ValueError):
                Seq.back_transcribe(s)

            if isinstance(s, Seq.Seq):
                with self.assertRaises(ValueError):
                    s.back_transcribe()

    def test_translation_on_proteins(self):
        """Test translation shouldn't work on a protein!"""
        for s in protein_seqs:
            with self.assertRaises(ValueError):
                Seq.translate(s)

            if isinstance(s, Seq.Seq):
                with self.assertRaises(ValueError):
                    s.translate()

 def rc_kmers(self, kmers):
     res={}
     keys=[]
     for s in kmers:
         if Seq.reverse_complement(s) in keys:
             res[s]=Seq.reverse_complement(s)
         else:
             keys.append(s)
             res[s]=s
     return keys,res

def delGene(geneName, cutsite):

# This function asks user for a chromosomal locus, a region to be deleted, a suitable CRIPSR cutsite 
# and outputs oligos for cloning of a pL308 Cas9-gRNA vector, and ones for generating a donor DNA
# to delete the unwanted chromosomal region. Primers Lup+Rdown produce a 1kb band if deletion was
# successful. 
# part of yCRISPRv3 by [email protected]

    #GeneName=input("Name, using quotes: ")
    #cutsite=input("20-mer cut sequence, using quotes: ").upper()
    locus = genomicData[geneName][0]    
    deletion = genomicData[geneName][1]

    deletion = Seq(deletion)
     
    if deletion.find(cutsite)==-1:
        if deletion.reverse_complement().find(cutsite)==-1:
            print ("WARNING: Guide 20-mer sequence not found in deletion region.")
        
    locus=Seq(locus)
    
    index=locus.find(deletion)                  
   
    # index gives the start position within locus of the string deletion.
    # now we delete the deletion region to redefine a newlocus:
    
    newlocus=locus[0:index]+locus[index+len(deletion):]

    # note that since index starts at 0, a value of n points to, in the newlocus,
    # the first nt after the deletion. So we define the newlocus as above. Note too
    # that a string of len=40 ends at an index of 39--so we pick up at index+len-1. 
    
    Lup=newlocus[index-500:index-470]
    Rdown=newlocus[index+469:index+499].reverse_complement()

    Rtemp1 = newlocus[:index].reverse_complement()
    Rtemp2 = newlocus[index:].reverse_complement()

    rPrimer, rLength = getPrimer(Rtemp1)
    lPrimer, lLength = getPrimer(newlocus[index:])

    Rup = getOverhang(Rtemp2, rLength) + rPrimer
    Ldown = getOverhang(newlocus[:index], lLength) + lPrimer

    cutSequence=Seq("cgggtggcgaatgggacttt")+cutsite+Seq("gttttagagctagaaatagc")
    seqprimer=Seq("gacttt")+cutsite
    
    print("cut" + GeneName + "  " + cutSequence)
    print("seq" + GeneName + "  " + seqprimer)
    print("Lup" + GeneName + "del" + " " + Lup)
    print("Rup" + GeneName + "del" + " " + Rup)
    print("Ldown" + GeneName + "del" + " " + Ldown)
    print("Rdown" + GeneName + "del" + " " + Rdown)

    return Ldown, Rup

def calc_total_subst(start_codon, end_codon):
    """
    Returns total synonymous substitutions, nonsynonymous substitutions.
    If there are multiple positions that differ between codons, then returns the average synonynous substitutions,
    average nonsynonymous substitutions across all possible pathways from codon1 to codon2
    where each stage in a pathway is separated by 1 position mutation.
    :param Bio.Seq.Seq start_codon:  3bp codon
    :param Bio.Seq.Seq end_codon:  3bp codon
    :return tuple (int, int):  (average point mutations that yield same amino acid across all pathways, average point mutations that yield different amino acid across all pathways)
    """
    total_syn = 0.0
    total_nonsyn = 0.0
    total_subs = 0.0

    upper_start_codon = start_codon.upper()
    upper_end_codon = end_codon.upper()

    # find positions where the codons differ
    diff_pos = []
    for pos, nucstr1 in enumerate(str(upper_start_codon)):
        nucstr2 = str(upper_end_codon[pos])
        if nucstr1 != nucstr2:
            diff_pos.extend([pos])

    # Traverse all possible pathways from start_codon to end_codon where
    # each stage of a pathway mutates by 1 base.
    last_codon = upper_start_codon
    last_aa = Seq.translate(last_codon)
    for pathway in itertools.permutations(diff_pos):
        print str(upper_start_codon) + " " + str(upper_end_codon) + " " + ",".join([str(x) for x in pathway])
        for mut_pos in pathway:
            mut_nuc = upper_end_codon[mut_pos]
            mut_codon =  last_codon[:mut_pos] + mut_nuc + last_codon[mut_pos+1:]
            mut_aa = Seq.translate(mut_codon)

            total_subs += 1
            if str(last_aa) == str(mut_aa):
                total_syn += 1
            else:
                total_nonsyn += 1

            last_codon = mut_codon
            last_aa = mut_aa

        if str(last_codon) != str(upper_end_codon):
            raise ValueError("Pathway does not yield end codon " + str(last_codon))

    if total_subs:
        ave_syn = total_syn/total_subs
        ave_nonsyn = total_nonsyn/total_subs
    else:
        ave_syn = 0.0
        ave_nonsyn = 0.0
    return ave_syn, ave_nonsyn

    def test_reverse_complement(self):
        test_seqs_copy = copy.copy(test_seqs)
        test_seqs_copy.pop(21)

        for nucleotide_seq in test_seqs_copy:
            if not isinstance(nucleotide_seq.alphabet, Alphabet.ProteinAlphabet) and \
                    isinstance(nucleotide_seq, Seq.Seq):
                expected = Seq.reverse_complement(nucleotide_seq)
                self.assertEqual(repr(expected), repr(nucleotide_seq.reverse_complement()))
                self.assertEqual(repr(expected[::-1]), repr(nucleotide_seq.complement()))
                self.assertEqual(str(nucleotide_seq.complement()),
                                 str(Seq.reverse_complement(nucleotide_seq))[::-1])
                self.assertEqual(str(nucleotide_seq.reverse_complement()),
                                 str(Seq.reverse_complement(nucleotide_seq)))

def translate(config, rc=False):
    table = 1
    if mycoplasma(config):
        # table 4 is for mycoplasma ala:
        # http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi
        table = 4
    fd, fmap = None, None
    try:
        log.debug("Doing translation with table %d, rc: %s", table, rc)
        fd = os.open(ddna(config), os.O_RDONLY)
        fmap = mmap.mmap(fd, 0, mmap.MAP_SHARED, mmap.PROT_READ)
        # By convention (e.g. from the C or NCBI) the DNA is is 1
        # indexed; our DDNA is a c style array that is 0 indexed
        startIdx = config['startBase'] - 1
        # The end index here is inclusive but array.slice isn't so we
        # don't need to subtract 1
        endIdx = config['endBase']
        seq = Seq.Seq(fmap[startIdx:endIdx])

        if rc:
            seq = seq.reverse_complement()
        return {
            'seq': str(seq),
            'trans': str(Seq.translate(seq, table))
        }
    finally:
        if fmap:
            fmap.close
        if fd:
            os.close(fd)

    def export(self, path = '', extra_attr = ['aa_muts']):
        from Bio import Seq
        from itertools import izip
        timetree_fname = path+'tree.json'
        sequence_fname = path+'sequences.json'
        tree_json = tree_to_json(self.tree.root, extra_attr=extra_attr)
        write_json(tree_json, timetree_fname, indent=None)
        elems = {}
        elems['root'] = {}
        elems['root']['nuc'] = "".join(self.tree.root.sequence)
        for prot in self.proteins:
            tmp = str(self.proteins[prot].extract(Seq.Seq(elems['root']['nuc'])))
            #elems['root'][prot] = str(Seq.translate(tmp.replace('---', 'NNN'))).replace('X','-')
            elems['root'][prot] = str(Seq.translate(tmp.replace('-', 'N'))).replace('X','-')


        for node in self.tree.find_clades():
            if hasattr(node, "clade") and hasattr(node, "sequence"):
                elems[node.clade] = {}
                elems[node.clade]['nuc'] = {pos:state for pos, (state, ancstate) in
                                enumerate(izip(node.sequence, self.tree.root.sequence)) if state!=ancstate}
        for node in self.tree.find_clades():
            if hasattr(node, "clade") and hasattr(node, "translations"):
                for prot in self.proteins:
                    elems[node.clade][prot] = {pos:state for pos, (state, ancstate) in
                                    enumerate(izip(node.translations[prot], elems['root'][prot])) if state!=ancstate}

        write_json(elems, sequence_fname, indent=None)

def translateDNAtoAA(input_fasta, output_fasta, remove_lower_case = False):
    with open(input_fasta, 'r') as f:
        with open(output_fasta, 'w+') as g:
            for line in f.readlines():
                if line[0] == '>':
                    g.write(line)
                    continue
                else:
                    if line[-2:] == '\r\n':
                        assert(len(line) %3 == 2)
                    elif line[-1:] == '\n':
                        assert(len(line) %3 == 1)
                    if remove_lower_case:
                        g.write(Seq.translate(line.translate(None, string.ascii_lowercase)[:-1], to_stop = True) + '\n')
                    else:
                        g.write(Seq.translate(line[:-1], to_stop = True) + '\n')

def translationBio(data):
    '''Uses Biopython translate '''
    proteinSeq = ''
    for line in data:
        proteinSeq += Seq.translate(line, table='Standard', stop_symbol='', to_stop=False)
        #proteinSeq += Seq.translate(line)
    print proteinSeq