Analyzing the results

Filter and plot INPHARMA output

The following analysis depends on the R package spinpharmaR. See SpinpharmaR for more information.

## Load required libraries
library(spinpharmaR)

## Read raw data
data <- read.qual("out/qual.dat")

## Print a summary
data

## Sort based on Inpharma score (R)
inds <- order(data$qual[,"R"], decreasing=T)

## Display results
head( data$qual[inds, ] )
               File_1             File_2        R       Qf       Sf
9  glide_dock-001.pdb glide_dock-009.pdb 0.816929 0.633009 4.398266
3  glide_dock-001.pdb glide_dock-003.pdb 0.815495 0.618362 4.460892
10 glide_dock-001.pdb glide_dock-010.pdb 0.806422 0.593947 4.766797
17 glide_dock-001.pdb glide_dock-017.pdb 0.796391 0.614154 4.227981
8  glide_dock-001.pdb glide_dock-008.pdb 0.795783 0.597894 4.752138
32 glide_dock-002.pdb glide_dock-003.pdb 0.784610 0.630993 3.604029

The column-wise output corresponds to:

File_1 / File_2

Filenames of the protein-ligand complexes

R

Pearson correlation coefficient (experimental vs theoretical NOEs)

Qf

Quality factor

Sf

Scaling factor

For validation, it can be useful to plot the INPHARMA results vs the structural deviation to the X-ray structures. The following code will use Bio3D to calculate the root mean squared deviation (RMSD) of all the docking modes (as found in the .nam files), and add the RMSD data to the spINPHARMA output:

## Calculating RMSD requires package 'bio3d'
library(bio3d)

## Calculate RMSD
rmsd.a <- calc.rmsd("modes_LL1.nam.all", "../../02_structures/3DNE_LL1.pdb", 299)
rmsd.b <- calc.rmsd("modes_M77.nam.all", "../../02_structures/3DNE_M77.pdb", 299)

## Join INPHARMA and RMSD data
data <- join.rmsd(data, rmsd.a, rmsd.b)

The function filter.hits and plot.qual provides functionality to filter and visualize the INPHARMA results:

## Filter based on INPHARMA Score
hits <- filter.hits(data, r.top=50, qf = FALSE, sf = FALSE)

## Print summary of hits
hits

## Print sorted results
head( data$qual[hits$good.inds,] )

## Plot Inpharma Score (R) vs RMSD to X-ray structures
pdf("scatterplot_R.pdf", w=7, h=7)
plot.qual(data, hits,
          highlights=c("background", "good", "top"),
          ylim=c(-0.3, 1))

## and quality factor...
plot.qual(data, hits, plot.y = "Qf", ylab="Quality factor",
          highlights=c("background", "good", "top"))

dev.off()

Plot individual NOEs

To plot the individual NOEs the functions read.noes and plot.noes can be used:

## Read individual NOEs
noes <- read.noes("out/corr.dat", mixing.times=3)

## Print summary of NOE data
noes

## Plot NOEs of combination 1 
cf <- plot.noes(noes[[1]]$noes, noes$exp, plot.labels=FALSE)
mtext(noes[[1]]$names, side=3, at=0.0, adj=0, line=0.5)
mtext(paste("R", cf$R, sep="="), side=3, at=0.0, adj=0, line=-1.5)
abline(cf$lm)

## With labels
cf <- plot.noes(noes[[1]]$noes, noes$exp, plot.labels=TRUE)

## Plot second set of NOEs
plot.noes(noes[[2]]$noes, noes$exp)
mtext(noes[[2]]$names, side=3, at=0.0, adj=0, line=0.5)

## Plot 9th set of NOEs
plot.noes(noes[[9]]$noes, noes$exp)
mtext(noes[[9]]$names, side=3, at=0.0, adj=0, line=0.5)
dev.off()