# -*- tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
# vi: set ts=2 noet:
#
# (c) Copyright Rosetta Commons Member Institutions.
# (c) This file is part of the Rosetta software suite and is made available under license.
# (c) The Rosetta software is developed by the contributing members of the Rosetta Commons.
# (c) For more information, see http://www.rosettacommons.org. Questions about this can be
# (c) addressed to University of Washington UW TechTransfer, email: license@u.washington.edu.
require(polynom)
valid_rosetta_database_path <- function(rosetta_database_path) {
if(is.null(rosetta_database_path)){
rosetta_database_path <- Sys.getenv("ROSETTA3_DB")
}
if(is.null(rosetta_database_path) || !file.exists(rosetta_database_path)){
stop("Please specify the rosetta_database path via the --database command line argument or set the $ROSETTA3_DB environment variable.")
}
rosetta_database_path
}
#It's necessary to actually go into the folder because schema.sql
#tries to import the .csv files with in the current directory.
convert_parameter_set_into_database <- function(
rosetta_database_path, parameter_set, debug=FALSE){
parameter_set_path <- paste(rosetta_database_path, "scoring", "score_functions", "hbonds", parameter_set, sep="/")
if(!file.exists(parameter_set_path)){
stop(paste("Attempting to convert the hbond parameter set '",
parameter_set_path, "' ",
"set into a database, but it does not exist.", sep=""))
}
cwd <- getwd()
setwd(parameter_set_path)
system("sqlite3 params.db3 < schema.sql")
setwd(cwd)
parameter_set_path
}
################# Polynomial Functions ###################
# return the HBPoly1D table
get_polynomials <- function(rosetta_database_path, parameter_set, debug=FALSE){
con <- dbConnect("SQLite",
paste(rosetta_database_path,
"scoring/score_functions/hbonds", parameter_set, "params.db3", sep="/"))
poly_params <- dbGetQuery(con, "SELECT * FROM HBPoly1D;")
dbDisconnect(con)
if(debug){
cat("poly_params:\n")
print(poly_params)
}
poly_params
}
#return a polynom::polynomial object that represents the polynomial
#with the specified name in the specified polynomials data.frame.
get_1d_polynomial <- function(polynomials, name){
df <- polynomials[polynomials$name==name,]
if (nrow(df) !=1){
stop(paste("Polynomial, ", name,", could not be found.", sep=""))
}
#we don't know the order in which the rows columns selected will be returned!
switch(df$degree,
p <- df[1,c("c_a")],
p <- df[1,c("c_b","c_a")],
p <- df[1,c("c_c","c_b","c_a")],
p <- df[1,c("c_d","c_c","c_b","c_a",)],
p <- df[1,c("c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_j","c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_k","c_j","c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")])
polynomial(p)
}
get_1d_polynomial_from_types <- function(polynomials, sample_source, acc_chem_type, don_chem_type, separation, dimension){
df <- polynomials[
polynomials$sample_source == as.character(sample_source) &
polynomials$don_chem_type == as.character(don_chem_type) &
polynomials$acc_chem_type == as.character(acc_chem_type) &
polynomials$separation == as.character(separation) &
polynomials$dimension == as.character(dimension),]
if(nrow(df) != 1){
print("ERROR: in get_1d_polynomial_from_types")
print(paste(sample_source, don_chem_type, acc_chem_type, separation, dimension))
print(polynomials)
print(df)
}
switch(df$degree,
p <- df[1,c("c_a")],
p <- df[1,c("c_b","c_a")],
p <- df[1,c("c_c","c_b","c_a")],
p <- df[1,c("c_d","c_c","c_b","c_a",)],
p <- df[1,c("c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_j","c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")],
p <- df[1,c("c_k","c_j","c_i","c_h","c_g","c_f","c_e","c_d","c_c","c_b","c_a")])
polynomial(p)
}
polynomial_range <- function(polynomials, name){
if(!(name %in% polynomials$name)){
stop(paste("Polynomial name, ", name,", is not recognized.", sep=""))
}
c(polynomials[polynomials$name==name, "xmin"], polynomials[polynomials$name==name, "xmax"])
}
display_polynomials <- function(parameter_set_path, polynomials){
polynomials <- adply(polynomials, 1, function(polynomial){
name <- polynomial[1, "name"]
polynomial$equation <- as.character(signif(get_1d_polynomial(polynomials,name), digits=4), decreasing=TRUE)
polynomial
})
cat("Plynomials defined in the parameter set '", parameter_set_path,"'\n", sep="")
print(polynomials[,c("name", "equation")], right=FALSE)
}
evaluate_polynomial <- function(params, name){
print(params)
print(class(params))
poly <- params[params$name == name,]
x <- seq(poly$xmin, poly$xmax, length.out=100)
p <- get_1d_polynomial(params, name)
data.frame(name=name, x=x, y=aaply(x, 1, as.function(p)))
}
################# Fade Interval Functions ###################
# return the HBFadeIntervals table
get_fade_intervals <- function(parameter_set_path, debug=FALSE){
con <- dbConnect("SQLite", paste(parameter_set_path, "/params.db3", sep=""))
fade_params <- dbGetQuery(con, "SELECT * FROM HBFadeIntervals;")
dbDisconnect(con)
if(debug){
cat("fade_params:\n")
print(fade_params)
}
fade_params
}
# the defining points of the piecewise linear for the fade interval
fade_interval_knots <- function(fade_intervals, name, lim){
df <- fade_intervals[fade_intervals$name==name,]
if (nrow(df) != 1){
print(fade_intervals)
stop(paste("Fade interval, ", name, ", could not be found.", sep=""))
}
if(name=="fade_zero"){
return(data.frame(
name=name,
x=c(lim[1],(3*lim[1]+lim[2])/4,(lim[1]+3*lim[2])/4,lim[2]),
y=c(0,0,0,0)))
}
data.frame(
name=name,
x=c(df[1,"min0"], df[1,"fmin"], df[1,"fmax"], df[1,"max0"]),
y=c(0,1,1,0))
}
# return function that evaluates the fade value
# see core::scoring::hbonds::FadeInterval::value()
fade_function <- function(fade_interval_knots, x, smooth=FALSE){
knots <- fade_interval_knots[,"x"]
y = NA
if(fade_interval_knots[1,"name"]=="fade_zero") return(0)
if (x <= knots[3]){
if(x <= knots[1]){
y = 0
} else if(x >= knots[2]) {
y = 1
} else {
z= (x - knots[1])/(knots[2]-knots[1])
if(smooth){
y = -2*z^3 + 3*z^2
} else {
y = z
}
}
} else {
if(x >= knots[4]){
y = 0
} else {
z = (x-knots[3])/(knots[4]-knots[3])
if(smooth){
y = 2*z^3 - 3*z^2 + 1
} else {
y = 1-z
}
}
}
# cat("fade_function: (", knots[1], " ", knots[2], " ", knots[3], " ", knots[4], ") x=", x, " y=", y, "\n")
y
}
evaluate_fade_interval <- function(params, xmin, xmax, name){
x <- seq(xmin, xmax, length.out=100)
cat("evaluating fade interval for ", name,"\n")
print(params)
p <- data.matrix(params[,c("smooth", "min0", "fmin", "fmax", "max0")])
data.frame(name=name, x=x, y=aaply(x, 1, eval_fade, p))
}
################# HBond Parameter Functions ###################
get_polynomial_parameters <- function(
parameter_set_path, don_chem_type, acc_chem_type, separation, debug=FALSE){
con <- dbConnect("SQLite", paste(parameter_set_path, "/params.db3", sep=""))
poly_params <- function(which_poly){
p <- dbGetQuery(con, paste("
SELECT
p.degree, p.xmin, p.xmax,
p.c_a, p.c_b, p.c_c, p.c_d, p.c_e, p.c_f, p.c_g, p.c_h, p.c_j, p.c_k
FROM
HBPoly1D as p,
HBEval as e
WHERE
e.don_chem_type='", don_chem_type, "' AND
e.acc_chem_type='", acc_chem_type, "' AND
e.separation='", separation, "' AND
e.", which_poly, "= p.name;", sep=""))
p$name <- which_poly
p
}
params <- rbind(
poly_params("AHdist"),
poly_params("cosBAH_short"),
poly_params("cosBAH_long"),
poly_params("cosAHD_short"),
poly_params("cosAHD_long"))
dbDisconnect(con)
params
}
get_fade_parameters <- function(
parameter_set_path, don_chem_type, acc_chem_type, separation, debug=FALSE){
con <- dbConnect("SQLite", paste(parameter_set_path, "/params.db3", sep=""))
fade_params <- function(which_fade){
f <- dbGetQuery(con, paste("
SELECT
CASE f.junction_type WHEN 'smooth' THEN 1 ELSE 0 END AS smooth,
f.min0, f.fmin, f.fmax, f.max0
FROM
HBFadeIntervals as f,
HBEval as e
WHERE
e.don_chem_type='", don_chem_type, "' AND
e.acc_chem_type='", acc_chem_type, "' AND
e.separation='", separation, "' AND
e.", which_fade, "= f.name;", sep=""))
f$which_fade = which_fade
f
}
params <- rbind(
fade_params("AHdist_short_fade"),
fade_params("AHdist_long_fade"),
fade_params("cosBAH_fade"),
fade_params("cosAHD_fade"))
dbDisconnect(con)
params
}
# the donor chemical type, acceptor chemical type and sequence
# separation uniquely identify hydrogen bond evaluation classes.
# return a row from HBEval for the specified class.
get_evaluation_type <- function(
parameter_set_path, don_chem_type, acc_chem_type, separation){
con <- dbConnect("SQLite", paste(parameter_set_path, "/params.db3", sep=""))
eval_type <- dbGetQuery(con, paste("
SELECT * FROM HBEval
WHERE
don_chem_type='", don_chem_type, "' AND
acc_chem_type='", acc_chem_type, "' AND
separation='", separation, "';", sep=""))
# now do error checking...
if(nrow(eval_type) != 1){
# check which part of the key is causing problems
don_row <- dbGetQuery(con, paste("
SELECT * FROM HBDonChemType
WHERE don_chem_type='",don_chem_type,"';", sep=""))
if(nrow(don_row) != 1){
stop(paste("The donor chemical type,",
don_chem_type, ", is not recognized."))
}
acc_row <- dbGetQuery(con, paste("
SELECT * FROM HBAccChemType
WHERE acc_chem_type='",acc_chem_type,"';", sep=""))
if(nrow(acc_row) != 1){
stop(paste("The acceptor chemical type,",
acc_chem_type, ", is not recognized."))
}
sep_row <- dbGetQuery(con, paste("
SELECT * FROM HBSeqSep
WHERE name='",separation,"';", sep=""))
if(nrow(sep_row) != 1){
stop(paste("The sequence separation type,",
separation, ", is not recognized."))
}
stop(paste("Unable to find evaluation type for:
don_chem_type: '", don_chem_type, "'
acc_chem_type: '", acc_chem_type, "'
separation: '", separation, "'", sep=""))
}
dbDisconnect(con)
eval_type
}
get_evaluation_function <- function(
polynomials, fade_intervals, eval_type){
AHdist_poly <- get_1d_polynomial(polynomials, eval_type[1,"AHdist"])
cosBAH_short_poly <- get_1d_polynomial(polynomials, eval_type[1,"cosBAH_short"])
cosBAH_long_poly <- get_1d_polynomial(polynomials, eval_type[1,"cosBAH_long"])
cosAHD_short_poly <- get_1d_polynomial(polynomials, eval_type[1,"cosAHD_short"])
cosAHD_long_poly <- get_1d_polynomial(polynomials, eval_type[1,"cosAHD_long"])
AHdist_range = polynomial_range(polynomials, eval_type[1, "AHdist"])
cosBAH_range = polynomial_range(polynomials, eval_type[1, "cosBAH_short"])
cosAHD_range = polynomial_range(polynomials, eval_type[1, "cosAHD_short"])
AHdist_short_fade <- fade_interval_knots(
fade_intervals, eval_type[1,"AHdist_short_fade"], AHdist_range)
AHdist_long_fade <- fade_interval_knots(
fade_intervals, eval_type[1,"AHdist_long_fade"], AHdist_range)
cosBAH_fade <- fade_interval_knots(
fade_intervals, eval_type[1, "cosBAH_fade"], cosBAH_range)
cosAHD_fade <- fade_interval_knots(
fade_intervals, eval_type[1, "cosAHD_fade"], cosAHD_range)
function(AHdist, cosBAH, cosAHD){
result = 0
if(AHdist > AHdist_range[1] && AHdist < AHdist_range[2]){
result = result +
as.function(AHdist_poly)(AHdist) *
fade_function(cosBAH_fade, cosBAH, opt$smooth) *
fade_function(cosAHD_fade, cosAHD, opt$smooth)
}
if(cosBAH > cosBAH_range[1] && cosBAH < cosBAH_range[2]){
result = result +
fade_function(AHdist_short_fade, AHdist, opt$smooth) *
as.function(cosBAH_short_poly)(cosBAH) *
fade_function(cosAHD_fade, cosAHD, opt$smooth)
result = result +
fade_function(AHdist_long_fade, AHdist, opt$smooth) *
as.function(cosBAH_long_poly)(cosBAH) *
fade_function(cosAHD_fade, cosAHD, opt$smooth)
}
if(cosAHD > cosAHD_range[1] && cosAHD < cosAHD_range[2]){
result = result +
fade_function(AHdist_short_fade, AHdist, opt$smooth) *
fade_function(cosBAH_fade, cosBAH, opt$smooth) *
as.function(cosAHD_short_poly)(cosAHD)
result = result +
fade_function(AHdist_long_fade, AHdist, opt$smooth) *
fade_function(cosBAH_fade, cosBAH, opt$smooth) *
as.function(cosAHD_long_poly)(cosAHD)
}
result
}
}
evaluate_function <- function(evaluation_function, AHdist, cosBAH, cosAHD){
ddply(expand.grid(list(AHdist=AHdist, cosBAH=cosBAH, cosAHD=cosAHD)),
.variables=c("AHdist", "cosBAH", "cosAHD"), function(df){
data.frame(energy=evaluation_function(
df[1,"AHdist"], df[1,"cosBAH"], df[1, "cosAHD"]))
})
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.