<jupyter_start><jupyter_text>"fpocket 4.0 - Towards biologics"> "A quick demo of all the new fpocket features with special focus on protein-peptide interfaces"- toc: true- branch: master- badges: true- comments: true- author: <NAME>- categories: [fpocket, pockets, cavity, PPI, protein, epitope, DNA, RNA, antibody]- image: images/fpocket-logo.png I just released the last version of fpocket (4.0) which contains quite a few new features. Most of the work was done by <NAME>, an eager summer intern we took in at Discngine during last year's COVID-mess summer. He did a great job under these circumstances (remote work in a new company etc). So kudos to him for these new additions.  MMCIF supportFirst of all, fpocket now supports mmCIF. Yay ... yet another very painful task to integrate for yet another painfully crappy molecular file format (I admire the openbabel folks ...  ;) ).  Same as for the PDB format we forked essentially an integration from the VMD molfile plugin and extended that a bit further. Mael also integrated mmCIF writers so all output files can be written in this new crappy format as well. Essentially you now have the choice to write either only pdb or cif or both. So everybody should be happy (nah ... that's never going to be case anyway).  New command line flagsThe most interesting part (for me) of this release are additional command line flags and behind these several use cases that can be addressed now. These command line arguments enable mainly to work in a more controlled way with multi chain protein structures. This can be particularly useful if:1. you want to study only one domain of one of these novel gigantic cryoEM structures2. you want to study in detail a protein - biomolecule binding site (so other than small molecule)3. you want to learn key characteristics of peptide/protein binding sites vs RNA/DNA binding sites vs drug binding sites```    -c char : (Default is none): Use this flag to choose which chains you want to delete before running fpocket. The selected chains can be specified with ',' or ':' delimiters, for example you can use it '-c B,D' or '-c B:D'. You can delete up to 20 different chains.``` This flag essentially allows you to explicitly delete particular chain(s) before doing the fpocket run. This allows you to identify (without any bias) clefts that might occur on the PPI that were inaccessible before for pocket prediction. ```    -k char : (Default is none): Use this flag to choose which chains you want to keep before running fpocket. The selected chains can be specified with ',' or ':' delimiters, for example you can use it '-k A,B,E' or '-k A:B:E'. You can keep up to 20 different chains.```This is essentially the inverse operator of the previous flag. If you do not want to list n chains, you can decide to keep only the relevant ones here. ```    -a char : (Default is none): With this flag you can select a chain you want to be considered as a ligand. Works the same way as the "-r" flag but with a whole chain. Only a single chain can be chosen, for example '-a D'.```This flag essentially allows you to target for only one particular binding epitope. If you want to extract descriptors for only that binding site where the select chain is located then you can use the `-a` flag. Fpocket will hide the chain specified via the flag. Run vertice detection and then will choose all vertices overlapping with the chain you specified in `-a`. These will then be clustered together into a final "pocket".  This is rather generic and can be applied on protein chains, RNA, DNA ... whatever you like given that it is defined in a seperate chain. It is so generic, that in theory you can do even very nasty stuff, but let's see if you can find that out yourselves.This argument is particularly interesting when you want to extract examples (3D patches) or numerical descriptors (polar/apolar ASA, curvature, buriedness ...) for other binding sites than classical small molecule binding sites. It essentially lies out the basics for training & learning what such other binding epitopes look like. You can use that for a fancy model, or to characterize binding epitopes ultimately.  MiscellaneousSeperately, several things were added to the [github repo of fpocket](https://github.com/Discngine/fpocket), namely the old documentation was ported to markdown (that was painful).There's also a tiny testing environment, CI/CD pipes finally set up to check compilation and unit tests upon PR's. Last, an official fpocket docker image is also available on [dockerhub](https://hub.docker.com/r/fpocket/fpocket/).   Detecting binding epitopes Enough text now. Let's check what you can do with that now. Here we have a nice example of a [TEAD4 / YAP complex](https://www.rcsb.org/structure/6HIK). I'm interested in studying the clefts covered by YAP on the TEAD4 structure. Before I essentially had to drop YAP from the PDB run fpocket and put back in YAP to see where it is located. Here YAP corresponds to chain L (the red one) and TEAD4 (the grey one).<jupyter_code>#collapse-hide

import py3Dmol
viewer = py3Dmol.view(query='pdb:6hik',width=300, height=300)
viewer.setStyle({'chain':'A'},{'cartoon':{'color':'grey'}})
viewer.setStyle({'chain':'L'},{'cartoon':{'color':'red'}})
viewer.zoomTo()<jupyter_output><empty_output><jupyter_text>We can now for instance launch fpocket on TEAD4 by either dropping YAP or by keeping TEAD4 explicitly. Whatever mindset you prefer ;) Let's be in a positive mood and keep TEAD4: ```fpocket -f 6hik.pdb -k A```<jupyter_code>#collapse-hide

import py3Dmol

viewer = py3Dmol.view(query='pdb:6hik',width=500, height=500)
pdbdata=open('../data/6hik_out.pdb', 'r').read()
viewer.addModel(pdbdata,'pdb')
viewer.setStyle({'chain':'A'},{'cartoon':{'color':'grey'}})
viewer.setStyle({'chain':'L'},{'cartoon':{'color':'red'},'line':{}})
viewer.setStyle({'chain':'C'},{'sphere':{'colorscheme':{'prop':'resi','gradient':'roygb','min':1,'max':15}}})
#viewer.setStyle({'chain':'C'},{'line':{'colorfunc':colorAsSnake}})

viewer.zoomTo()<jupyter_output><empty_output><jupyter_text>When doing this in an automatic way, so with classical fpocket parameters, you'll now be able to identify pockets that overlap with YAP. That's kind of nice, but what if I want to study really in detail all the binding interface itself between YAP and TEAD4?In order to instruct fpocket to define a pocket only on the interface itself you can now do something like that: ```bashfpocket -f 6hik.pdb -a L```<jupyter_code>#collapse-hide

import py3Dmol

viewer = py3Dmol.view(query='pdb:6hik',width=500, height=500)
pdbdata=open('../data/6hik_out_explicit.pdb', 'r').read()
viewer.addModel(pdbdata,'pdb')
viewer.setStyle({'chain':'A'},{'cartoon':{'color':'grey'}})
viewer.setStyle({'chain':'L'},{'cartoon':{'color':'red'},'line':{}})
viewer.setStyle({'chain':'C'},{'sphere':{'color':'orange'}})<jupyter_output><empty_output><jupyter_text>You can now get the full orange blob as the interacting epitope. Together with that you get the usual statistics etc. The cool thing is, this can be very well used on protein RNA interfaces as well.  A word on druggabilityIf you want to assess druggability, please use the default fpocket parameters and not a guided pocket detection or other pocket detection parameters. The druggability assessment intends to estimate a pockets tractability for small molecule binding sites. Applying this to such larger surfaces is out of the applicability domain here in my opinion.   Extracting descriptors with fpocketOk running these things on individual examples is nice. You can automate things with fpocket quite easily and extract descriptors on a larger scale using the fancy `-d` flag (oh yes, yet another flag) if the only thing you're after is descriptor extraction.  Why would I ever want to extract descriptors?Calculating descriptors on binding epitopes is generating the basis for using these to derive several potential applications afterwards. For instance, to derive the druggability prediction in fpocket, first we extracted descriptors of known druggable and supposed non-druggable (no religious debate here) pockets, second we determined the most relevant descriptors and last we trained a scoring function based on these descriptors. You can very well imagine to do the same on other types of binding epitopes (there are a few papers out there on that already I guess), like protein protein interfaces, antibody / antigene interfaces more specifically, crystal contacts, protein DNA/RNA interfaces etc...Ultimately these characterisations will allow you to train functions or super fancy __deep__ learning models (if you really need that).  In order to do that on a larger scale you could use dpocket (a less well-known sibling of fpocket), but it currently still doesn't fully integrate the logic with selecting chains as ligands (still a bit of work needed here). So let's stick to fpocket on a larger scale example for descriptor extraction which will work just fine: ```data/peptide_data/3uqp.pdb Bdata/peptide_data/3uri.pdb Bdata/peptide_data/4rcp.pdb Bdata/peptide_data/4tpg.pdb Edata/peptide_data/5jxh.pdb H```Let's suppose we have a csv file like this one above. It defines a set of PDB files and the chain we want to consider as a ligand explicitly during the fpocket run to assess the binding epitope this chain is making with the "receptor". You could run this in a very geeky way using this here: `awk '{ print "fpocket -f " $1 " -a "$2 " -d"}' list.txt` | shThis will output a messy thing like this: ```csvcav_id drug_score volume nb_asph inter_chain apol_asph_proportion mean_asph_radius as_density mean_asph_solv_acc mean_loc_hyd_dens flex hydrophobicity_score volume_score charge_score polarity_score a0_apol a0_pol af_apol af_pol n_abpa ala cys asp glu phe gly his ile lys leu met asn pro gln arg ser thr val trp tyr chain_1_type chain_2_type num_res_chain_1 num_res_chain_2 lig_het_tag name_chain_1 name_chain_21 0.0021 4091.2471 213 0 0.1831 4.3624 10.8343 0.5065 10.3590 0.0000 17.6957 4.1304 3 31 176.3137 207.7781 91.7446 196.2307 21 1 4 2 3 0 3 1 4 0 4 1 3 0 1 1 6 4 2 3 3 0 0 457 457 PSA A Acav_id drug_score volume nb_asph inter_chain apol_asph_proportion mean_asph_radius as_density mean_asph_solv_acc mean_loc_hyd_dens flex hydrophobicity_score volume_score charge_score polarity_score a0_apol a0_pol af_apol af_pol n_abpa ala cys asp glu phe gly his ile lys leu met asn pro gln arg ser thr val trp tyr chain_1_type chain_2_type num_res_chain_1 num_res_chain_2 lig_het_tag name_chain_1 name_chain_21 0.1847 4801.9624 204 0 0.3529 4.4814 10.8531 0.5450 22.0000 0.0000 27.5882 3.4706 -8 26 183.5595 226.1442 124.2374 226.0680 11 1 0 0 7 0 2 1 7 0 8 4 0 0 4 0 9 5 0 2 1 0 0 328 328 NULL A Acav_id drug_score volume nb_asph inter_chain apol_asph_proportion mean_asph_radius as_density mean_asph_solv_acc mean_loc_hyd_dens flex hydrophobicity_score volume_score charge_score polarity_score a0_apol a0_pol af_apol af_pol n_abpa ala cys asp glu phe gly his ile lys leu met asn pro gln arg ser thr val trp tyr chain_1_type chain_2_type num_res_chain_1 num_res_chain_2 lig_het_tag name_chain_1 name_chain_21 1.0000 2998.7053 187 0 0.5241 4.3207 11.5341 0.5341 48.6735 0.0000 34.3462 5.0769 3 16 107.4789 84.2598 47.7836 96.3200 8 1 2 1 2 0 0 1 1 1 1 3 3 0 3 0 1 0 1 4 1 0 0 238 238 NULL A Acav_id drug_score volume nb_asph inter_chain apol_asph_proportion mean_asph_radius as_density mean_asph_solv_acc mean_loc_hyd_dens flex hydrophobicity_score volume_score charge_score polarity_score a0_apol a0_pol af_apol af_pol n_abpa ala cys asp glu phe gly his ile lys leu met asn pro gln arg ser thr val trp tyr chain_1_type chain_2_type num_res_chain_1 num_res_chain_2 lig_het_tag name_chain_1 name_chain_21 0.9998 1552.9973 110 0 0.4182 4.5201 6.7627 0.4934 36.8261 0.0000 34.2593 4.8519 0 17 105.0637 48.2599 82.1878 27.8238 12 1 1 2 0 0 3 2 0 0 2 1 1 0 3 1 2 0 0 6 2 0 0 594 594 NULL A Acav_id drug_score volume nb_asph inter_chain apol_asph_proportion mean_asph_radius as_density mean_asph_solv_acc mean_loc_hyd_dens flex hydrophobicity_score volume_score charge_score polarity_score a0_apol a0_pol af_apol af_pol n_abpa ala cys asp glu phe gly his ile lys leu met asn pro gln arg ser thr val trp tyr chain_1_type chain_2_type num_res_chain_1 num_res_chain_2 lig_het_tag name_chain_1 name_chain_21 0.0044 3875.7029 152 0 0.1513 4.4958 9.8519 0.5972 9.5652 0.0000 -7.5946 3.4865 -7 27 114.7247 187.8280 70.7198 141.5900 15 2 3 2 8 0 0 4 4 2 0 1 0 0 0 2 2 4 1 1 1 0 0 499 499 NULL A A``` Let's clean this up the geeky way: `awk '{ print "fpocket -f " $1 " -a "$2 " -d"}' list.txt | sh | awk '{if(NR%2==0)print}'`This should only give the descriptors: ```1 0.0021 3972.3301 213 0 0.1831 4.3624 10.8343 0.5065 10.3590 0.0000 17.6957 4.1304 3 31 176.3137 207.7781 91.7446 196.2307 21 1 4 2 3 0 3 1 4 0 4 1 3 0 1 1 6 4 2 3 3 0 0 457 457 PSA A A1 0.1847 4604.5068 204 0 0.3529 4.4814 10.8531 0.5450 22.0000 0.0000 27.5882 3.4706 -8 26 183.5595 226.1442 124.2374 226.0680 11 1 0 0 7 0 2 1 7 0 8 4 0 0 4 0 9 5 0 2 1 0 0 328 328 NULL A A1 1.0000 2895.7646 187 0 0.5241 4.3207 11.5341 0.5341 48.6735 0.0000 34.3462 5.0769 3 16 107.4789 84.2598 47.7836 96.3200 8 1 2 1 2 0 0 1 1 1 1 3 3 0 3 0 1 0 1 4 1 0 0 238 238 NULL A A1 0.9998 1574.4559 110 0 0.4182 4.5201 6.7627 0.4934 36.8261 0.0000 34.2593 4.8519 0 17 105.0637 48.2599 82.1878 27.8238 12 1 1 2 0 0 3 2 0 0 2 1 1 0 3 1 2 0 0 6 2 0 0 594 594 NULL A A1 0.0044 3794.7029 152 0 0.1513 4.4958 9.8519 0.5972 9.5652 0.0000 -7.5946 3.4865 -7 27 114.7247 187.8280 70.7198 141.5900 15 2 3 2 8 0 0 4 4 2 0 1 0 0 0 2 2 4 1 1 1 0 0 499 499 NULL A A``` I did this on a larger list of peptide binding structures and this is typically a good way to start studying some of the properties of these epitopes compared to drug binding sites or others.<jupyter_code>import pandas as pd
import altair as alt
d=pd.read_table('../data/peptide_out.txt',delim_whitespace=True,header=None,names=["cav_id","drug_score","volume","nb_asph","inter_chain","apol_asph_proportion","mean_asph_radius","as_density","mean_asph_solv_acc","mean_loc_hyd_dens","flex","hydrophobicity_score","volume_score","charge_score","polarity_score","a0_apol","a0_pol","af_apol","af_pol","n_abpa","ala","cys","asp","glu","phe","gly","his","ile","lys","leu","met","asn","pro","gln","arg","ser","thr","val","trp","tyr","chain_1_type","chain_2_type","num_res_chain_1","num_res_chain_2","lig_het_tag","name_chain_1","name_chain_2"])
alt.Chart(d).mark_bar().encode(
    alt.X("volume", bin=True),
    y='count()',
)

alt.Chart(d).transform_fold( 
    ['a0_pol', 'a0_apol'], 
    as_=['Columns', 'Values'] 
).mark_area( 
    opacity=0.5, 
    interpolate='step'
).encode( 
    alt.X('Values:Q', bin=alt.Bin(maxbins=20)), 
    alt.Y('count()', stack=None), 
    alt.Color('Columns:N') 
)<jupyter_output><empty_output><jupyter_text>Here we can see the distribution of polar versus apolar surface areas of these peptide binding sites under consideration. What is interesting already is that we can clearly distinguish this epitope with these descriptors alone from typical drug binding sites (more hydrophobic in general).<jupyter_code>alt.Chart(d).mark_circle(size=60).encode(
    x='n_abpa',
    y='a0_apol',
    color='volume',
    tooltip=['drug_score', 'mean_asph_solv_acc', 'a0_apol','a0_pol',"ala","cys","phe","gly","ile","leu","met","pro","val","trp","tyr"]
).interactive()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08326853291861068, 'fraction_numerical': 0.08540694660445827, 'mean_word_length': 4.767189835575485, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 24, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '505135', 'license': 'Apache-2.0', 'n_tokens_mistral': 5922, 'n_tokens_neox': 4750, 'n_words': 2180, 'path': '_notebooks/2021-02-02-fpocket-towards-biologics.ipynb', 'repo_name': 'pschmidtke/blog', 'size': '15447'}