do plot for ProMod3/ProModII comparison

modelling/do_plots_promodII.py

+import matplotlib.pyplot as plt
+import json
+import numpy as np
+import math
+
+
+promod_label = 'ProMod3'
+promod2_label = 'ProModII'
+
+promod_data_file = 'promod_scores.json'
+promod2_data_file = 'promodII_scores.json'
+
+
+plot_name = 'promod3_vs_promodII.png'
+
+cred = (128.0/255,0.0,0.0)
+cblue = (102.0/255,153.0/255,204.0/255)
+cgreen = (102.0/255,148.0/255,0.0)
+cpurple = (100.0/255,0.0,200.0/255)
+corange = (255.0/255,123.0/255,0.0)
+
+with open(promod_data_file) as fh:
+    promod_data = json.load(fh)
+with open(promod2_data_file) as fh:
+    promod2_data = json.load(fh)
+
+lddt_values_promod = list()
+lddt_values_promod2 = list()
+probity_values_promod = list()
+probity_values_promod2 = list()
+lddt_diffs = list()
+probity_diffs = list()
+keys = list()
+
+for key in promod_data:
+    if key in promod2_data:
+        lddt_values_promod.append(promod_data[key]['lddt'] * 100)
+        lddt_values_promod2.append(promod2_data[key]['lddt'] * 100)
+        lddt_diffs.append(lddt_values_promod[-1] - lddt_values_promod2[-1])
+        probity_values_promod.append(promod_data[key]['MolProbity score'])
+        probity_values_promod2.append(promod2_data[key]['MolProbity score'])
+        probity_diffs.append(probity_values_promod[-1] - probity_values_promod2[-1])
+        keys.append(key)
+
+
+fig, axs = plt.subplots(3, 2, figsize=(7,10.5))
+hist_ax = axs[0,0]
+
+# plot both in the same plot
+xs = np.linspace(-7.0, 7.0, 300)
+n_lddt, bins_lddt, patches_lddt = hist_ax.hist(lddt_diffs, 50, range=(-7.0,7.0), 
+                                               facecolor=cred, alpha=0.75, 
+                                               label='lDDT', linewidth=2.0, edgecolor='k')
+n_probity, bins_probity, patches_probity = hist_ax.hist(probity_diffs, 50,
+                                                        range=(-7.0,7.0), 
+                                                        facecolor=cblue, alpha=0.75, 
+                                                        label='MolProbity',
+                                                        linewidth=2.0,
+                                                        edgecolor='k')
+hist_ax.axvline(x=0.0,  linewidth=2, color='k', linestyle='--')
+
+hist_ax.set_title('a) Modelling Benchmark', loc='left', y=1.08, x=-0.11, fontsize='x-large')
+hist_ax.set_xlabel(r'$\Delta$ score (ProMod3 - ProModII)',fontsize='large')
+hist_ax.set_ylabel('N',fontsize='large')
+hist_ax.legend(frameon=False)
+
+probity_clash_promod = list()
+probity_clash_promod2 = list()
+probity_rotamer_outliers_promod = list()
+probity_rotamer_outliers_promod2 = list()
+probity_ramachandran_outliers_promod = list()
+probity_ramachandran_outliers_promod2 = list()
+keys = list()
+for key in promod_data:
+    if key in promod2_data:
+        probity_clash_promod.append(promod_data[key]['Clashscore'])
+        probity_clash_promod2.append(promod2_data[key]['Clashscore'])
+        probity_rotamer_outliers_promod.append(promod_data[key]['Rotamer outliers'])
+        probity_rotamer_outliers_promod2.append(promod2_data[key]['Rotamer outliers'])
+        probity_ramachandran_outliers_promod.append(promod_data[key]['Ramachandran outliers'])
+        probity_ramachandran_outliers_promod2.append(promod2_data[key]['Ramachandran outliers'])
+        keys.append(key)
+
+
+
+def DoThingsWithAxes(ax, x_values, y_values, title, xlabel, ylabel):
+    
+    ax.plot(x_values, y_values, '.', color = (128.0/255,0.0,0.0))
+    ax.plot([-1.0,1000.0], [-1.0,1000], color = 'k',linestyle='--')
+    ax.set_title(title, loc='left', y=1.08, x=-0.11, fontsize='x-large')
+    ax.set_xlabel(xlabel, fontsize='large')
+    ax.set_ylabel(ylabel, fontsize='large')
+    max_val = math.ceil(max([max(x_values), max(y_values)]))
+    ax.set_xlim([0, max_val])
+    ax.set_ylim([0, max_val])
+    ax.set_aspect('equal', 'box')
+
+    tick_locations = list()
+    step_size = None
+    if max_val <= 5:
+      step_size = 1
+    elif max_val <= 14:
+      step_size = 2
+    elif max_val <= 30:
+      step_size = 5
+    elif max_val <= 100:
+      step_size = 20
+    else:
+      step_size = 50
+    for i in range(0, int(max_val)+step_size, step_size):
+      tick_locations.append(i)
+
+    ax.set_xticks(tick_locations) 
+    ax.set_yticks(tick_locations)
+   
+
+lddt_ax = axs[0, 1]
+probity_overall_ax = axs[1, 0]
+probity_clash_ax = axs[1, 1]
+probity_rotamer_ax = axs[2, 0]
+probity_ramachandran_ax = axs[2, 1]
+
+
+DoThingsWithAxes(lddt_ax, lddt_values_promod, 
+                 lddt_values_promod2, 'b) lDDT',
+                 promod_label, promod2_label)
+
+DoThingsWithAxes(probity_overall_ax, probity_values_promod, 
+                 probity_values_promod2, 'c) MolProbity Overall',
+                 promod_label, promod2_label)
+
+DoThingsWithAxes(probity_clash_ax, probity_clash_promod, 
+                 probity_clash_promod2, 'd) MolProbity Clash',
+                 promod_label, promod2_label)
+
+DoThingsWithAxes(probity_rotamer_ax, probity_rotamer_outliers_promod, 
+                 probity_rotamer_outliers_promod2, 'e) MolProbity Rot. Outliers',
+                 promod_label, promod2_label)
+
+DoThingsWithAxes(probity_ramachandran_ax, probity_ramachandran_outliers_promod, 
+                 probity_ramachandran_outliers_promod2, 
+                 'f) MolProbity Ram. Outliers', promod_label, promod2_label)
+
+plt.tight_layout(pad=1.2, h_pad=1.5, w_pad=1.5, rect=None)
+plt.savefig(plot_name, dpi=300)
+
+print('avg. lddt value', promod_label, np.mean(lddt_values_promod))
+print('avg. lddt value', promod2_label, np.mean(lddt_values_promod2))
+print('diff avg lddt value',np.mean(lddt_values_promod)-np.mean(lddt_values_promod2))
+print('avg. probity value', promod_label, np.mean(probity_values_promod))
+print('avg. probity value', promod2_label, np.mean(probity_values_promod2))
+print('diff avg probity value', np.mean(probity_values_promod)-np.mean(probity_values_promod2))
+print('avg. Molprobity clash', promod_label, np.mean(probity_clash_promod))
+print('avg. Molprobity clash', promod2_label, np.mean(probity_clash_promod2))
+print('avg. Molprobity rotamer outliers', promod_label, 
+      np.mean(probity_rotamer_outliers_promod))
+print('avg. Molprobity rotamer outliers', promod2_label, 
+      np.mean(probity_rotamer_outliers_promod2))
+print('avg. Ramachandran outliers', promod_label, 
+      np.mean(probity_ramachandran_outliers_promod))
+print('avg. Ramachandran outliers', promod2_label, 
+      np.mean(probity_ramachandran_outliers_promod2))
+