#------------------------------------------------------------------------------- # Name: module1 # Purpose: # # Author: oeffner # # Created: 20/05/2013 # Copyright: (c) oeffner 2013 # Licence: #------------------------------------------------------------------------------- from iotbx.bioinformatics import seq_sequence_parse from iotbx import pdb import ccp4io_adaptbx from iotbx import mtz import re, sys, os, math from phaser.utils2 import * def GetSSMSequenceIdentity(pdb1name, pdb2name, chainid1="", chainid2=""): try: pdbxtal = pdb.input(file_name= pdb1name) pdbtest = pdb.input(file_name= pdb2name) hierarchy1 = pdbxtal.construct_hierarchy() hierarchy2 = pdbtest.construct_hierarchy() chain1 = None if chainid1=="": chain1 = hierarchy1.models()[0].chains()[0] else: for chain in hierarchy1.models()[0].chains(): if chainid1[0] == chain.id: chain1 = chain break chain2 = None if chainid2=="": chain2 = hierarchy2.models()[0].chains()[0] else: for chain in hierarchy2.models()[0].chains(): if chainid2[0] == chain.id: chain2 = chain break ssm = ccp4io_adaptbx.SecondaryStructureMatching( moving = chain1, reference = chain2) seqids =[] #count any other alignments so we can iterate and calculate seqid of them Qvalsindex = [(i,e) for (i,e) in enumerate(ssm.qvalues)] # make list with the best SSM Q-score first sortQvals = sorted(Qvalsindex, key= lambda Qvindx: Qvindx[1], reverse=True) for qi in sortQvals: ssm.AlignSelectedMatch(qi[0]) alignment = ccp4io_adaptbx.SSMAlignment.residue_groups( match = ssm ) s = alignment.GetSSMSequenceIdentity() seqids.append((Roundoff(s[0],3), Roundoff(s[1],3), Roundoff(s[2],3), Roundoff(s[3],3), Roundoff(s[4],3), Roundoff(s[5],3), s[6])) except Exception, m: stdoutstr = "SSM alignment of %s and %s failed." %(pdb1name, pdb2name) return [(-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1)] return seqids def GetResolutionFromPDBfile(pdbfname): mypdb = pdb.input(file_name= pdbfname) remark2str = mypdb.extract_remark_iii_records(2)[0] # remark 2 looks like 'REMARK 2 RESOLUTION. 2.80 ANGSTROMS.' match = re.search('REMARK\s+2\s+RESOLUTION\.\s+([-+]?[0-9]*\.?[0-9]+)', remarkstr, re.MULTILINE) resolution = "%3.2f" % float(match.group(1)) return resolution def GetMassOfModelInPDBfile(pdbfname): nO = 0 nC = 0 nN = 0 nS = 0 mypdb = pdb.input(file_name= pdbfname) root = mypdb.construct_hierarchy() # reject any waters present # dryroot = root.select(root.atom_selection_cache().selection( # "not (resname HOH or resname 0WO or resname WAT)"), True) # only count atoms in first model if more models are present for chain in root.models()[0].chains(): for atom in chain.atoms(): if atom.hetero: # ignore HETATM records such as waters or heme groups continue elm = atom.element.strip().upper() if elm == "O": nO = nO + 1 if elm == "C": nC = nC +1 if elm == "N": nN = nN + 1 if elm == "S": nS = nS + 1 mass = nO*16.0 + nC*12.01 + nN*14.01 + nS*32.07 mstr = "Atomic weight of %s is %f" %(pdbfname,mass) return mass, mstr # currently not used def GetMassFromSequence(seq): # copied losely from phaser/lib/scattering.cc seq_error = False seq = seq.upper() H={} C={} N={} O={} S={} #atomic composition of amino acids in protein H['A']=5; C['A']=3; N['A']=1; O['A']=1 H['C']=5; C['C']=3; N['C']=1; O['C']=1; S['C']=1 H['D']=4; C['D']=4; N['D']=1; O['D']=3 H['E']=6; C['E']=5; N['E']=1; O['E']=3 H['F']=9; C['F']=9; N['F']=1; O['F']=1 H['G']=3; C['G']=2; N['G']=1; O['G']=1 H['H']=8; C['H']=6; N['H']=3; O['H']=1 H['I']=11; C['I']=6; N['I']=1; O['I']=1 H['K']=13; C['K']=6; N['K']=2; O['K']=1 H['L']=11; C['L']=6; N['L']=1; O['L']=1 H['M']=9; C['M']=5; N['M']=1; O['M']=1; S['M']=1 H['N']=6; C['N']=4; N['N']=2; O['N']=2 H['P']=7; C['P']=5; N['P']=1; O['P']=1 H['Q']=8; C['Q']=5; N['Q']=2; O['Q']=2 H['R']=13; C['R']=6; N['R']=4; O['R']=1 H['S']=5; C['S']=3; N['S']=1; O['S']=2 H['T']=7; C['T']=4; N['T']=1; O['T']=2 H['V']=9; C['V']=5; N['V']=1; O['V']=1 H['W']=10; C['W']=11; N['W']=2; O['W']=1 H['Y']=9; C['Y']=9; N['Y']=1; O['Y']=2 Hydrogens = 2 # One water for unbonded N- and C-termini Carbons = 0 Nitrogens = 0 Oxygens = 1 Sulphurs = 0 for res in seq: if res not in set(['A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y']): continue Hydrogens += H[res] Carbons += C[res] Nitrogens += N[res] Oxygens += O[res] if res in S: Sulphurs += S[res] weight = Hydrogens + 12.01*Carbons + 14.01*Nitrogens + 16.0*Oxygens + 32.07*Sulphurs return weight def GetNumberofResiduesinPDBfile(pdbfname): nres = 0 mypdb = pdb.input(file_name= pdbfname) root = mypdb.construct_hierarchy() for chain in root.models()[0].chains(): for rg in chain.residue_groups(): # don't count HETATOMs or waters if rg.atoms()[0].parent().resname in pdb.common_residue_names_amino_acid: nres += 1 return nres def GetNumberofResiduesFromSequence(sequence, raw=False): if not raw: seq = seq_sequence_parse(sequence)[0][0].sequence else: seq = sequence return len(seq) def GetScatteringOfModelInPDBfile(pdbfname): scattering = 0 mypdb = pdb.input(file_name= pdbfname) root = mypdb.construct_hierarchy() # reject any waters present # dryroot = root.select(root.atom_selection_cache().selection( # "not (resname HOH or resname 0WO or resname WAT)"), True) # only count atoms in first model if more models are present # for chain in dryroot.models()[0].chains(): for chain in root.models()[0].chains(): nO = 0 nC = 0 nN = 0 nS = 0 nH = 0 for atom in chain.atoms(): if atom.hetero: # ignore HETATM records such as waters or heme groups continue elm = atom.element.strip().upper() if elm == "H": nH += atom.occ if elm == "O": nO += atom.occ if elm == "C": nC += atom.occ if elm == "N": nN += atom.occ if elm == "S": nS += atom.occ # scattering is proportional to the square of the numbers of electrons as we are dealing with intensities scattering += (nH + nO*8.0*8.0 + nC*6.0*6.0 + nN*7.0*7.0 + nS*16.0*16.0) mstr = "Atomic scattering of %s is %f" %(pdbfname, scattering) return scattering, mstr def GetScatteringFromSequence(sequence, raw=False): # copied losely from phaser/lib/scattering.cc if not raw: seq = seq_sequence_parse(sequence)[0][0].sequence else: seq = sequence seq_error = False seq = seq.upper() H={} C={} N={} O={} S={} #atomic composition of amino acids in protein H['A']=5; C['A']=3; N['A']=1; O['A']=1 H['C']=5; C['C']=3; N['C']=1; O['C']=1; S['C']=1 H['D']=4; C['D']=4; N['D']=1; O['D']=3 H['E']=6; C['E']=5; N['E']=1; O['E']=3 H['F']=9; C['F']=9; N['F']=1; O['F']=1 H['G']=3; C['G']=2; N['G']=1; O['G']=1 H['H']=8; C['H']=6; N['H']=3; O['H']=1 H['I']=11; C['I']=6; N['I']=1; O['I']=1 H['K']=13; C['K']=6; N['K']=2; O['K']=1 H['L']=11; C['L']=6; N['L']=1; O['L']=1 H['M']=9; C['M']=5; N['M']=1; O['M']=1; S['M']=1 H['N']=6; C['N']=4; N['N']=2; O['N']=2 H['P']=7; C['P']=5; N['P']=1; O['P']=1 H['Q']=8; C['Q']=5; N['Q']=2; O['Q']=2 H['R']=13; C['R']=6; N['R']=4; O['R']=1 H['S']=5; C['S']=3; N['S']=1; O['S']=2 H['T']=7; C['T']=4; N['T']=1; O['T']=2 H['V']=9; C['V']=5; N['V']=1; O['V']=1 H['W']=10; C['W']=11; N['W']=2; O['W']=1 H['Y']=9; C['Y']=9; N['Y']=1; O['Y']=2 # Hydrogens = 2 # One water for unbonded N- and C-termini Carbons = 0 Nitrogens = 0 Oxygens = 1 Sulphurs = 0 for res in seq: if res not in set(['A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y']): continue Hydrogens += H[res] Carbons += C[res] Nitrogens += N[res] Oxygens += O[res] if res in S: Sulphurs += S[res] # scattering is proportional to the square of the numbers of electrons as we are dealing with intensities scattering = Hydrogens + 6.0*6.0*Carbons + 7.0*7.0*Nitrogens + 8.0*8.0*Oxygens + 16.0*16.0*Sulphurs mstr = "Scattering of sequence is %f" %scattering return scattering, mstr def GetMtzColumns(fname): #get columns labels from mtz file Fcol = ""; Sigcol = ""; Icol = ""; Rfreecol = "" avgF = 0.0; avgSig = 0.0; avgI = 0.0 exceptmsg = "" try: obj = mtz.object(file_name = fname) columns = obj.column_labels() for i in range(0,len(columns)): columntype = obj.get_column(columns[i]).type() values = obj.get_column(columns[i]).extract_values() if (Fcol == ''): # to be assigned below if (columntype=='F'): Fcol = columns[i] #if min(values) > 0.0: avgF = math.fsum(values)/len(values) if (Sigcol == ""): if (columntype=='Q'): Sigcol = columns[i] #if min(values) > 0.0: avgSig = math.fsum(values)/len(values) if (Icol == ""): if (columntype=='J'): Icol = columns[i] if min(values) > 0.0: avgI = math.fsum(values)/len(values) if (Rfreecol == ""): if (columntype=='I'): Rfreecol = columns[i] except Exception, m: exceptmsg = m return Fcol, Sigcol, Icol, Rfreecol, avgF, avgSig, avgI, exceptmsg def mymain(args): if len(args) < 3: execmd= "phenix.python " + os.path.basename(args[0]) print "usage: %s