|
| 1 | +import math |
| 2 | +import os |
| 3 | +import sys |
| 4 | + |
| 5 | +import matplotlib.pyplot as plt |
| 6 | +import numpy as np |
| 7 | +from Bio import PDB |
| 8 | +from matplotlib import colors |
| 9 | + |
| 10 | + |
| 11 | +def ramachandran(file_name_list): |
| 12 | + """ |
| 13 | + Main calculation and plotting definition |
| 14 | + :param file_name_list: List of PDB files to plot |
| 15 | + :return: Nothing |
| 16 | + """ |
| 17 | + # General variable for the background preferences |
| 18 | + rama_preferences = { |
| 19 | + "General": { |
| 20 | + "file": "data/pref_general.data", |
| 21 | + "cmap": colors.ListedColormap(['#FFFFFF', '#B3E8FF', '#7FD9FF']), |
| 22 | + "bounds": [0, 0.0005, 0.02, 1], |
| 23 | + }, |
| 24 | + "GLY": { |
| 25 | + "file": "data/pref_glycine.data", |
| 26 | + "cmap": colors.ListedColormap(['#FFFFFF', '#FFE8C5', '#FFCC7F']), |
| 27 | + "bounds": [0, 0.002, 0.02, 1], |
| 28 | + }, |
| 29 | + "PRO": { |
| 30 | + "file": "data/pref_proline.data", |
| 31 | + "cmap": colors.ListedColormap(['#FFFFFF', '#D0FFC5', '#7FFF8C']), |
| 32 | + "bounds": [0, 0.002, 0.02, 1], |
| 33 | + }, |
| 34 | + "PRE-PRO": { |
| 35 | + "file": "data/pref_preproline.data", |
| 36 | + "cmap": colors.ListedColormap(['#FFFFFF', '#B3E8FF', '#7FD9FF']), |
| 37 | + "bounds": [0, 0.002, 0.02, 1], |
| 38 | + } |
| 39 | + } |
| 40 | + |
| 41 | + # Read in the expected torsion angles |
| 42 | + __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__))) |
| 43 | + rama_pref_values = {} |
| 44 | + for key, val in rama_preferences.items(): |
| 45 | + rama_pref_values[key] = np.full((360, 360), 0, dtype=np.float64) |
| 46 | + with open(os.path.join(__location__, val["file"])) as fn: |
| 47 | + for line in fn: |
| 48 | + if not line.startswith("#"): |
| 49 | + # Preference file has values for every second position only |
| 50 | + rama_pref_values[key][int(float(line.split()[1])) + 180][int(float(line.split()[0])) + 180] = float( |
| 51 | + line.split()[2]) |
| 52 | + rama_pref_values[key][int(float(line.split()[1])) + 179][int(float(line.split()[0])) + 179] = float( |
| 53 | + line.split()[2]) |
| 54 | + rama_pref_values[key][int(float(line.split()[1])) + 179][int(float(line.split()[0])) + 180] = float( |
| 55 | + line.split()[2]) |
| 56 | + rama_pref_values[key][int(float(line.split()[1])) + 180][int(float(line.split()[0])) + 179] = float( |
| 57 | + line.split()[2]) |
| 58 | + |
| 59 | + normals = {} |
| 60 | + outliers = {} |
| 61 | + for key, val in rama_preferences.items(): |
| 62 | + normals[key] = {"x": [], "y": []} |
| 63 | + outliers[key] = {"x": [], "y": []} |
| 64 | + |
| 65 | + # Calculate the torsion angle of the inputs |
| 66 | + for inp in file_name_list: |
| 67 | + if not os.path.isfile(inp): |
| 68 | + print("{} not found!".format(inp)) |
| 69 | + continue |
| 70 | + structure = PDB.PDBParser().get_structure('input_structure', inp) |
| 71 | + for model in structure: |
| 72 | + for chain in model: |
| 73 | + polypeptides = PDB.PPBuilder().build_peptides(chain) |
| 74 | + for poly_index, poly in enumerate(polypeptides): |
| 75 | + phi_psi = poly.get_phi_psi_list() |
| 76 | + for res_index, residue in enumerate(poly): |
| 77 | + res_name = "{}".format(residue.resname) |
| 78 | + res_num = residue.id[1] |
| 79 | + phi, psi = phi_psi[res_index] |
| 80 | + if phi and psi: |
| 81 | + if str(poly[res_index + 1].resname) == "PRO": |
| 82 | + aa_type = "PRE-PRO" |
| 83 | + elif res_name == "PRO": |
| 84 | + aa_type = "PRO" |
| 85 | + elif res_name == "GLY": |
| 86 | + aa_type = "GLY" |
| 87 | + else: |
| 88 | + aa_type = "General" |
| 89 | + if rama_pref_values[aa_type][int(math.degrees(psi)) + 180][int(math.degrees(phi)) + 180] < \ |
| 90 | + rama_preferences[aa_type]["bounds"][1]: |
| 91 | + print("{} {} {} {}{} is an outlier".format(inp, model, chain, res_name, res_num)) |
| 92 | + outliers[aa_type]["x"].append(math.degrees(phi)) |
| 93 | + outliers[aa_type]["y"].append(math.degrees(psi)) |
| 94 | + else: |
| 95 | + normals[aa_type]["x"].append(math.degrees(phi)) |
| 96 | + normals[aa_type]["y"].append(math.degrees(psi)) |
| 97 | + |
| 98 | + # Generate the plots |
| 99 | + for idx, (key, val) in enumerate(sorted(rama_preferences.items(), key=lambda x: x[0].lower())): |
| 100 | + plt.subplot(2, 2, idx + 1) |
| 101 | + plt.title(key) |
| 102 | + plt.imshow(rama_pref_values[key], cmap=rama_preferences[key]["cmap"], |
| 103 | + norm=colors.BoundaryNorm(rama_preferences[key]["bounds"], rama_preferences[key]["cmap"].N), |
| 104 | + extent=(-180, 180, 180, -180)) |
| 105 | + plt.scatter(normals[key]["x"], normals[key]["y"]) |
| 106 | + plt.scatter(outliers[key]["x"], outliers[key]["y"], color="red") |
| 107 | + plt.xlim([-180, 180]) |
| 108 | + plt.ylim([-180, 180]) |
| 109 | + plt.plot([-180, 180], [0, 0], color="black") |
| 110 | + plt.plot([0, 0], [-180, 180], color="black") |
| 111 | + plt.locator_params(axis='x', nbins=7) |
| 112 | + plt.xlabel(r'$\phi$') |
| 113 | + plt.ylabel(r'$\psi$') |
| 114 | + plt.grid() |
| 115 | + |
| 116 | + plt.tight_layout() |
| 117 | + # plt.savefig("asd.png", dpi=300) |
| 118 | + plt.show() |
| 119 | + |
| 120 | + |
| 121 | +if __name__ == '__main__': |
| 122 | + if len(sys.argv) < 2: |
| 123 | + sys.exit("Usage: ramachandran_calc.py my_pdb_file.pdb") |
| 124 | + ramachandran_calc(sys.argv[1:]) |
0 commit comments