<jupyter_start><jupyter_text>generatePocketPDBsWe want to attempt to run our SOAP pipeline just on the pockets, to eliminate the influence of the surface. For each file in the given xyz file, run fpocket, get the residues touching the pocket, then write them all backout to a new xyz file<jupyter_code>import quippy
import subprocess
import os
import glob
import ase
input_xyz_file = "intermediates/samples_filtered.xyz"
input_file = quippy.AtomsReader(input_xyz_file, format="xyz")
labels = []
for protein in input_file:
    label = protein.params["tag"]
    labels.append(label[:5])

os.mkdir("temp")

for label in labels:
    pdb_ref = label[:4]
    chain_ref = label[4]
    subprocess.call(["wget","http://www.rcsb.org/pdb/files/{}.pdb".format(pdb_ref)])
    try:
        with open("{}.pdb".format(pdb_ref)) as flines:
            pdb_data = flines.readlines()
    except IOError:
        print("IO Error")
    single_chain_data = []
    for line in pdb_data:
        if line.startswith("ATOM"):
            if line[21] == chain_ref:
                single_chain_data.append(line.strip())

    with open("temp/{}.pdb".format(label), "w") as flines:
        flines.write("\n".join(single_chain_data))
    
# For each pdb file in the temporary directory, run fpocket (stored in ~/.local/bin here)
for path in os.listdir("temp"):
    print(path)
    # For whatever reason, standard IPython syntax isn't working
    command="cd temp && fpocket -f {} && cat {}_out/pockets/*_atm.pdb > {}_pockets.pdb".format(
    path, path[:5], path[:5])
    subprocess.call(command, shell=True)
# Get the indices for the residues touching the pockets
def getIndicesOfPocketResidues(pathname):
    with open(pathname) as flines:
        residueIds = set([int(line[22:26].strip()) for line in flines if line[:4] == "ATOM"])
    return residueIds
# Given that a _pockets file has a source pdb file in the same directory, open that, and extract all atoms in the pocket.
def getPocketAtoms(pathname):
    residueIds = getIndicesOfPocketResidues(pathname)
    sourcePdb = pathname[:-12] + ".pdb" # BRITTLE
    print(sourcePdb)
    with open(sourcePdb) as flines:
        data = [line.strip() for line in flines if int(line[22:26].strip()) in residueIds]
    outPdb = sourcePdb[:-4] + "_fullpocket.pdb"
    print(outPdb)
    with open(outPdb, "w") as outflines:
        outflines.write("\n".join(data))
for path in glob.glob("temp/*_pockets.pdb"):
    getPocketAtoms(path)
!cat temp/2biyA_fullpocket.pdb<jupyter_output>ATOM     25  N   LYS A  76      33.606  57.689  35.415  1.00 42.38           N
ATOM     26  CA  LYS A  76      32.547  58.689  35.429  1.00 40.48           C
ATOM     27  C   LYS A  76      32.844  59.716  36.511  1.00 37.69           C
ATOM     28  O   LYS A  76      33.215  59.360  37.630  1.00 36.99           O
ATOM     29  CB  LYS A  76      31.180  58.044  35.649  1.00 40.61           C
ATOM     30  CG  LYS A  76      30.054  58.727  34.877  1.00 43.50           C
ATOM     31  CD  LYS A  76      30.020  58.357  33.382  1.00 43.96           C
ATOM     32  CE  LYS A  76      28.711  58.804  32.728  1.00 44.97           C
ATOM     33  NZ  LYS A  76      28.840  59.255  31.312  1.00 44.29           N
ATOM     34  N   LYS A  77      32.735  60.990  36.161  1.00 34.18           N
ATOM     35  CA  LYS A  77      32.946  62.052  37.128  1.00 32.65           C
ATOM     36  C   LYS A  77      31.739  62.152  38.077  1.00 31.20           C
ATOM     37  O   LYS A  77      30.689  [...]<jupyter_text>So we have all the pockets. Tag them all again and stick into a .xyz file<jupyter_code>for 

input_file = quippy.AtomsReader(input_xyz_file, format="xyz")
full_labels = []
for protein in input_file:
    label = protein.params["tag"]
    full_labels.append(label)

protein_list = []
labels = []
for path in glob.glob("temp/*_fullpocket.pdb"):
    temp = quippy.Atoms(ase.io.read(path, format='proteindatabank')) # done like this to catch warning
    protein_list.append(temp)
    reference = path[5:10]
    tag = None
    for full_label in full_labels:
        if full_label.startswith(reference):
            tag=full_label
    
    if tag is not None:
        labels.append(tag)
    else:
        raise IOError
testFamily = quippy.AtomsList(protein_list)
print(len(testFamily), "proteins converted to xyz")
print(len(labels), "labels listed")
testFamily.write("intermediates/samples_filtered_pockets.xyz")
xyzFileName= "intermediates/samples_filtered_pockets.xyz"
with open(xyzFileName) as flines:
    data = [line.strip() for line in flines]
newdata = []
counter = 0 
for line in data:
    if line.startswith("cutoff"):
        label = labels[counter]
        line += " tag=\"{}\" ".format(label)
        counter += 1
    newdata.append(line)
assert len(newdata) == len(data)
assert len(labels) == counter
with open("intermediates/samples_filtered_pockets_tagged.xyz", "w") as flines:
    flines.write("\n".join(newdata))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11090982286634461, 'fraction_numerical': 0.06964573268921095, 'mean_word_length': 2.661753868828298, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, '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}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '201430', 'license': 'MIT', 'n_tokens_mistral': 1999, 'n_tokens_neox': 1759, 'n_words': 499, 'path': 'deepdrugdataset/generatePocketPDBs.ipynb', 'repo_name': 'ProteinNetworks/protein-soap', 'size': '4968'}