<commit_before>import tempfile, os, sys, commands, re, math

"""
Routines to rename atoms and residues from a 'labeled' step2_out.pdb file  ("0","+", or "-" tacked on the end) 

REVISION LOG
- 8-21-2007: VAV; spun these functions into their own python file 
"""

def MCCELabeledPDB_to_Amber(labeledPDBfile, outPDBfile):
    """Reads a 'labeled' PDB file written by MCCE (with "0","+", or "-" indicating the protonation state - see mcce.py)
    Example:
    ATOM      1  CA  LEU _0001_000  -8.847   7.195  16.727   1.700       0.000      BK____M000
    0
    ATOM      2  C   LEU _0001_000  -9.019   7.821  15.386   1.700       0.000      BK____M000
    0
    ....

    Writes an AMBER-residue-formatted PDB file, outPDBfile.
    """

    # Read in the pdb lines
    fin = open(labeledPDBfile,'r')
    lines = fin.readlines()
    fin.close()
    # Format so the even-numbered lines, "0", "+", "-", etc are supposed to be tacked onto the ends of each atom's line
    for i in range(0, len(lines)):
        asdfghjk
 
    pdbarr=rename.rename_residues(pdbarr)
    pdbarr=rename.pdb_cleanup(pdbarr)
    return pdbarr


def renumber_atoms(pdbarr):
    """Renumbers atom entries (list items) of a PDB file so they are in sequence starting from 1.
    
    ARGUMENTS
        pdbarr - list of PDB atoms to be renumbered; mutated by call
    """
    
    # Renumber the atom entries in place.
    for i in range(len(pdbarr)):
        pdbarr[i]=pdbarr[i][0:6]+str(i+1).rjust(5)+pdbarr[i][11:]        

def rename_termini(npdb):
    """Renames the 'NTR' and 'CTR' residues generated by MCCE.
    Renames the O and OXT atoms at the C-terminus to OC1 and OC2 for AMBER.
    
    ARGUMENTS
        npdb - 
    
    """

    # Get three-letter residue names of N- and C-termini.
    ntrname = npdb[1][0][17:20]
    ctrname = npdb[-2][0][17:20]

    # Rename N-terminal residue 'XXX' to 'NXXX'.
    for j in [0,1]:
        for i in range(len(npdb[j])):
            npdb[j][i] = npdb[j][i][0:17]+'N'+ntrname+npdb[j][i][21:]
            
    # Rename C-terminal residue from 'XXX' to 'CXXX', with some exceptions (noted below).
    # Also rename terminal oxygen atoms.
    for j in [-1,-2]:
        for i in range(len(npdb[j])):
            # VAV 5/14/2007: gmx ffamber does not have CLYN, only CLYP 
            if ctrname=='LYN':      
                npdb[j][i] = npdb[j][i][0:17]+'CLYP'+npdb[j][i][21:]
            else:   
                npdb[j][i] = npdb[j][i][0:17]+'C'+ctrname+npdb[j][i][21:]
                            
            # DLM 5/9/2007: Though this says that it changes the O and OXT names, it doesn't
            # Adding the below to do the renaming.
            if npdb[j][i][13:16].split()[0]=='O':
                npdb[j][i] = npdb[j][i][0:13]+'OC1'+npdb[j][i][16:]
            elif npdb[j][i][13:16].split()[0]=='OXT':
                npdb[j][i] = npdb[j][i][0:13]+'OC2'+npdb[j][i][16:]
                
            #VAV 6/5/2007: gmx ffamber doesn't have "CD1" as an ILE/CILE atomname, only "CD"
            if ctrname=='ILE':      
                if npdb[j][i][13:16].split()[0]=='CD1':
                    npdb[j][i] = npdb[j][i][0:13]+'CD '+npdb[j][i][16:]
        
def nest_pdb(pdbarr):
    """Collect lines of the PDB file by residue.
    
    ARGUMENTS:
        pdbarr - list of lines from PDB file
        
    RETURNS
        nestedpdb - nested PDB file, such that nestedpdb[i][j] will be line j from residue i.
    """
    
    nestedpdb=[]
    residue=[]
    for line in pdbarr:
        if (len(residue) == 0):
            residue.append(line)
        else:
            if (line[17:27] == residue[-1][17:27]):
                residue.append(line)
            else:
                nestedpdb.append(residue)
                residue=[line]
    nestedpdb.append(residue)
    
    # Return the nexted PDB file.
    return nestedpdb

def unnest_pdb(npdb):
    """Expand the lines from the "nested" PDB file into a flat list of lines.
    
    ARGUMENTS
        npdb - "nested" PDB file; mutated in place
        
    RETURNS
        pdbarr - un-nested PDB file, such that pdbarr[i] will be line i from the PDB file
    """
    
    pdbarr=[]
    for res in npdb:
        for atm in res:
            pdbarr.append(atm)
    return pdbarr
    
def rename_residues(pdbarr):
    """Convert residue names (and some heavy atom names) to AMBER (ffamber) format.
    
    ARGUMENTS
        pdbarr - array of PDB atoms
        
    RETURNS
        pdbarr - updated array of PDB atoms with AMBER appropriate residue names
    """

    # DEBUG check.
    #print "At start of renaming, pdbarr has",len(pdbarr),"items"

    # Transform the PDB file into a nested format that is easier to manipulate.
    npdb = nest_pdb(pdbarr)
    
    # DEBUG check.
    #print "Nest/unnest leads to",len(unnest_pdb(npdb)),"items"

    # Convert residue and atom names to be appropriate for AMBER.
    rename_charged(npdb)
    histidine_search(npdb)
    disulfide_search(npdb)
    rename_termini(npdb)

    # Un-nest PDB
    pdbarr = unnest_pdb(npdb)
    
    # DEBUG check.
    #print "At end of renaming, pdbarr has ", len(pdbarr), "items"

    # Return the PDB file.
    return pdbarr

def rename_charged(npdb):
    """Generate AMBER-specific residue names for charged residues from MCCE residue names.
        
    ARGUMENTS
        npdb - "nested" PDB data structure generated by nest_pdb. Modified to reflect AMBER names.
    """
        
    resname_and_state = map(lambda x:x[-1][17:20]+x[-1][-1],npdb)
    #print resname_and_state
    for i in range(len(npdb)):
        if (resname_and_state[i]=='HIS+'):
            npdb[i]=map(lambda x:x.replace('HIS','HIP'),npdb[i])
        elif (resname_and_state[i]=='LYS0'):
            npdb[i]=map(lambda x:x.replace('LYS','LYN'),npdb[i])
        elif (resname_and_state[i]=='LYS+'):
            npdb[i]=map(lambda x:x.replace('LYS','LYP'),npdb[i])
        elif (resname_and_state[i]=='CYS0'):
            npdb[i]=map(lambda x:x.replace('CYS','CYN'),npdb[i])
        elif (resname_and_state[i]=='CYS-'):
            npdb[i]=map(lambda x:x.replace('CYS','CYM'),npdb[i])
        elif (resname_and_state[i]=='ASP0'):
            npdb[i]=map(lambda x:x.replace('ASP','ASH'),npdb[i])
        #Aspartate requires no sub
        elif (resname_and_state[i]=='GLU0'):
            npdb[i]=map(lambda x:x.replace('GLU','GLH'),npdb[i])
        # Glutamate requires no sub

    return

def get_coords(atomname,residue):
    """Extract coordinates from given atom in list of PDB lines (such as all those lines for a residue in a nested PDB file).
    
    ARGUMENTS
        atomname - the name of the atom
        residue - list of lines from PDB file

    RETURNS
        (x,y,z) - tuple of atom coordinates, as floats
    
    """
    for n in range(len(residue)):
        if (residue[n][12:16].strip() == atomname.strip()):
            iX=float(residue[n][30:38])
            iY=float(residue[n][38:46])
            iZ=float(residue[n][46:54])
            return (iX,iY,iZ)

    raise "Atom not found!"        

def disulfide_search(npdb, min_dist = 1.8, max_dist = 2.2):
    """Rename CYS to CYX if it participates in a disulfide bond, as judged by distance range (inclusive).
    
    ARGUMENTS
        npdb - nested PDB file
        
    OPTIONAL ARGUMENTS
        min_dist - minium distance cutoff for perceiving disulfide bond (default 1.8 A)
        max_dist - maximum distance cutoff for perceiving disulfide bond (default 2.2 A)
    """
        
    residues_to_rename=set([])
    for i in range(len(npdb)):
        if (npdb[i][0][17:20] != 'CYS'):
            continue
        # Found a cysteine, now track down the sulfur
        iX,iY,iZ=get_coords('SG',npdb[i])
                
        for j in range(i+1,len(npdb)):
            if (npdb[j][0][17:20]!= 'CYS'):
                continue
            (jX,jY,jZ) = get_coords('SG',npdb[j])
                
            dX=iX-jX
            dY=iY-jY
            dZ=iZ-jZ

            distance = math.sqrt(dX*dX+dY*dY+dZ*dZ)
            if (distance >= min_dist and distance <= max_dist):
                residues_to_rename = residues_to_rename | set([i,j])
        
    # Rename the residues we selected
    for i in residues_to_rename:
        npdb[i]=map(lambda x:x.replace('CYS','CYS2'),npdb[i])           
        
    return

def atom_is_present(pdblines, atomname):
    """Returns TRUE if the given atom is present in the given PDB atom lines.
    
    ARGUMENTS
        pdblines - list of PDB lines
        atomname - the name of the atom to check the existence of

    RETURNS
        is_present - True if the given atom name is present, False otherwise

    """

    is_present = False
    for pdbline in pdblines:
        if (pdbline[13:16] == atomname):
            is_present = True

    return is_present

def histidine_search(npdb):
    """Rename HIS residues to HID, HIE, or HIP by examining which protons are present.
    
    ARGUMENTS
        npdb - nested  PDB
    """
    
    for i in range(len(npdb)):
        resname = npdb[i][0][17:20]
        if (resname == 'HIS'):
            HE_present = atom_is_present(npdb[i], 'HE2') # bonded to NE2
            HD_present = atom_is_present(npdb[i], 'HD1') # bonded to ND1
            
            if (HD_present and HE_present):
                npdb[i]=map(lambda x:x.replace('HIS','HIP'),npdb[i])                
            elif (HE_present):
                npdb[i]=map(lambda x:x.replace('HIS','HIE'),npdb[i])
            elif (HD_present):
                npdb[i]=map(lambda x:x.replace('HIS','HID'),npdb[i])
            else:
                raise "No protons found for histidine."

    return

def pdb_cleanup(pdbarr):
    """Remove extraneous entries from MCCE PDB files.
    
    ARGUMENTS
        pdbarr - list of PDB lines
        
    RETURNS
        updated_pdbarr - updated list of PDB lines
    """
    # Nuke everything after the coordinates
    # Use a default occupancy of 1 and a default B-factor of 0
    for i in range(len(pdbarr)):
        atom=pdbarr[i][0:54]
        # The atom symbol will be the first character after stripping
        # whitespace and numbers
        # NOTE: this will fail for atoms with more than one letter in their symbol
        #            eg, counterions
        atomsymbol = atom[12:16].strip(" 0123456789")[0]
        atom = atom + "1.00".rjust(6) + "0.00".rjust(6) + " "*10 + atomsymbol.rjust(2)+ " "*2
        pdbarr[i]=atom

    npdb=nest_pdb(pdbarr)

    for i in range(len(npdb)):
        # Set the residue index numbers
        # NOTE: We set chain identifier to blank, so this won't work for multi-chain PDBs
        for j in range(len(npdb[i])):
            atom = npdb[i][j]
            atom = atom[:21]+ " " + str(i+1).rjust(4)+ " "*4 + atom[30:]
            #print len(atom)
            npdb[i][j] = atom
            
    updated_pdbarr = unnest_pdb(npdb)
        
    return updated_pdbarr

            
<commit_msg>Committed most current versions of rename.py and mcce.py - have separated protonation state step from renaming steps<commit_after>import tempfile, os, sys, commands, re, math

"""
Routines to rename atoms and residues from a 'labeled' step2_out.pdb file  ("0","+", or "-" tacked on the end) 

REVISION LOG
- 8-21-2007: VAV; spun these functions into their own python file 
"""
def rename_residues(pdbarr, renameTermini=True):
    """Convert residue names (and some heavy atom names) to AMBER (ffamber) format.
    
    ARGUMENTS
        pdbarr - array of PDB atoms
        
    RETURNS
        pdbarr - updated array of PDB atoms with AMBER appropriate residue names
    """

    
    # DEBUG check.
    print "At start of renaming, pdbarr has",len(pdbarr),"items"
     
    # Transform the PDB file into a nested format that is easier to manipulate.
    npdb = nest_pdb(pdbarr)
    
    # DEBUG check.
    print "Nest/unnest leads to",len(unnest_pdb(npdb)),"items"

    # Convert residue and atom names to be appropriate for AMBER.
    rename_charged(npdb)
    histidine_search(npdb)
    disulfide_search(npdb)
    if (renameTermini):
        rename_termini(npdb)

    # Un-nest PDB
    pdbarr = unnest_pdb(npdb)
    
    # DEBUG check.
    #print "At end of renaming, pdbarr has ", len(pdbarr), "items"

    # Return the PDB file.
    return pdbarr


def rename_charged(npdb):
    """Generate AMBER-specific residue names for charged residues from MCCE residue names.
        
    ARGUMENTS
        npdb - "nested" PDB data structure generated by nest_pdb. Modified to reflect AMBER names.
    """
        
    resname_and_state = map(lambda x:x[-1][17:20]+x[-1][-1],npdb)
    #print resname_and_state
    for i in range(len(npdb)):
        if (resname_and_state[i]=='HIS+'):
            npdb[i]=map(lambda x:x.replace('HIS','HIP'),npdb[i])
        elif (resname_and_state[i]=='LYS0'):
            npdb[i]=map(lambda x:x.replace('LYS','LYN'),npdb[i])
        elif (resname_and_state[i]=='LYS+'):
            npdb[i]=map(lambda x:x.replace('LYS','LYP'),npdb[i])
        elif (resname_and_state[i]=='CYS0'):
            npdb[i]=map(lambda x:x.replace('CYS','CYN'),npdb[i])
        elif (resname_and_state[i]=='CYS-'):
            npdb[i]=map(lambda x:x.replace('CYS','CYM'),npdb[i])
        elif (resname_and_state[i]=='ASP0'):
            npdb[i]=map(lambda x:x.replace('ASP','ASH'),npdb[i])
        #Aspartate requires no sub
        elif (resname_and_state[i]=='GLU0'):
            npdb[i]=map(lambda x:x.replace('GLU','GLH'),npdb[i])
        # Glutamate requires no sub

    return


def labeledPDB_to_AmberPDB(labeledPDBfile, outPDBfile, renameResidues=True):
    """Reads a 'labeled' PDB file written by MCCE (with "0","+", or "-" indicating the protonation state - see mcce.py)
    Example:
    ATOM      1  CA  LEU _0001_000  -8.847   7.195  16.727   1.700       0.000      BK____M000
    0
    ATOM      2  C   LEU _0001_000  -9.019   7.821  15.386   1.700       0.000      BK____M000
    0
    ....

    Writes an AMBER-residue-formatted PDB file, outPDBfile.
    """

    # Read in the pdb lines
    fin = open(labeledPDBfile,'r')
    lines = fin.readlines()
    fin.close()
    # Format so the even-numbered lines, "0", "+", "-", etc are supposed to be tacked onto the ends of each atom's line
    pdbarr = []
    while len(lines) > 0:
        pdbline = lines.pop(0)
        label = lines.pop(0)
        pdbarr.append( pdbline + label )
 
    if renameResidues:
        pdbarr=rename_residues(pdbarr)
    pdbarr=pdb_cleanup(pdbarr)
    
    # write AmberPDB to file
    fout = open(outPDBfile,'w')
    for line in pdbarr:
        fout.write(line)
    fout.close()


def renumber_atoms(pdbarr):
    """Renumbers atom entries (list items) of a PDB file so they are in sequence starting from 1.
    
    ARGUMENTS
        pdbarr - list of PDB atoms to be renumbered; mutated by call
    """
    
    # Renumber the atom entries in place.
    for i in range(len(pdbarr)):
        pdbarr[i]=pdbarr[i][0:6]+str(i+1).rjust(5)+pdbarr[i][11:]        

def rename_termini(npdb):
    """Renames the 'NTR' and 'CTR' residues generated by MCCE.
    Renames the O and OXT atoms at the C-terminus to OC1 and OC2 for AMBER.
    
    ARGUMENTS
        npdb - 
    
    """

    # Get three-letter residue names of N- and C-termini.
    ntrname = npdb[1][0][17:20]
    ctrname = npdb[-2][0][17:20]

    # Rename N-terminal residue 'XXX' to 'NXXX'.
    for j in [0,1]:
        for i in range(len(npdb[j])):
            npdb[j][i] = npdb[j][i][0:17]+'N'+ntrname+npdb[j][i][21:]
            
    # Rename C-terminal residue from 'XXX' to 'CXXX', with some exceptions (noted below).
    # Also rename terminal oxygen atoms.
    for j in [-1,-2]:
        for i in range(len(npdb[j])):
            # VAV 5/14/2007: gmx ffamber does not have CLYN, only CLYP 
            if ctrname=='LYN':      
                npdb[j][i] = npdb[j][i][0:17]+'CLYP'+npdb[j][i][21:]
            else:   
                npdb[j][i] = npdb[j][i][0:17]+'C'+ctrname+npdb[j][i][21:]
                            
            # DLM 5/9/2007: Though this says that it changes the O and OXT names, it doesn't
            # Adding the below to do the renaming.
            if npdb[j][i][13:16].split()[0]=='O':
                npdb[j][i] = npdb[j][i][0:13]+'OC1'+npdb[j][i][16:]
            elif npdb[j][i][13:16].split()[0]=='OXT':
                npdb[j][i] = npdb[j][i][0:13]+'OC2'+npdb[j][i][16:]
                
            #VAV 6/5/2007: gmx ffamber doesn't have "CD1" as an ILE/CILE atomname, only "CD"
            if ctrname=='ILE':      
                if npdb[j][i][13:16].split()[0]=='CD1':
                    npdb[j][i] = npdb[j][i][0:13]+'CD '+npdb[j][i][16:]
        
def nest_pdb(pdbarr):
    """Collect lines of the PDB file by residue.
    
    ARGUMENTS:
        pdbarr - list of lines from PDB file
        
    RETURNS
        nestedpdb - nested PDB file, such that nestedpdb[i][j] will be line j from residue i.
    """
    
    nestedpdb=[]
    residue=[]
    for line in pdbarr:
        if (len(residue) == 0):
            residue.append(line)
        else:
            if (line[17:27] == residue[-1][17:27]):
                residue.append(line)
            else:
                nestedpdb.append(residue)
                residue=[line]
    nestedpdb.append(residue)
    
    # Return the nexted PDB file.
    return nestedpdb

def unnest_pdb(npdb):
    """Expand the lines from the "nested" PDB file into a flat list of lines.
    
    ARGUMENTS
        npdb - "nested" PDB file; mutated in place
        
    RETURNS
        pdbarr - un-nested PDB file, such that pdbarr[i] will be line i from the PDB file
    """
    
    pdbarr=[]
    for res in npdb:
        for atm in res:
            pdbarr.append(atm)
    return pdbarr
    

def get_coords(atomname,residue):
    """Extract coordinates from given atom in list of PDB lines (such as all those lines for a residue in a nested PDB file).
    
    ARGUMENTS
        atomname - the name of the atom
        residue - list of lines from PDB file

    RETURNS
        (x,y,z) - tuple of atom coordinates, as floats
    
    """
    for n in range(len(residue)):
        if (residue[n][12:16].strip() == atomname.strip()):
            iX=float(residue[n][30:38])
            iY=float(residue[n][38:46])
            iZ=float(residue[n][46:54])
            return (iX,iY,iZ)

    raise "Atom not found!"        

def disulfide_search(npdb, min_dist = 1.8, max_dist = 2.2):
    """Rename CYS to CYX if it participates in a disulfide bond, as judged by distance range (inclusive).
    
    ARGUMENTS
        npdb - nested PDB file
        
    OPTIONAL ARGUMENTS
        min_dist - minium distance cutoff for perceiving disulfide bond (default 1.8 A)
        max_dist - maximum distance cutoff for perceiving disulfide bond (default 2.2 A)
    """
        
    residues_to_rename=set([])
    for i in range(len(npdb)):
        if (npdb[i][0][17:20] != 'CYS'):
            continue
        # Found a cysteine, now track down the sulfur
        iX,iY,iZ=get_coords('SG',npdb[i])
                
        for j in range(i+1,len(npdb)):
            if (npdb[j][0][17:20]!= 'CYS'):
                continue
            (jX,jY,jZ) = get_coords('SG',npdb[j])
                
            dX=iX-jX
            dY=iY-jY
            dZ=iZ-jZ

            distance = math.sqrt(dX*dX+dY*dY+dZ*dZ)
            if (distance >= min_dist and distance <= max_dist):
                residues_to_rename = residues_to_rename | set([i,j])
        
    # Rename the residues we selected
    for i in residues_to_rename:
        npdb[i]=map(lambda x:x.replace('CYS','CYS2'),npdb[i])           
        
    return

def atom_is_present(pdblines, atomname):
    """Returns TRUE if the given atom is present in the given PDB atom lines.
    
    ARGUMENTS
        pdblines - list of PDB lines
        atomname - the name of the atom to check the existence of

    RETURNS
        is_present - True if the given atom name is present, False otherwise

    """

    is_present = False
    for pdbline in pdblines:
        if (pdbline[13:16] == atomname):
            is_present = True

    return is_present

def histidine_search(npdb):
    """Rename HIS residues to HID, HIE, or HIP by examining which protons are present.
    
    ARGUMENTS
        npdb - nested  PDB
    """
    
    for i in range(len(npdb)):
        resname = npdb[i][0][17:20]
        if (resname == 'HIS'):
            HE_present = atom_is_present(npdb[i], 'HE2') # bonded to NE2
            HD_present = atom_is_present(npdb[i], 'HD1') # bonded to ND1
            
            if (HD_present and HE_present):
                npdb[i]=map(lambda x:x.replace('HIS','HIP'),npdb[i])                
            elif (HE_present):
                npdb[i]=map(lambda x:x.replace('HIS','HIE'),npdb[i])
            elif (HD_present):
                npdb[i]=map(lambda x:x.replace('HIS','HID'),npdb[i])
            else:
                raise "No protons found for histidine."

    return

def pdb_cleanup(pdbarr):
    """Remove extraneous entries from MCCE PDB files.
    
    ARGUMENTS
        pdbarr - list of PDB lines
        
    RETURNS
        updated_pdbarr - updated list of PDB lines
    """
    # Nuke everything after the coordinates
    # Use a default occupancy of 1 and a default B-factor of 0
    for i in range(len(pdbarr)):
        atom=pdbarr[i][0:54]
        # The atom symbol will be the first character after stripping
        # whitespace and numbers
        # NOTE: this will fail for atoms with more than one letter in their symbol
        #            eg, counterions
        atomsymbol = atom[12:16].strip(" 0123456789")[0]
        atom = atom + "1.00".rjust(6) + "0.00".rjust(6) + " "*10 + atomsymbol.rjust(2)+ " "*2
        pdbarr[i]=atom

    npdb=nest_pdb(pdbarr)

    for i in range(len(npdb)):
        # Set the residue index numbers
        # NOTE: We set chain identifier to blank, so this won't work for multi-chain PDBs
        for j in range(len(npdb[i])):
            atom = npdb[i][j]
            atom = atom[:21]+ " " + str(i+1).rjust(4)+ " "*4 + atom[30:]
            #print len(atom)
            npdb[i][j] = atom
            
    updated_pdbarr = unnest_pdb(npdb)
        
    return updated_pdbarr

            

{'fraction_non_alphanumeric': 0.12398364853330937, 'fraction_numerical': 0.025245945824536185, 'mean_word_length': 2.137260428410372, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7206173', 'n_tokens_mistral': 8226, 'n_tokens_neox': 7619, 'n_words': 2066}