In-depth analysis of the sidechain separation distance in a set of PDB files

analyze_peaks.py

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
from marshal import dump, load
from random import sample
from glob import glob
from pdb import *
import sys



def run_experiment(name, residues, interests):
  left, right = [], []
  matching_names = map(lambda i: i[0], interests)
  for residue in residues:
    for interest in interests:
      res_name, a1n, a2n, a3n, a4n, divide = interest
      if name != res_name: continue
      if res_name == residue['resName']:
        distance = residue_distance(residue, a3n, a4n)
        if distance == None:
          continue
        elif distance < divide:
          left.append(distance)
        elif distance >= divide:
          right.append(distance)
  return [name, left, right]

def analyze_result(results, mind=2.0, maxd=3.0, nbins=50):
  name, left, right  = results
  data, labels = [], []
  left         = filter(lambda d: mind < d < maxd, left)
  right        = filter(lambda d: mind < d < maxd, right)
  print(len(left), ' | ', len(right))
  nleft = len(left)
  if nleft == 0:
    med_left = 0.0
  else:
    med_left = median(left)
  med_right    = median(right)
  delta        = med_right - med_left
  print("Left Peak: %.2f Angstrom, Right Peak: %.2f Angstrom, Delta: %.2f Angstrom" %
    (med_left, med_right, delta))
  if nleft == 0:
    figure(figsize=[11, 8.5])
    hist(right, histtype='bar', bins=nbins, label="right")
  else:
    data.append(left)
    data.append(right)
    labels.append("left")
    labels.append("right")
    figure(figsize=[11, 8.5])
    hist(data, histtype='bar', bins=nbins, label=labels)
  legend()
  grid()
  savefig(name + '.png', dpi=300)


interests = [
  ['ASP', 'CB', 'CG', 'OD1', 'OD2', 2.16],
  ['GLU', 'CG', 'CD', 'OE1', 'OE2', 2.16],
  ['ASN', 'CB', 'CG', 'OD1', 'ND2', 2.20],
  ['GLN', 'CG', 'CD', 'OE1', 'NE2', 2.20],
  ['VAL', 'CA', 'CB', 'CG1', 'CG2', 0.00]
]

try:
  __IPYTHON__
  is_ipy = True
except: is_ipy = False

if is_ipy or len(sys.argv) == 1: num_samples = 256
else: num_samples = int(sys.argv[1])
residue_data = 'r%d.pydata' % num_samples

if is_ipy:
  residues_file = open(residue_data, 'rb')
  all_residues  = load(residues_file)
  residues_file.close()

asp = run_experiment('ASP', all_residues, interests)
asn = run_experiment('ASN', all_residues, interests)
glu = run_experiment('GLU', all_residues, interests)
gln = run_experiment('GLN', all_residues, interests)
val = run_experiment('VAL', all_residues, interests)