Examine the relationship between the sidechain separation and local pKa

analyze_separations_and_pka.py

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
from __future__ import division
from __future__ import print_function
from marshal import dump, load
from random import sample
from glob import glob
from pdb import *
import subprocess as sub
import sys
import re

def run_experiment(name, residues, interests):
  return None

def analyze_result(left, right, mind=2.1, maxd=2.3, nbins=50):
  return None




try:
  __IPYTHON__
  is_ipy = True
except: is_ipy = False

mind, maxd  = 2.1, 2.3
asp_results = []
glu_results = []
num_samples = 64
propka_cmd  = "propka %s | grep '^  *[AG][SL][PU]'"
nrpdb       = sample(glob('nrpdb_ph/*'), k=num_samples)

nprocd = 0
for pdb_fname in nrpdb:
  try:
    pka_data = sub.check_output(propka_cmd % pdb_fname, shell=True)
  except: continue

  interests = {}
  for line in pka_data.split("\n"):
    line = re.compile('\s+').split(line.strip())
    if len(line) != 5: continue
    name, seq, chain_id, pka, model_pka = line
    name      = name.strip()
    seq       = int(seq)
    chain_id  = chain_id.strip()
    pka       = float(pka)
    model_pka = float(model_pka)
    interests[seq] = [name, seq, pka, model_pka]
  interesting_seqs = interests.keys()

  last_seq = None
  for residue in extract_models(pdb_fname)[0]['residues']:
    if residue['resSeq'] in interesting_seqs:
      name, seq, pka, model_pka = interests[residue['resSeq']]
      a3n, a4n = None, None
      if name == residue['resName'] and name == 'ASP':
        a3n, a4n = 'OD1', 'OD2'
      if name == residue['resName'] and name == 'GLU':
        a3n, a4n = 'OE1', 'OE2'
      distance = residue_distance(residue, a3n, a4n)
      if distance == None: continue
      if mind <= distance <= maxd:
        if name == 'ASP': asp_results.append([pka, distance])
        if name == 'GLU': glu_results.append([pka, distance])
    last_seq = residue['resSeq']

  nprocd += 1
  if nprocd % 25 == 0: print("%d processed" % nprocd)

sub.call('rm -rf *.propka_input *.pka', shell=True)
asp_xs = map(lambda a: a[0], asp_results)
asp_ys = map(lambda a: a[1], asp_results)
glu_xs = map(lambda a: a[0], glu_results)
glu_ys = map(lambda a: a[1], glu_results)