Analyzing the results

From INPHARMA wiki
Jump to: navigation, search

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()


Scatterplot R LL1-M77.png

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)

Noes LL1-M77.png


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

Noes LL1-M77-labels.png


## 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()