pKa/get_muts.py
def main(options,args):
    """Extract mutations from a multiple sequence alignment"""
    import PEATDB.sequence_alignment as SA
    if not options.fasta:
        alignment=SA.read_clustal_aln_file(args[0])
    else:
        alignment=SA.read_fasta_aln_file(args[0])
    print sorted(alignment.keys())

    HEWL_seq=alignment[options.wt]


    fd=open('mutations.txt','w')
    fd2=open('frequencies.csv','w')

    aas=convert.keys()
    aas.sort()
    import string
    fd2.write('WT Residue number, %s\n' %(string.join(aas,',')))
    #
    real_pos=0
    lines = []
    PDB_mutations={}
    #
    for position in range(0,len(HEWL_seq)):
        res_observed={}
        if HEWL_seq[position]=='-':
            continue
        real_pos=real_pos+1
        #print 'Now looking at position',real_pos
        for seq in alignment.keys():
            res_pos=alignment[seq][position]
            if not res_observed.has_key(res_pos):
                res_observed[res_pos]=0
            res_observed[res_pos]=res_observed[res_pos]+1
        #
        # Calculate frequencies of observation
        #
        total=sum(res_observed.values())
        text='%3d' %real_pos
        #print res_observed.keys()
        for aa in aas:
            if res_observed.has_key(aa):
                text=text+',%5.1f' %(float(res_observed[aa])/total*100.0)
            else:
                text=text+',%5.1f' %(0)
        fd2.write(text+'\n')
        #
        # -----
        #
        lines += ['%3d   %d\n' %(real_pos,len(res_observed.keys()))]
        for mut in res_observed.keys():
            if mut=='-':
                continue
            if mut==HEWL_seq[position]:
                continue
            #
            org_res=HEWL_seq[position]
            new_res=mut
            import string
            if org_res=='X' or new_res=='X':
                continue
            #
            # Within the PDB file?
            #
            if real_pos<options.start_aa or real_pos>options.end_aa:
                pass
            else:
                PDB_residue='%s:%s' %(options.CID,string.zfill(real_pos+options.offset-options.noffset,4))
                if not PDB_mutations.has_key(PDB_residue):
                    PDB_mutations[PDB_residue]=[]
                PDB_mutations[PDB_residue].append(convert[new_res])
                muttext='%s:%s:%s:%s' %(options.CID,string.zfill(real_pos+options.offset-options.noffset,4),convert[org_res],convert[new_res])
                fd.write('%s,%s\n' %(muttext,muttext))
                #print muttext
    fd.close()
    fd2.close()
    #
    # Read PDB file?
    #
    if options.pdbfile:
        import Protool
        PI=Protool.structureIO()
        PI.readpdb(options.pdbfile)
    #
    # Plot the figure?
    #
    if options.plotfigure:
        xs=[]
        ys=[]
        zs=[]
        for residue in sorted(PDB_mutations.keys()):
            resnum=int(residue.split(':')[1])
            xs.append(resnum)
            ys.append(len(PDB_mutations[residue]))
            zs.append(PI.dist(options.atom,residue+':CA'))
        import pylab
        pylab.plot(zs,ys,'ro')
        pylab.xlabel('Distance from %s' %(options.atom))
        pylab.ylabel('Number of mutations')
        pylab.show()
    return
Extract mutations from a multiple sequence alignment
Extract mutations from a multiple sequence alignment
Extract mutations from a multiple sequence alignment
def main(options,args):
    import PEATDB.sequence_alignment as SA
    if not options.fasta:
        alignment=SA.read_clustal_aln_file(args[0])
    else:
        alignment=SA.read_fasta_aln_file(args[0])
    print sorted(alignment.keys())
    HEWL_seq=alignment[options.wt]
    fd=open('mutations.txt','w')
    fd2=open('frequencies.csv','w')
    aas=convert.keys()
    aas.sort()
    import string
    fd2.write('WT Residue number, %s\n' %(string.join(aas,',')))
    real_pos=0
    lines = []
    PDB_mutations={}
    for position in range(0,len(HEWL_seq)):
        res_observed={}
        if HEWL_seq[position]=='-':
            continue
        real_pos=real_pos+1
        print 'Now looking at position',real_pos
        for seq in alignment.keys():
            res_pos=alignment[seq][position]
            if not res_observed.has_key(res_pos):
                res_observed[res_pos]=0
            res_observed[res_pos]=res_observed[res_pos]+1
         Calculate frequencies of observation
        total=sum(res_observed.values())
        text='%3d' %real_pos
        print res_observed.keys()
        for aa in aas:
            if res_observed.has_key(aa):
                text=text+',%5.1f' %(float(res_observed[aa])/total*100.0)
            else:
                text=text+',%5.1f' %(0)
        fd2.write(text+'\n')
         -----
        lines += ['%3d   %d\n' %(real_pos,len(res_observed.keys()))]
        for mut in res_observed.keys():
            if mut=='-':
                continue
            if mut==HEWL_seq[position]:
                continue
            org_res=HEWL_seq[position]
            new_res=mut
            import string
            if org_res=='X' or new_res=='X':
                continue
             Within the PDB file?
            if real_pos<options.start_aa or real_pos>options.end_aa:
                pass
            else:
                PDB_residue='%s:%s' %(options.CID,string.zfill(real_pos+options.offset-options.noffset,4))
                if not PDB_mutations.has_key(PDB_residue):
                    PDB_mutations[PDB_residue]=[]
                PDB_mutations[PDB_residue].append(convert[new_res])
                muttext='%s:%s:%s:%s' %(options.CID,string.zfill(real_pos+options.offset-options.noffset,4),convert[org_res],convert[new_res])
                fd.write('%s,%s\n' %(muttext,muttext))
                print muttext
    fd.close()
    fd2.close()
     Read PDB file?
    if options.pdbfile:
        import Protool
        PI=Protool.structureIO()
        PI.readpdb(options.pdbfile)
     Plot the figure?
    if options.plotfigure:
        xs=[]
        ys=[]
        zs=[]
        for residue in sorted(PDB_mutations.keys()):
            resnum=int(residue.split(':')[1])
            xs.append(resnum)
            ys.append(len(PDB_mutations[residue]))
            zs.append(PI.dist(options.atom,residue+':CA'))
        import pylab
        pylab.plot(zs,ys,'ro')
        pylab.xlabel('Distance from %s' %(options.atom))
        pylab.ylabel('Number of mutations')
        pylab.show()
    return
def main ( options , args ) : import PEATDB . sequence_alignment as SA if not options . fasta : alignment = SA . read_clustal_aln_file ( args [ 0 ] ) else : alignment = SA . read_fasta_aln_file ( args [ 0 ] ) print sorted ( alignment . keys ( ) ) HEWL_seq = alignment [ options . wt ] fd = open ( 'mutations.txt' , 'w' ) fd2 = open ( 'frequencies.csv' , 'w' ) aas = convert . keys ( ) aas . sort ( ) import string fd2 . write ( 'WT Residue number, %s\n' % ( string . join ( aas , ',' ) ) ) real_pos = 0 lines = [ ] PDB_mutations = { } for position in range ( 0 , len ( HEWL_seq ) ) : res_observed = { } if HEWL_seq [ position ] == '-' : continue real_pos = real_pos + 1 for seq in alignment . keys ( ) : res_pos = alignment [ seq ] [ position ] if not res_observed . has_key ( res_pos ) : res_observed [ res_pos ] = 0 res_observed [ res_pos ] = res_observed [ res_pos ] + 1 total = sum ( res_observed . values ( ) ) text = '%3d' % real_pos for aa in aas : if res_observed . has_key ( aa ) : text = text + ',%5.1f' % ( float ( res_observed [ aa ] ) / total * 100.0 ) else : text = text + ',%5.1f' % ( 0 ) fd2 . write ( text + '\n' ) lines += [ '%3d   %d\n' % ( real_pos , len ( res_observed . keys ( ) ) ) ] for mut in res_observed . keys ( ) : if mut == '-' : continue if mut == HEWL_seq [ position ] : continue org_res = HEWL_seq [ position ] new_res = mut import string if org_res == 'X' or new_res == 'X' : continue if real_pos < options . start_aa or real_pos > options . end_aa : pass else : PDB_residue = '%s:%s' % ( options . CID , string . zfill ( real_pos + options . offset - options . noffset , 4 ) ) if not PDB_mutations . has_key ( PDB_residue ) : PDB_mutations [ PDB_residue ] = [ ] PDB_mutations [ PDB_residue ] . append ( convert [ new_res ] ) muttext = '%s:%s:%s:%s' % ( options . CID , string . zfill ( real_pos + options . offset - options . noffset , 4 ) , convert [ org_res ] , convert [ new_res ] ) fd . write ( '%s,%s\n' % ( muttext , muttext ) ) fd . close ( ) fd2 . close ( ) if options . pdbfile : import Protool PI = Protool . structureIO ( ) PI . readpdb ( options . pdbfile ) if options . plotfigure : xs = [ ] ys = [ ] zs = [ ] for residue in sorted ( PDB_mutations . keys ( ) ) : resnum = int ( residue . split ( ':' ) [ 1 ] ) xs . append ( resnum ) ys . append ( len ( PDB_mutations [ residue ] ) ) zs . append ( PI . dist ( options . atom , residue + ':CA' ) ) import pylab pylab . plot ( zs , ys , 'ro' ) pylab . xlabel ( 'Distance from %s' % ( options . atom ) ) pylab . ylabel ( 'Number of mutations' ) pylab . show ( ) return
Extract mutations from a multiple sequence alignment
Extract mutations from a multiple sequence alignment
"""Extract mutations from a multiple sequence alignment""" # # #print 'Now looking at position',real_pos # # Calculate frequencies of observation # #print res_observed.keys() # # ----- # # # # Within the PDB file? # #print muttext # # Read PDB file? # # # Plot the figure? #
options
args