R/meltR.a.R
In JPSieg/MeltR: Automated fitting of RNA/DNA absorbance melting curves and fluorescence binding isotherms in R
Defines functions
meltR.A
Documented in
meltR.A
#'Fits melting curves to obtain thermodynamic parameters
#'
#'Automates the trivial but time-consuming tasks associated with non-linear regression analysis of melting curves.
#'Calculates extinction coefficients, subtracts out the baseline buffer readings, and
#'calculates the strand concentration, Ct, in each sample. Then uses three non-linear regression
#'methods to calculate thermodynamic parameters. Method 1 fits each melting curve individually
#'then reports the average dH and dS from all of the curves. Method 2 calculates the Tm for each
#'melting curve, and calculates thermodynamic parameters by fitting the relationship between Tm
#'and Ct. Method 3 calculates thermodynamic parameters with a global fit, where H and S are constant
#'between isotherms and the baselines are allowed to float.
#'
#'@param data_frame data_frame containing absorbance melting data. Requires Sample, Pathlength, Temperature, and Absorbance columns.
#'@param blank The blank sample for background subtraction, or a list of blanks to apply to different samples for background subtraction. "none" to turn off background subtraction. If there is a single blank in the data set, the identity of the blank, for example, blank = 1 or blank = "blank". If there are multiple blanks in the data, blank = list(c("Sample 1", "Blank 1"), c("Sample 2", "Blank 2")) and so on. Sample identifiers should be what they are in the data frame. If you need to figure out what the sample identifiers are, run unique(df$Sample), where df is the name of the R data frame you are using.
#'@param NucAcid A vector containing the Nucleic acid type and the sequences you are fitting for calculating extinction coefficients. Examples: c("RNA", "UUUUUU", "AAAAAA"), c("DNA", "GCTAGC"), etc... . For a custom extinction coefficient enter "Custom" followed by the molar extinction coefficients for every nucleic acid in the sample. For example, c("Custom", 10000, 20000). For non-absorbance melts, one can supply concentrations directly instead of extinction coefficients. For example, c("Concentration", 6.90e-06, 1.15e-05, 1.81e-05, 2.86e-05, 4.99e-05, 7.62e-05, 1.35e-04, 2.19e-04, 3.50e-04).
#'@param wavelength The wavelength you are using in the data set in nm. Options for RNA: 300, 295, 290, 285, 280, 275, 270, 265, 260, 255, 250, 245, 240, 235, and 230 nm. Options for DNA 260 nm. Most accurate at pH 7.0.
#'@param Mmodel The molecular model you want to fit. Options: "Monomolecular.2State", "Monomolecular.3State", "Heteroduplex.2State", "Homoduplex.2State".
#'@param Tmodel The thermodynamic model you want to fit. Options: "VantHoff". Default = "VantHoff".
#'@param concT The temperature used to calculate the concentration with Beers law. Default = 90.
#'@param outliers A vector containing the identifiers of the outlier samples that you want to remove. If you need to figure out what the sample identifiers are, run unique(df$Sample), where df is the name of the R data frame you are using.
#'@param fitTs Option to only fit certain temperature ranges for melting curves. Either a vector or a list. If this is set to a vector, meltR.A will only fit temperatures in this range for all melting curves, example = c(17, 75). If set to a list of vectors, meltR.A will change what values are fit for each curve. Example = list(c(0,100), c(17,75), .... , c(0,100)). The length of this list has to be the equal to the number of samples that will be fit. The list should not include vectors for blanks.
#'@param methods What methods do you want to use to fit data. Default = c(TRUE, TRUE, TRUE). Note, method 1 must be set to TRUE or the subsequent steps will not work. Set to c(TRUE, FALSE, FALSE) to only use method 1.
#'@param Tm_method Either "nls" to use the Tms from the fits in Method 1, "lm" to use a numeric method based on linear regression of fraction unfolded calculated with method 1, or "polynomial" to calculate Tms using the first derivative of a polynomial that approximates each curve.
#'@param Weight_Tm_M2 If you use the "nls" method for Tm method, this option can turn on weighted non-linear regression for method 2. If TRUE, method 2 will use the port algorithm to weight the regression in method 2 to standard errors in the Tm determined with method 1. Set to FALSE by default.
#'@param Save_results What results to save. Options: "all" to save PDF plots and ".csv" formatted tables of parameters, "some" to save ".csv" formatted tables of parameters, or "none" to save nothing.
#'@param file_prefix Prefix that you want on the saved files.
#'@param file_path Path to the directory you want to save results in.
#'@param auto.trimmed Ignore this argument unless you are writing auto baseline trimmers
#'@param Silent TRUE to not print data in your console. Default = FALSE.
#'@return A meltR.A fit object containing a list of data objects containing raw data, data transformations, fit objects, and statistics from the fits for plotting, exporting, and other advanced analysis.
#' \itemize{
#' \item 1. Summary - A data frame containing the thermodynamic parameters from each method.
#' \item 2. Method.1.indvfits - A data frame containing the thermodynamic parameters from the individual fits.
#' \item 3. Range - A data frame containing maximum %error between Method 1, 2, and 3 for each thermodynamic parameter.
#' \item 4. Derivatives.data - A data frame containing the first and second derivatives for each sample.
#' \item 5. Method.1.data - A data frame containing the raw data from method 1 and the model.
#' \item 6. Method.1.fit - A list of nls objects containing the fits obtained from fitting melting curves individually. Fit statistics, such as residuals and covariance, can be extracted here.
#' \item 7. Method.2.data - A data frame containing the raw data from method 2 and the model.
#' \item 8. Method.2.fit - A nls object containing the fit obtained from fitting the relationship of Tm and Ct. Fit statistics, such as residuals and covariance, can be extracted here.
#' \item 9. Method.3.data - A data frame containing the raw data from method 3 and the model.
#' \item 10. Method.3.fit - A nls object containing the fit obtained from fitting the raw data. Fit statistics, such as residuals and covariance, can be extracted here.
#' }
#' @export
meltR.A = function(data_frame,
blank = "none",
NucAcid,
wavelength = 260,
concT = 90,
outliers = NA,
fitTs = NULL,
methods = c(TRUE, TRUE, TRUE),
Tm_method = "nls",
Weight_Tm_M2 = F,
Mmodel,
Tmodel = "VantHoff",
Save_results = "none",
file_prefix = "Fit",
file_path = getwd(),
auto.trimmed = FALSE,
Silent = FALSE) {
####List of molecular models to fit####
Mmodel_names <- c("Monomolecular.2State",
"Monomolecular.3State",
"Heteroduplex.2State",
"Homoduplex.2State")
Mmodels <- list(function(K){ (K/(1+K)) },
function(K1, K2, Ct){ 1/(1 + K1 + (K1*K2)) },
function(K, Ct){ ( (2/(K*Ct)) + 2 - sqrt(((((2/(K*Ct)) + 2)^2) - 4)))/2 },
function(K, Ct){ ((1/(2*K*Ct)) + 2 - sqrt(((((1/(2*K*Ct)) + 2)^2) - 4)))/2 })
names(Mmodels) <- Mmodel_names
####List of thermodynamics models to fit to####
Tmodel_names <- c("VantHoff")
Tmodels <- list(function(H, Tm, Temperature){(H/0.0019872)*((1/(Tm + 273.15)) - (1/(Temperature + 273.15)))})
names(Tmodels) <- Tmodel_names
####Assemble the models####
G_VH = function(H, S, Temperature){exp((S/0.0019872) - ((1/((Temperature + 273.15)*0.0019872))*H))}
if (Mmodel == "Monomolecular.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K = K)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Monomolecular.2State(K = K)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm){ (H/(273.15 + Tm)) }
calcG = function(H, Tm){H - (310*((H/(273.15 + Tm))))}
calcTM = function(H, S, Ct){(H/(S)) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(1)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
if (Mmodel == "Monomolecular.3State"){
if (Tmodel == "VantHoff"){
Model = function(H1, Tm1, H2, Tm2, mED, bED, EI, mESS, bESS, Temperature, Ct){
K1 <- exp(Tmodels$VantHoff(H = H1, Tm = Tm1, Temperature = Temperature))
K2 <- exp(Tmodels$VantHoff(H = H2, Tm = Tm2, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K1 = K1, K2 = K2, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*K1*K2*ED) + (f*K1*EI) + (f)*ESS
return(model)
}
GModel = function(H1, Tm1, H2, Tm2, mED, bED, EI, mESS, bESS, Sample, Temperature, Ct){
K1 <- exp(Tmodels$VantHoff(H = H1, Tm = Tm1, Temperature = Temperature))
K2 <- exp(Tmodels$VantHoff(H = H2, Tm = Tm2, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K1 = K1, K2 = K2, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*K1*K2*ED) + (f*K1*EI[Sample]) + (f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) }
calcTM = function(H, S, Ct){NA}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
NA
}
}
}
if (Mmodel == "Heteroduplex.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature, Ct){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature) + log(4/Ct))
f <- Mmodels$Heteroduplex.2State(K = K, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
TmModel = function(H, S, lnCt){
((0.0019872/H)*lnCt) + ((S - 0.0019872*log(4))/H)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Heteroduplex.2State(K = K, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) + (0.0019872*log(4/Ct)) }
calcG = function(H, Tm, Ct){ H - (310.15*((H/(273.15 + Tm)) + (0.0019872*log(4/Ct)))) }
calcTM = function(H, S, Ct){(H/((S) - 0.0019872*log(4/Ct))) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(4/Ct)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
if (Mmodel == "Homoduplex.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature, Ct){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature) + log(1/Ct))
f <- Mmodels$Homoduplex.2State(K = K, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
TmModel = function(H, S, lnCt){
((0.0019872/H)*lnCt) + (S/H)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Homoduplex.2State(K = K, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) + (0.0019872*log(1/Ct)) }
calcG = function(H, Tm, Ct){ H - (310.15*((H/(273.15 + Tm)) + (0.0019872*log(1/Ct)))) }
calcTM = function(H, S, Ct){(H/((S) - 0.0019872*log(1/Ct))) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(1/Ct)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
####Subtract out the blank####
if (is.list(blank)){ #routine for multiple blanks in a data set
df.samples <- {}
df.blank <- {}
for (i in c(1:length(blank))){
df.samples[[i]] = subset(data_frame, Sample == blank[[i]][1])
df.blank[[i]] = subset(data_frame, Sample == blank[[i]][2])
df.samples[[i]]$Absorbance = df.samples[[i]]$Absorbance - df.blank[[i]]$Absorbance
}
k <- df.samples[[1]]
if (length(blank) > 1){
for (i in c(2:length(df.samples))){
k <- rbind(k, df.samples[[i]])
}
}
no.background = k
}else{
if (blank == "none"){ #routine for no blank in the data set
no.background = data_frame
}else{ #routine for a single blank in the data set
samples <- {}
for (i in c(1:length(unique(data_frame$Sample)))){
samples[[i]] <- subset(data_frame, Sample == unique(data_frame$Sample)[i])
df.blank = subset(data_frame, Sample == blank)
samples[[i]]$Absorbance = samples[[i]]$Absorbance - (df.blank$Absorbance*samples[[i]]$Pathlength/df.blank$Pathlength)
}
k <- samples[[1]]
for (i in c(2:length(samples))){
k <- rbind(k, samples[[i]])
}
no.background <- subset(k, Sample != blank)
}
}
raw.df = no.background
####Calculate extinction coefficients####
if (NucAcid[1] == "Concentration"){
}else{
if (NucAcid[1] == "Custom"){
extcoef = c(sum(as.numeric(NucAcid[2:length(NucAcid)])), as.numeric(NucAcid[2:length(NucAcid)]))
names(extcoef) <- c("Total", as.character(2:length(NucAcid)))
}else{
tryCatch({
extcoef = MeltR::calc.extcoeff(NucAcid, wavelength, T)
},
error = function(e){print("There is no nucleotide extinction coefficient at this wavelength for a nucleotide you specified in your sequence. You will need to provide your own custom extinction coefficient")})
}
}
####Remove values not in fitTs####
if (is.null(fitTs) == FALSE){
if (is.list(fitTs) == FALSE){
ranges <- rep(list(fitTs), length(unique(no.background$Sample)))
}else{
ranges <- fitTs
}
a <- {}
for (i in c(1:length(unique(no.background$Sample)))){
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
if(length(which(a[[i]]$Temperature > ranges[[i]][2])) > 0){
a[[i]] <- a[[i]][ -which(a[[i]]$Temperature > ranges[[i]][2]) ,]
}
if(length(which(a[[i]]$Temperature < ranges[[i]][1])) > 0){
a[[i]] <- a[[i]][ -which(a[[i]]$Temperature < ranges[[i]][1]) ,]
}
}
b <- a[[1]]
if (length(a) > 1){
for (i in c(2:length(unique(no.background$Sample)))){
b <- rbind(b, a[[i]])
}
}
no.background <- b
}
####Calculate Ct for each curve####
if (NucAcid[1] == "Concentration"){
ct <- as.numeric(NucAcid[2:length(NucAcid)])
samples = {}
for (i in c(1:length(unique(no.background$Sample)))){
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
samples[[i]]$Ct <- ct[i]
}
k <- samples[[1]]
if (length(samples) > 1){
for (i in c(2:length(samples))){
#print(i)
k <- rbind(k, samples[[i]])
}
}
no.background <- k
}else{
samples <- {}
if (length(extcoef) == 3){
ct <- c()
for (i in c(1:length(unique(no.background$Sample)))){
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
df.raw = subset(raw.df, Sample == unique(no.background$Sample)[i])
if (is.atomic(extcoef)){
ct[i] <- 2*(df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef[[1]]*df.raw$Pathlength[1]))
}else{
ct[i] <- 2*(df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef$Total*df.raw$Pathlength[1]))
}
samples[[i]]$Ct <- ct[i]
plot(samples[[i]]$Temperature, 2*samples[[i]]$Absorbance/(samples[[i]]$Pathlength*samples[[i]]$Ct))
abline(h = extcoef$Total)
}
}
if (length(extcoef) == 2){
ct <- c()
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
df.raw = subset(data_frame, Sample == unique(no.background$Sample)[i])
if (is.atomic(extcoef)){
ct[i] <- (df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef[[1]]*df.raw$Pathlength[1]))
}else{
ct[i] <- (df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef$Total*df.raw$Pathlength[1]))
}
samples[[i]]$Ct <- ct[i]
}
}
k <- samples[[1]]
if (length(samples) > 1){
for (i in c(2:length(samples))){
k <- rbind(k, samples[[i]])
}
}
no.background = k
}
####Remove outliers####
if (is.na(outliers)){
}else{
for (i in 1:length(outliers)){
no.background = subset(no.background, Sample != outliers[i])
}
}
####Calculate starting thermo parameters for nls####
first.derive <- {}
T0.5 <- c()
T0.75 <- c()
startH <- c()
for (i in c(1:length(unique(no.background$Sample)))){
first.derive[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
fit.Em = lm(Absorbance ~ poly(Temperature, 20, raw=TRUE),
data = first.derive[[i]])
string.dE.dT = gsub("NA", "0", gsub(", ", " + ", toString(paste(0:20, "*", coef(fit.Em), "*x^", -1:19, sep =""))))
string.dE.dT2 = gsub("NA", "0", gsub(", ", " + ", toString(paste(0:20, "*", c(0,0:19),"*", coef(fit.Em), "*x^", -2:18, sep =""))))
dE.dT = function(x){eval(parse(text = string.dE.dT))}
dE.dT2 = function(x){eval(parse(text = string.dE.dT2))}
T0.5[i] = seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)[which.max(dE.dT(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)))]
T0.75[i] = seq(min(first.derive[[i]]$Temperature) + 10, max(first.derive[[i]]$Temperature) - 10, length.out = 1000)[which.min(dE.dT2(seq(min(first.derive[[i]]$Temperature) + 10, max(first.derive[[i]]$Temperature) - 10, length.out = 1000)))]
if (Mmodel == "Monomolecular.2State"){
startH[i] <- -0.0032/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
if (Mmodel == "Heteroduplex.2State"){
startH[i] <- -0.007/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
if (Mmodel == "Homoduplex.2State"){
startH[i] <- -0.0044/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
#svg("~/Desktop/Melt_curve_derivative.svg")
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$Absorbance,
# xlab = Temperature ~ (degree ~ C), ylab = "Absorbance")
#lines(first.derive[[i]]$Temperature, predict(fit.Em), col = "red")
#lines(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000),
# 1.4 +10*dE.dT(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)), col = "red")
#lines(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000),
# 1.475+10*dE.dT2(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)), col = "blue")
#abline(v = T0.75[i])
#abline(v = T0.5[i])
#dev.off()
first.derive[[i]]$dA.dT = dE.dT(first.derive[[i]]$Temperature) + dE.dT(first.derive[[i]]$Temperature)*((first.derive[[i]]$Absorbance - predict(fit.Em))/first.derive[[i]]$Absorbance)
first.derive[[i]]$dA.dT2 = dE.dT2(first.derive[[i]]$Temperature) + dE.dT2(first.derive[[i]]$Temperature)*((first.derive[[i]]$Absorbance - predict(fit.Em))/first.derive[[i]]$Absorbance)
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$dA.dT)
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$dA.dT2)
}
####Calculate starting baseline values for nls####
a <-{}
uppbl_fit <- {}
lowbl_fit <- {}
bl.start.plot <- {}
for (i in c(1:length(unique(no.background$Sample)))){
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
upperbl <- data.frame(
"Absorbance" = a[[i]]$Absorbance[which(a[[i]]$Absorbance >= quantile(a[[i]]$Absorbance, probs = 0.75))],
"Temperature" = a[[i]]$Temperature[which(a[[i]]$Absorbance >= quantile(a[[i]]$Absorbance, probs = 0.75))]
)
lowbl <- data.frame(
"Absorbance" = a[[i]]$Absorbance[which(a[[i]]$Absorbance <= quantile(a[[i]]$Absorbance, probs = 0.25))],
"Temperature" = a[[i]]$Temperature[which(a[[i]]$Absorbance <= quantile(a[[i]]$Absorbance, probs = 0.25))]
)
uppbl_fit[[i]] <- lm(Absorbance ~ Temperature, data = upperbl)
lowbl_fit[[i]] <- lm(Absorbance ~ Temperature, data = lowbl)
}
####Method 1 fit each curve individually####
if (methods[1] == TRUE){
a <-{}
fit <- {}
indvfits.H <- c()
indvfits.S <- c()
indvfits.G <- c()
indvfits.Tm <- c()
indvfits.SE.Tm <- c()
bED <- c()
mED <- c()
bSS <- c()
mSS <- c()
Ind.model = Model
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
tryCatch({
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
if (is.list(auto.trimmed)){
fitstart = auto.trimmed[[i]]
}else{
fitstart <- list(H = startH[i], Tm = T0.5[i],
mED = lowbl_fit[[i]]$coefficients[2], bED = lowbl_fit[[i]]$coefficients[1],
mESS = uppbl_fit[[i]]$coefficients[2], bESS = uppbl_fit[[i]]$coefficients[1])
}
if (Mmodel == "Monomolecular.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Heteroduplex.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature, Ct),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Homoduplex.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature, Ct),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
mED[i] <- coef(fit[[i]])[3]
bED[i] <- coef(fit[[i]])[4]
mSS[i] <- coef(fit[[i]])[5]
bSS[i] <- coef(fit[[i]])[6]
indvfits.H[i] <- coef(fit[[i]])[1]
if (Mmodel == "Monomolecular.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2])
}
if (Mmodel == "Heteroduplex.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
}
if (Mmodel == "Homoduplex.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
}
indvfits.Tm[i] <- coef(fit[[i]])[2]
indvfits.SE.Tm[i] <- coef(summary(fit[[i]]))[2,2]
}, error = function(e){print(paste("There was a problem fitting Sample", a[[i]]$Sample[1], "with method 1") )})
}
indvfits <- data.frame("Sample" = unique(no.background$Sample),
"Ct" = unique(no.background$Ct),
"H" = round(indvfits.H, 2),
"S" = round(indvfits.S, 5),
"G" = round(indvfits.G, 2),
"Tm" = round(indvfits.Tm, 2))
indvfits.mean <- list(round(mean(indvfits$H), 2), round(sd(indvfits$H), 2), round(1000*mean(indvfits$S), 2), round(1000*sd(indvfits$S), 2), round(mean(indvfits.G), 2), round(sd(indvfits.G), 2))
names(indvfits.mean) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
indvfits.mean <- data.frame(indvfits.mean)
if (Mmodel == "Heteroduplex.2State"){
no.background$Ext <- 2*no.background$Absorbance/(no.background$Pathlength*no.background$Ct)
}else{
no.background$Ext <- no.background$Absorbance/(no.background$Pathlength*no.background$Ct)
}
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_1_raw_fit_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(no.background$Temperature, no.background$Absorbance,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbance",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, predict(fit[[i]]), col = "red")
}
dev.off()
pdf(paste(file_path, "/", file_prefix, "_method_1_normalized_fit_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(no.background$Temperature, no.background$Ext,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbtivity (1/M*cm)",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
if (Mmodel == "Heteroduplex.2State"){
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, (2*predict(fit[[i]])/(a[[i]]$Ct[1]*a[[i]]$Pathlength[1])), col = "red")
}
}else{
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, (predict(fit[[i]])/(a[[i]]$Ct[1]*a[[i]]$Pathlength[1])), col = "red")
}
}
dev.off()
}
}
if (methods[1] == FALSE){
indvfits <- data.frame("Sample" = NA,
"Ct" = NA,
"H" = NA,
"S" = NA,
"G" = NA,
"Tm" = NA)
indvfits.mean <- list(NA, NA, NA, NA, NA, NA)
names(indvfits.mean) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
indvfits.mean <- data.frame(indvfits.mean)
}
####Method 2 Tm vs Ct####
if (methods[2] == TRUE){
#Decide what Tm_method to use
if (Tm_method == "lm"){
a <- {}
Tm_range <- {}
Tm_fit <- {}
Tm <- c()
lnCt <- c()
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i]) #Pull out data
#calculate f at each temp
a[[i]]$f <- (a[[i]]$Absorbance - (coef(fit[[i]])[5]*a[[i]]$Temperature + coef(fit[[i]])[6]))/((coef(fit[[i]])[3]*a[[i]]$Temperature + coef(fit[[i]])[4]) - (coef(fit[[i]])[5]*a[[i]]$Temperature + coef(fit[[i]])[6]))
Tm_range[[i]] <- data.frame( #Pull out data where f >= 0.4 and f <= 0.6
"Temperature" = a[[i]]$Temperature,
"f" = a[[i]]$f
)
Tm_fit[[i]] <- lm(f~Temperature, data = Tm_range[[i]]) #fit data where f >= 0.4 and f <= 0.6
Tm[i] <- (0.5 - coef(Tm_fit[[i]])[1])/coef(Tm_fit[[i]])[2]
lnCt[i] <- log(a[[i]]$Ct[1])
}
Tm_data <- data.frame("lnCt" = lnCt, "Tm" = Tm)
Weight_Tm_M2 = FALSE
}
if (Tm_method == "nls"){
lnCt = log(indvfits$Ct)
Tm = indvfits$Tm
Tm_data = data.frame(lnCt, Tm)
}
if (Tm_method == "polynomial"){
lnCt = log(indvfits$Ct)
Tm = T0.5
Tm_data = data.frame(lnCt, Tm)
Weight_Tm_M2 = FALSE
}
Tm_data$invT <- 1/(Tm_data$Tm + 273.15)
if (Mmodel != "Monomolecular.2State"){
if (Mmodel != "Monomolecular.3State"){
Tm.model = TmModel
if (Weight_Tm_M2){ #Weighted non-linear regression
propagated.weights = indvfits.SE.Tm/((Tm_data$Tm+273.15)^2) #Propagate Tm SE
Tm_vs_lnCt_fit <- nls(invT ~ TmModel(H, S, lnCt),
data = Tm_data,
weights = 1/propagated.weights,
start = list(H = mean(indvfits$H), S = mean(indvfits$S)),
control = nls.control(warnOnly = TRUE))
}else{
Tm_vs_lnCt_fit <- nls(invT ~ TmModel(H, S, lnCt),
data = Tm_data,
start = list(H = mean(indvfits$H), S = mean(indvfits$S)),
control = nls.control(warnOnly = TRUE))
}
Tm_vs_lnCt <- list(round(coef(Tm_vs_lnCt_fit)[1], 2), round(summary(Tm_vs_lnCt_fit)$coefficients[1,2], 2),
round(1000*coef(Tm_vs_lnCt_fit)[2], 2), round(1000*summary(Tm_vs_lnCt_fit)$coefficients[2,2], 2),
round(coef(Tm_vs_lnCt_fit)[1] - (310.15*coef(Tm_vs_lnCt_fit)[2]), 2), round(sqrt((summary(Tm_vs_lnCt_fit)$coefficients[1,2])^2 + (310*summary(Tm_vs_lnCt_fit)$coefficients[2,2])^2 - 2*310*summary(Tm_vs_lnCt_fit)$cov.unscaled[1,2]*(summary(Tm_vs_lnCt_fit)$sigma^2)),2))
names(Tm_vs_lnCt) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Tm_vs_lnCt)
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_2_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
if (Weight_Tm_M2){
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8,
ylim = c(min(Tm_data$invT) - max(propagated.weights), max(Tm_data$invT) + max(propagated.weights)))
arrows(x = Tm_data$lnCt, x0 = Tm_data$lnCt, length = 0.05,
y = Tm_data$invT - propagated.weights, y0 = Tm_data$invT + propagated.weights,
code = 3, lwd=2, angle = 90)
lines(x = Tm_data$lnCt, predict(Tm_vs_lnCt_fit), col = "red")
}else{
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
lines(x = Tm_data$lnCt, predict(Tm_vs_lnCt_fit), col = "red")
}
dev.off()
}
}else{
Tm.model = NA
}
}else{
Tm.model = NA
}
if (Mmodel == "Monomolecular.2State"){
if (Mmodel == "Monomolecular.3State"){
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_2_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
lines(Tm_data$lnCt, predict(lm(Tm_data$invT ~ Tm_data$lnCt)), col = "red")
dev.off()
}
}}
}
if (methods[2] == FALSE){
Tm_vs_lnCt <- list(NA, NA,
NA, NA,
NA, NA)
names(Tm_vs_lnCt) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Tm_vs_lnCt)
}
####Method 3 Global fitting####
if (methods[3] == TRUE){
b <- data.frame("Sample" = c(), "Pathlength" = c(), "Temperature" = c(),
"Absorbance" = c(), "Ct" = c(), "Ext" = c())
for (i in c(1:length(unique(no.background$Sample)))){
a <- subset(no.background, Sample == unique(no.background$Sample)[i])
a$Sample <- i
b <- rbind(b, a)
}
gfit_data <- b
gfit_start = list(H = mean(indvfits$H), S = mean(indvfits$S), mED = mED, bED = bED, mESS = mSS, bESS = bSS)
if (Mmodel == "Monomolecular.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Heteroduplex.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Homoduplex.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_3_Gfit_raw_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(gfit_data$Temperature, gfit_data$Absorbance,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbance",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(unique(gfit_data$Sample)))){
a <- subset(gfit_data, Sample == unique(gfit_data$Sample)[i])
if (Mmodel == "Monomolecular.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature)))),
col = "red")
}
if (Mmodel == "Heteroduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1]),
col = "red")
}
if (Mmodel == "Homoduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1]),
col = "red")
}
}
dev.off()
pdf(paste(file_path, "/", file_prefix, "_method_3_Gfit_normalized_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(gfit_data$Temperature, gfit_data$Ext,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbtivity (1/M*cm)",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(unique(gfit_data$Sample)))){
a <- subset(gfit_data, Sample == unique(gfit_data$Sample)[i])
if (Mmodel == "Monomolecular.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))))/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
if (Mmodel == "Heteroduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), 2*GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1])/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
if (Mmodel == "Homoduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1])/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
}
dev.off()
}
Gfit_summary <- list(round(coef(gfit)[1], 2), round(summary(gfit)$coefficients[1,2], 2),
round(1000*coef(gfit)[2], 2), round(1000*summary(gfit)$coefficients[2,2], 2),
round(coef(gfit)[1] - (310.15)*coef(gfit)[2], 2), round(sqrt((summary(gfit)$coefficients[1,2])^2 + (310.15*summary(gfit)$coefficients[2,2])^2 - 2*310.15*summary(gfit)$cov.unscaled[1,2]*(summary(gfit)$sigma^2)), 2))
names(Gfit_summary) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Gfit_summary <- data.frame(Gfit_summary)
}
if (methods[3] == FALSE){
Gfit_summary <- list(NA, NA,
NA, NA,
NA, NA)
names(Gfit_summary) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Gfit_summary)
}
####Assemble Results####
indvfits$Ct = as.numeric(formatC(indvfits$Ct, format = "e", digits = 2))
indvfits$S = 1000*indvfits$S
colnames(indvfits) = c("Sample", "Ct", "dH", "dS", "dG", "Tm")
if (Mmodel == "Monomolecular.2State"){
comparison <- rbind(indvfits.mean, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "3 Global fit")), comparison)
row.names(comparison) <- c(1:2)
}
if (Mmodel == "Heteroduplex.2State"){
comparison <- rbind(indvfits.mean, Tm_vs_lnCt, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "2 Tm versus ln[Ct]", "3 Global fit")), comparison)
row.names(comparison) <- c(1:3)
}
if (Mmodel == "Homoduplex.2State"){
comparison <- rbind(indvfits.mean, Tm_vs_lnCt, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "2 Tm versus ln[Ct]", "3 Global fit")), comparison)
row.names(comparison) <- c(1:3)
}
####Expectation maximized Tm calculation####
Tm_at_0.1mM = c()
SE.Tm_at_0.1mM = c()
for (i in 1:nrow(comparison)){
if (i == 1){
df.covar = indvfits
udH = df.covar$dH - mean(df.covar$dH)
udS = (df.covar$dS - mean(df.covar$dS))/1000
sigma = sum(udH*udS)/nrow(df.covar)
dH = mean(df.covar$dH)
dS = mean(df.covar$dS)/1000
SE.dH = sd(df.covar$dH)
SE.dS = sd(df.covar$dS)/1000
}
if (i == 2){
if (Mmodel == "Monomolecular.2State"){
SE.dH = NA
SE.dS = NA
}else{
if (methods[2] == F){
sigma = NA
dH = NA
dS = NA
SE.dH = NA
SE.dS = NA
}else{
sigma = summary(Tm_vs_lnCt_fit)$cov.unscaled[1,2]*(summary(Tm_vs_lnCt_fit)$sigma^2)
dH = coef(Tm_vs_lnCt_fit)[1]
dS = coef(Tm_vs_lnCt_fit)[2]
SE.dH = coef(summary(Tm_vs_lnCt_fit))[1,2]
SE.dS = coef(summary(Tm_vs_lnCt_fit))[2,2]
}
}
}
if (i == 3){
if (methods[2] == F){
sigma = NA
dH = NA
dS = NA
SE.dH = NA
SE.dS = NA
}else{
sigma = summary(gfit)$cov.unscaled[1,2]*(summary(gfit)$sigma^2)
dH = coef(gfit)[1]
dS = coef(gfit)[2]
SE.dH = coef(summary(gfit))[1,2]
SE.dS = coef(summary(gfit))[2,2]
}
}
Tm_at_0.1mM[i] = round(calcTM(H = dH, S = dS,
Ct = 10^-4), digits = 2)
SE.Tm_at_0.1mM[i] = round(calcTM.SE(H = dH, S = dS,
SE.H = SE.dH, SE.S = SE.dS,
Ct = 10^-4, covar = sigma), digits = 2)
}
comparison$Tm_at_0.1mM = Tm_at_0.1mM
comparison$SE.Tm_at_0.1mM = SE.Tm_at_0.1mM
colnames(comparison) = c("Method", "dH", "SE.dH", "dS", "SE.dS", "dG", "SE.dG", "Tm_at_0.1mM", "SE.Tm_at_0.1mM")
if (Save_results != "none"){
write.table(comparison, paste(file_path, "/", file_prefix, "_summary.csv", sep = ""), sep = ",", row.names = FALSE)
}
if (Save_results != "none"){
write.table(indvfits, paste(file_path, "/", file_prefix, "_method_1_individual_fits.csv", sep = ""), sep = ",", row.names = FALSE)
}
####Assemble data for custom derivative plots####
df.deriv = first.derive[[1]]
if (length(first.derive) > 1){
for (i in 2:length(first.derive)){
#print(i)
df.deriv = rbind(df.deriv, first.derive[[i]])
}
}
####Assemble data for custom individual fit plots####
df.method.1 = subset(no.background, Sample == unique(no.background$Sample)[1])
df.method.1$Model = predict(fit[[1]])
if (length(unique(no.background$Sample)) >1){
for (i in 2:length(unique(no.background$Sample))){
df.method.1.i = subset(no.background, Sample == unique(no.background$Sample)[i])
df.method.1.i$Model = predict(fit[[i]])
df.method.1 = rbind(df.method.1, df.method.1.i)
}
}
####Assemble data for custom 1/Tm versus lnCt plots####
if (Mmodel != "Monomolecular.2State"){
if (methods[2] == TRUE){
Tm_data$Model = predict(Tm_vs_lnCt_fit)
}
}
####Assemble global fit data for custom fit plots####
if (methods[3] == TRUE){
gfit_data$Model = predict(gfit)
}
####Compile data for the BLTrimmer####
BLTrimmer = list(raw.df, #background subtracted raw data
Mmodel,
Tmodel)
####Make a list of MeltR A settings####
settings = list(data_frame, #1
blank, #2
NucAcid, #3
concT, #4
fitTs, #5
methods, #6
Mmodel, #7
Tmodel, #8
Save_results, #9
file_prefix, #10
file_path, #11
Weight_Tm_M2) #12
####Assemble final output####
if (Silent){}else{
print("Individual curves")
print("dH and dG are in kcal/mol and dS is in cal/mol/K. Tms are in deg Celsius")
print(indvfits[order(indvfits$Ct),])
print("Summary")
print("dH and dG are in kcal/mol and dS is in cal/mol/K. Tms are in deg Celsius")
print(comparison)
range <- data.frame("dH" = round(100*abs((range(comparison$dH)[1]-range(comparison$dH)[2])/mean(comparison$dH)), digits = 1),
"dS" = round(100*abs((range(comparison$dS)[1]-range(comparison$dS)[2])/mean(comparison$dS)), digits = 1),
"dG" = round(100*abs((range(comparison$dG)[1]-range(comparison$dG)[2])/mean(comparison$dG)), digits = 1))
print("Maximum %error across methods")
print(range)
}
output <- list("Summary" = comparison,
"Method.1.indvfits" = indvfits,
"Range" = range,
"Derivatives.data" = df.deriv,
"Method.1.data" = df.method.1,
"Method.1.fit" = fit,
"BLTrimmer.data" = BLTrimmer,
"meltR.A.settings" = settings)
if (Mmodel != "Monomolecular.2State"){
if (methods[2] == TRUE){
output$Method.2.data = Tm_data
output$Method.2.fit <- Tm_vs_lnCt_fit
}
}
if (methods[3] == TRUE){
output$Method.3.data = gfit_data
output$Method.3.fit <- gfit
}
output <- output
}
JPSieg/MeltR documentation built on Feb. 4, 2024, 7:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions meltR.A
Documented in meltR.A
#'Fits melting curves to obtain thermodynamic parameters
#'
#'Automates the trivial but time-consuming tasks associated with non-linear regression analysis of melting curves.
#'Calculates extinction coefficients, subtracts out the baseline buffer readings, and
#'calculates the strand concentration, Ct, in each sample. Then uses three non-linear regression
#'methods to calculate thermodynamic parameters. Method 1 fits each melting curve individually
#'then reports the average dH and dS from all of the curves. Method 2 calculates the Tm for each
#'melting curve, and calculates thermodynamic parameters by fitting the relationship between Tm
#'and Ct. Method 3 calculates thermodynamic parameters with a global fit, where H and S are constant
#'between isotherms and the baselines are allowed to float.
#'
#'@param data_frame data_frame containing absorbance melting data. Requires Sample, Pathlength, Temperature, and Absorbance columns.
#'@param blank The blank sample for background subtraction, or a list of blanks to apply to different samples for background subtraction. "none" to turn off background subtraction. If there is a single blank in the data set, the identity of the blank, for example, blank = 1 or blank = "blank". If there are multiple blanks in the data, blank = list(c("Sample 1", "Blank 1"), c("Sample 2", "Blank 2")) and so on. Sample identifiers should be what they are in the data frame. If you need to figure out what the sample identifiers are, run unique(df$Sample), where df is the name of the R data frame you are using.
#'@param NucAcid A vector containing the Nucleic acid type and the sequences you are fitting for calculating extinction coefficients. Examples: c("RNA", "UUUUUU", "AAAAAA"), c("DNA", "GCTAGC"), etc... . For a custom extinction coefficient enter "Custom" followed by the molar extinction coefficients for every nucleic acid in the sample. For example, c("Custom", 10000, 20000). For non-absorbance melts, one can supply concentrations directly instead of extinction coefficients. For example, c("Concentration", 6.90e-06, 1.15e-05, 1.81e-05, 2.86e-05, 4.99e-05, 7.62e-05, 1.35e-04, 2.19e-04, 3.50e-04).
#'@param wavelength The wavelength you are using in the data set in nm. Options for RNA: 300, 295, 290, 285, 280, 275, 270, 265, 260, 255, 250, 245, 240, 235, and 230 nm. Options for DNA 260 nm. Most accurate at pH 7.0.
#'@param Mmodel The molecular model you want to fit. Options: "Monomolecular.2State", "Monomolecular.3State", "Heteroduplex.2State", "Homoduplex.2State".
#'@param Tmodel The thermodynamic model you want to fit. Options: "VantHoff". Default = "VantHoff".
#'@param concT The temperature used to calculate the concentration with Beers law. Default = 90.
#'@param outliers A vector containing the identifiers of the outlier samples that you want to remove. If you need to figure out what the sample identifiers are, run unique(df$Sample), where df is the name of the R data frame you are using.
#'@param fitTs Option to only fit certain temperature ranges for melting curves. Either a vector or a list. If this is set to a vector, meltR.A will only fit temperatures in this range for all melting curves, example = c(17, 75). If set to a list of vectors, meltR.A will change what values are fit for each curve. Example = list(c(0,100), c(17,75), .... , c(0,100)). The length of this list has to be the equal to the number of samples that will be fit. The list should not include vectors for blanks.
#'@param methods What methods do you want to use to fit data. Default = c(TRUE, TRUE, TRUE). Note, method 1 must be set to TRUE or the subsequent steps will not work. Set to c(TRUE, FALSE, FALSE) to only use method 1.
#'@param Tm_method Either "nls" to use the Tms from the fits in Method 1, "lm" to use a numeric method based on linear regression of fraction unfolded calculated with method 1, or "polynomial" to calculate Tms using the first derivative of a polynomial that approximates each curve.
#'@param Weight_Tm_M2 If you use the "nls" method for Tm method, this option can turn on weighted non-linear regression for method 2. If TRUE, method 2 will use the port algorithm to weight the regression in method 2 to standard errors in the Tm determined with method 1. Set to FALSE by default.
#'@param Save_results What results to save. Options: "all" to save PDF plots and ".csv" formatted tables of parameters, "some" to save ".csv" formatted tables of parameters, or "none" to save nothing.
#'@param file_prefix Prefix that you want on the saved files.
#'@param file_path Path to the directory you want to save results in.
#'@param auto.trimmed Ignore this argument unless you are writing auto baseline trimmers
#'@param Silent TRUE to not print data in your console. Default = FALSE.
#'@return A meltR.A fit object containing a list of data objects containing raw data, data transformations, fit objects, and statistics from the fits for plotting, exporting, and other advanced analysis.
#' \itemize{
#' \item 1. Summary - A data frame containing the thermodynamic parameters from each method.
#' \item 2. Method.1.indvfits - A data frame containing the thermodynamic parameters from the individual fits.
#' \item 3. Range - A data frame containing maximum %error between Method 1, 2, and 3 for each thermodynamic parameter.
#' \item 4. Derivatives.data - A data frame containing the first and second derivatives for each sample.
#' \item 5. Method.1.data - A data frame containing the raw data from method 1 and the model.
#' \item 6. Method.1.fit - A list of nls objects containing the fits obtained from fitting melting curves individually. Fit statistics, such as residuals and covariance, can be extracted here.
#' \item 7. Method.2.data - A data frame containing the raw data from method 2 and the model.
#' \item 8. Method.2.fit - A nls object containing the fit obtained from fitting the relationship of Tm and Ct. Fit statistics, such as residuals and covariance, can be extracted here.
#' \item 9. Method.3.data - A data frame containing the raw data from method 3 and the model.
#' \item 10. Method.3.fit - A nls object containing the fit obtained from fitting the raw data. Fit statistics, such as residuals and covariance, can be extracted here.
#' }
#' @export
meltR.A = function(data_frame,
blank = "none",
NucAcid,
wavelength = 260,
concT = 90,
outliers = NA,
fitTs = NULL,
methods = c(TRUE, TRUE, TRUE),
Tm_method = "nls",
Weight_Tm_M2 = F,
Mmodel,
Tmodel = "VantHoff",
Save_results = "none",
file_prefix = "Fit",
file_path = getwd(),
auto.trimmed = FALSE,
Silent = FALSE) {
####List of molecular models to fit####
Mmodel_names <- c("Monomolecular.2State",
"Monomolecular.3State",
"Heteroduplex.2State",
"Homoduplex.2State")
Mmodels <- list(function(K){ (K/(1+K)) },
function(K1, K2, Ct){ 1/(1 + K1 + (K1*K2)) },
function(K, Ct){ ( (2/(K*Ct)) + 2 - sqrt(((((2/(K*Ct)) + 2)^2) - 4)))/2 },
function(K, Ct){ ((1/(2*K*Ct)) + 2 - sqrt(((((1/(2*K*Ct)) + 2)^2) - 4)))/2 })
names(Mmodels) <- Mmodel_names
####List of thermodynamics models to fit to####
Tmodel_names <- c("VantHoff")
Tmodels <- list(function(H, Tm, Temperature){(H/0.0019872)*((1/(Tm + 273.15)) - (1/(Temperature + 273.15)))})
names(Tmodels) <- Tmodel_names
####Assemble the models####
G_VH = function(H, S, Temperature){exp((S/0.0019872) - ((1/((Temperature + 273.15)*0.0019872))*H))}
if (Mmodel == "Monomolecular.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K = K)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Monomolecular.2State(K = K)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm){ (H/(273.15 + Tm)) }
calcG = function(H, Tm){H - (310*((H/(273.15 + Tm))))}
calcTM = function(H, S, Ct){(H/(S)) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(1)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
if (Mmodel == "Monomolecular.3State"){
if (Tmodel == "VantHoff"){
Model = function(H1, Tm1, H2, Tm2, mED, bED, EI, mESS, bESS, Temperature, Ct){
K1 <- exp(Tmodels$VantHoff(H = H1, Tm = Tm1, Temperature = Temperature))
K2 <- exp(Tmodels$VantHoff(H = H2, Tm = Tm2, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K1 = K1, K2 = K2, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*K1*K2*ED) + (f*K1*EI) + (f)*ESS
return(model)
}
GModel = function(H1, Tm1, H2, Tm2, mED, bED, EI, mESS, bESS, Sample, Temperature, Ct){
K1 <- exp(Tmodels$VantHoff(H = H1, Tm = Tm1, Temperature = Temperature))
K2 <- exp(Tmodels$VantHoff(H = H2, Tm = Tm2, Temperature = Temperature))
f <- Mmodels$Monomolecular.2State(K1 = K1, K2 = K2, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*K1*K2*ED) + (f*K1*EI[Sample]) + (f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) }
calcTM = function(H, S, Ct){NA}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
NA
}
}
}
if (Mmodel == "Heteroduplex.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature, Ct){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature) + log(4/Ct))
f <- Mmodels$Heteroduplex.2State(K = K, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
TmModel = function(H, S, lnCt){
((0.0019872/H)*lnCt) + ((S - 0.0019872*log(4))/H)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Heteroduplex.2State(K = K, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) + (0.0019872*log(4/Ct)) }
calcG = function(H, Tm, Ct){ H - (310.15*((H/(273.15 + Tm)) + (0.0019872*log(4/Ct)))) }
calcTM = function(H, S, Ct){(H/((S) - 0.0019872*log(4/Ct))) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(4/Ct)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
if (Mmodel == "Homoduplex.2State"){
if (Tmodel == "VantHoff"){
Model = function(H, Tm, mED, bED, mESS, bESS, Temperature, Ct){
K <- exp(Tmodels$VantHoff(H = H, Tm = Tm, Temperature = Temperature) + log(1/Ct))
f <- Mmodels$Homoduplex.2State(K = K, Ct = Ct)
ED <- mED*Temperature + bED
ESS <- mESS*Temperature + bESS
model <- (f*ED) + (1-f)*ESS
return(model)
}
TmModel = function(H, S, lnCt){
((0.0019872/H)*lnCt) + (S/H)
}
GModel = function(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct){
K <- G_VH(H = H, S = S, Temperature = Temperature)
f <- Mmodels$Homoduplex.2State(K = K, Ct = Ct)
ED <- mED[Sample]*Temperature + bED[Sample]
ESS <- mESS[Sample]*Temperature + bESS[Sample]
model <- (f*ED) + (1-f)*ESS
return(model)
}
calcS = function(H, Tm, Ct){ (H/(273.15 + Tm)) + (0.0019872*log(1/Ct)) }
calcG = function(H, Tm, Ct){ H - (310.15*((H/(273.15 + Tm)) + (0.0019872*log(1/Ct)))) }
calcTM = function(H, S, Ct){(H/((S) - 0.0019872*log(1/Ct))) - 273.15}
calcTM.SE = function(H, S, SE.H, SE.S, Ct, covar){
C = 0.0019872*log(1/Ct)
x = sqrt((SE.H/(S-C))^2 + (H*SE.S/(S - C)^2)^2 - 2*H*covar/((S-C)^3))
}
}
}
####Subtract out the blank####
if (is.list(blank)){ #routine for multiple blanks in a data set
df.samples <- {}
df.blank <- {}
for (i in c(1:length(blank))){
df.samples[[i]] = subset(data_frame, Sample == blank[[i]][1])
df.blank[[i]] = subset(data_frame, Sample == blank[[i]][2])
df.samples[[i]]$Absorbance = df.samples[[i]]$Absorbance - df.blank[[i]]$Absorbance
}
k <- df.samples[[1]]
if (length(blank) > 1){
for (i in c(2:length(df.samples))){
k <- rbind(k, df.samples[[i]])
}
}
no.background = k
}else{
if (blank == "none"){ #routine for no blank in the data set
no.background = data_frame
}else{ #routine for a single blank in the data set
samples <- {}
for (i in c(1:length(unique(data_frame$Sample)))){
samples[[i]] <- subset(data_frame, Sample == unique(data_frame$Sample)[i])
df.blank = subset(data_frame, Sample == blank)
samples[[i]]$Absorbance = samples[[i]]$Absorbance - (df.blank$Absorbance*samples[[i]]$Pathlength/df.blank$Pathlength)
}
k <- samples[[1]]
for (i in c(2:length(samples))){
k <- rbind(k, samples[[i]])
}
no.background <- subset(k, Sample != blank)
}
}
raw.df = no.background
####Calculate extinction coefficients####
if (NucAcid[1] == "Concentration"){
}else{
if (NucAcid[1] == "Custom"){
extcoef = c(sum(as.numeric(NucAcid[2:length(NucAcid)])), as.numeric(NucAcid[2:length(NucAcid)]))
names(extcoef) <- c("Total", as.character(2:length(NucAcid)))
}else{
tryCatch({
extcoef = MeltR::calc.extcoeff(NucAcid, wavelength, T)
},
error = function(e){print("There is no nucleotide extinction coefficient at this wavelength for a nucleotide you specified in your sequence. You will need to provide your own custom extinction coefficient")})
}
}
####Remove values not in fitTs####
if (is.null(fitTs) == FALSE){
if (is.list(fitTs) == FALSE){
ranges <- rep(list(fitTs), length(unique(no.background$Sample)))
}else{
ranges <- fitTs
}
a <- {}
for (i in c(1:length(unique(no.background$Sample)))){
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
if(length(which(a[[i]]$Temperature > ranges[[i]][2])) > 0){
a[[i]] <- a[[i]][ -which(a[[i]]$Temperature > ranges[[i]][2]) ,]
}
if(length(which(a[[i]]$Temperature < ranges[[i]][1])) > 0){
a[[i]] <- a[[i]][ -which(a[[i]]$Temperature < ranges[[i]][1]) ,]
}
}
b <- a[[1]]
if (length(a) > 1){
for (i in c(2:length(unique(no.background$Sample)))){
b <- rbind(b, a[[i]])
}
}
no.background <- b
}
####Calculate Ct for each curve####
if (NucAcid[1] == "Concentration"){
ct <- as.numeric(NucAcid[2:length(NucAcid)])
samples = {}
for (i in c(1:length(unique(no.background$Sample)))){
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
samples[[i]]$Ct <- ct[i]
}
k <- samples[[1]]
if (length(samples) > 1){
for (i in c(2:length(samples))){
#print(i)
k <- rbind(k, samples[[i]])
}
}
no.background <- k
}else{
samples <- {}
if (length(extcoef) == 3){
ct <- c()
for (i in c(1:length(unique(no.background$Sample)))){
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
df.raw = subset(raw.df, Sample == unique(no.background$Sample)[i])
if (is.atomic(extcoef)){
ct[i] <- 2*(df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef[[1]]*df.raw$Pathlength[1]))
}else{
ct[i] <- 2*(df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef$Total*df.raw$Pathlength[1]))
}
samples[[i]]$Ct <- ct[i]
plot(samples[[i]]$Temperature, 2*samples[[i]]$Absorbance/(samples[[i]]$Pathlength*samples[[i]]$Ct))
abline(h = extcoef$Total)
}
}
if (length(extcoef) == 2){
ct <- c()
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
samples[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
df.raw = subset(data_frame, Sample == unique(no.background$Sample)[i])
if (is.atomic(extcoef)){
ct[i] <- (df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef[[1]]*df.raw$Pathlength[1]))
}else{
ct[i] <- (df.raw$Absorbance[which.min(abs(df.raw$Temperature - concT))]/(extcoef$Total*df.raw$Pathlength[1]))
}
samples[[i]]$Ct <- ct[i]
}
}
k <- samples[[1]]
if (length(samples) > 1){
for (i in c(2:length(samples))){
k <- rbind(k, samples[[i]])
}
}
no.background = k
}
####Remove outliers####
if (is.na(outliers)){
}else{
for (i in 1:length(outliers)){
no.background = subset(no.background, Sample != outliers[i])
}
}
####Calculate starting thermo parameters for nls####
first.derive <- {}
T0.5 <- c()
T0.75 <- c()
startH <- c()
for (i in c(1:length(unique(no.background$Sample)))){
first.derive[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
fit.Em = lm(Absorbance ~ poly(Temperature, 20, raw=TRUE),
data = first.derive[[i]])
string.dE.dT = gsub("NA", "0", gsub(", ", " + ", toString(paste(0:20, "*", coef(fit.Em), "*x^", -1:19, sep =""))))
string.dE.dT2 = gsub("NA", "0", gsub(", ", " + ", toString(paste(0:20, "*", c(0,0:19),"*", coef(fit.Em), "*x^", -2:18, sep =""))))
dE.dT = function(x){eval(parse(text = string.dE.dT))}
dE.dT2 = function(x){eval(parse(text = string.dE.dT2))}
T0.5[i] = seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)[which.max(dE.dT(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)))]
T0.75[i] = seq(min(first.derive[[i]]$Temperature) + 10, max(first.derive[[i]]$Temperature) - 10, length.out = 1000)[which.min(dE.dT2(seq(min(first.derive[[i]]$Temperature) + 10, max(first.derive[[i]]$Temperature) - 10, length.out = 1000)))]
if (Mmodel == "Monomolecular.2State"){
startH[i] <- -0.0032/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
if (Mmodel == "Heteroduplex.2State"){
startH[i] <- -0.007/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
if (Mmodel == "Homoduplex.2State"){
startH[i] <- -0.0044/((1/(273.15 + T0.5[i]))-(1/(273.15 + T0.75[i])))
}
#svg("~/Desktop/Melt_curve_derivative.svg")
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$Absorbance,
# xlab = Temperature ~ (degree ~ C), ylab = "Absorbance")
#lines(first.derive[[i]]$Temperature, predict(fit.Em), col = "red")
#lines(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000),
# 1.4 +10*dE.dT(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)), col = "red")
#lines(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000),
# 1.475+10*dE.dT2(seq(min(first.derive[[i]]$Temperature) + 2, max(first.derive[[i]]$Temperature) - 2, length.out = 1000)), col = "blue")
#abline(v = T0.75[i])
#abline(v = T0.5[i])
#dev.off()
first.derive[[i]]$dA.dT = dE.dT(first.derive[[i]]$Temperature) + dE.dT(first.derive[[i]]$Temperature)*((first.derive[[i]]$Absorbance - predict(fit.Em))/first.derive[[i]]$Absorbance)
first.derive[[i]]$dA.dT2 = dE.dT2(first.derive[[i]]$Temperature) + dE.dT2(first.derive[[i]]$Temperature)*((first.derive[[i]]$Absorbance - predict(fit.Em))/first.derive[[i]]$Absorbance)
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$dA.dT)
#plot(first.derive[[i]]$Temperature, first.derive[[i]]$dA.dT2)
}
####Calculate starting baseline values for nls####
a <-{}
uppbl_fit <- {}
lowbl_fit <- {}
bl.start.plot <- {}
for (i in c(1:length(unique(no.background$Sample)))){
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
upperbl <- data.frame(
"Absorbance" = a[[i]]$Absorbance[which(a[[i]]$Absorbance >= quantile(a[[i]]$Absorbance, probs = 0.75))],
"Temperature" = a[[i]]$Temperature[which(a[[i]]$Absorbance >= quantile(a[[i]]$Absorbance, probs = 0.75))]
)
lowbl <- data.frame(
"Absorbance" = a[[i]]$Absorbance[which(a[[i]]$Absorbance <= quantile(a[[i]]$Absorbance, probs = 0.25))],
"Temperature" = a[[i]]$Temperature[which(a[[i]]$Absorbance <= quantile(a[[i]]$Absorbance, probs = 0.25))]
)
uppbl_fit[[i]] <- lm(Absorbance ~ Temperature, data = upperbl)
lowbl_fit[[i]] <- lm(Absorbance ~ Temperature, data = lowbl)
}
####Method 1 fit each curve individually####
if (methods[1] == TRUE){
a <-{}
fit <- {}
indvfits.H <- c()
indvfits.S <- c()
indvfits.G <- c()
indvfits.Tm <- c()
indvfits.SE.Tm <- c()
bED <- c()
mED <- c()
bSS <- c()
mSS <- c()
Ind.model = Model
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
tryCatch({
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i])
if (is.list(auto.trimmed)){
fitstart = auto.trimmed[[i]]
}else{
fitstart <- list(H = startH[i], Tm = T0.5[i],
mED = lowbl_fit[[i]]$coefficients[2], bED = lowbl_fit[[i]]$coefficients[1],
mESS = uppbl_fit[[i]]$coefficients[2], bESS = uppbl_fit[[i]]$coefficients[1])
}
if (Mmodel == "Monomolecular.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Heteroduplex.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature, Ct),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Homoduplex.2State"){
fit[[i]] <- nls(Absorbance ~ Model(H, Tm, mED, bED, mESS, bESS, Temperature, Ct),
data = a[[i]],
start = fitstart,
trace = FALSE,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
mED[i] <- coef(fit[[i]])[3]
bED[i] <- coef(fit[[i]])[4]
mSS[i] <- coef(fit[[i]])[5]
bSS[i] <- coef(fit[[i]])[6]
indvfits.H[i] <- coef(fit[[i]])[1]
if (Mmodel == "Monomolecular.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2])
}
if (Mmodel == "Heteroduplex.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
}
if (Mmodel == "Homoduplex.2State"){
indvfits.S[i] <- calcS(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
indvfits.G[i] <- calcG(coef(fit[[i]])[1], coef(fit[[i]])[2], a[[i]]$Ct[1])
}
indvfits.Tm[i] <- coef(fit[[i]])[2]
indvfits.SE.Tm[i] <- coef(summary(fit[[i]]))[2,2]
}, error = function(e){print(paste("There was a problem fitting Sample", a[[i]]$Sample[1], "with method 1") )})
}
indvfits <- data.frame("Sample" = unique(no.background$Sample),
"Ct" = unique(no.background$Ct),
"H" = round(indvfits.H, 2),
"S" = round(indvfits.S, 5),
"G" = round(indvfits.G, 2),
"Tm" = round(indvfits.Tm, 2))
indvfits.mean <- list(round(mean(indvfits$H), 2), round(sd(indvfits$H), 2), round(1000*mean(indvfits$S), 2), round(1000*sd(indvfits$S), 2), round(mean(indvfits.G), 2), round(sd(indvfits.G), 2))
names(indvfits.mean) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
indvfits.mean <- data.frame(indvfits.mean)
if (Mmodel == "Heteroduplex.2State"){
no.background$Ext <- 2*no.background$Absorbance/(no.background$Pathlength*no.background$Ct)
}else{
no.background$Ext <- no.background$Absorbance/(no.background$Pathlength*no.background$Ct)
}
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_1_raw_fit_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(no.background$Temperature, no.background$Absorbance,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbance",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, predict(fit[[i]]), col = "red")
}
dev.off()
pdf(paste(file_path, "/", file_prefix, "_method_1_normalized_fit_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(no.background$Temperature, no.background$Ext,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbtivity (1/M*cm)",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
if (Mmodel == "Heteroduplex.2State"){
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, (2*predict(fit[[i]])/(a[[i]]$Ct[1]*a[[i]]$Pathlength[1])), col = "red")
}
}else{
for (i in c(1:length(a))){
lines(a[[i]]$Temperature, (predict(fit[[i]])/(a[[i]]$Ct[1]*a[[i]]$Pathlength[1])), col = "red")
}
}
dev.off()
}
}
if (methods[1] == FALSE){
indvfits <- data.frame("Sample" = NA,
"Ct" = NA,
"H" = NA,
"S" = NA,
"G" = NA,
"Tm" = NA)
indvfits.mean <- list(NA, NA, NA, NA, NA, NA)
names(indvfits.mean) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
indvfits.mean <- data.frame(indvfits.mean)
}
####Method 2 Tm vs Ct####
if (methods[2] == TRUE){
#Decide what Tm_method to use
if (Tm_method == "lm"){
a <- {}
Tm_range <- {}
Tm_fit <- {}
Tm <- c()
lnCt <- c()
for (i in c(1:length(unique(no.background$Sample)))){
#print(i)
a[[i]] <- subset(no.background, Sample == unique(no.background$Sample)[i]) #Pull out data
#calculate f at each temp
a[[i]]$f <- (a[[i]]$Absorbance - (coef(fit[[i]])[5]*a[[i]]$Temperature + coef(fit[[i]])[6]))/((coef(fit[[i]])[3]*a[[i]]$Temperature + coef(fit[[i]])[4]) - (coef(fit[[i]])[5]*a[[i]]$Temperature + coef(fit[[i]])[6]))
Tm_range[[i]] <- data.frame( #Pull out data where f >= 0.4 and f <= 0.6
"Temperature" = a[[i]]$Temperature,
"f" = a[[i]]$f
)
Tm_fit[[i]] <- lm(f~Temperature, data = Tm_range[[i]]) #fit data where f >= 0.4 and f <= 0.6
Tm[i] <- (0.5 - coef(Tm_fit[[i]])[1])/coef(Tm_fit[[i]])[2]
lnCt[i] <- log(a[[i]]$Ct[1])
}
Tm_data <- data.frame("lnCt" = lnCt, "Tm" = Tm)
Weight_Tm_M2 = FALSE
}
if (Tm_method == "nls"){
lnCt = log(indvfits$Ct)
Tm = indvfits$Tm
Tm_data = data.frame(lnCt, Tm)
}
if (Tm_method == "polynomial"){
lnCt = log(indvfits$Ct)
Tm = T0.5
Tm_data = data.frame(lnCt, Tm)
Weight_Tm_M2 = FALSE
}
Tm_data$invT <- 1/(Tm_data$Tm + 273.15)
if (Mmodel != "Monomolecular.2State"){
if (Mmodel != "Monomolecular.3State"){
Tm.model = TmModel
if (Weight_Tm_M2){ #Weighted non-linear regression
propagated.weights = indvfits.SE.Tm/((Tm_data$Tm+273.15)^2) #Propagate Tm SE
Tm_vs_lnCt_fit <- nls(invT ~ TmModel(H, S, lnCt),
data = Tm_data,
weights = 1/propagated.weights,
start = list(H = mean(indvfits$H), S = mean(indvfits$S)),
control = nls.control(warnOnly = TRUE))
}else{
Tm_vs_lnCt_fit <- nls(invT ~ TmModel(H, S, lnCt),
data = Tm_data,
start = list(H = mean(indvfits$H), S = mean(indvfits$S)),
control = nls.control(warnOnly = TRUE))
}
Tm_vs_lnCt <- list(round(coef(Tm_vs_lnCt_fit)[1], 2), round(summary(Tm_vs_lnCt_fit)$coefficients[1,2], 2),
round(1000*coef(Tm_vs_lnCt_fit)[2], 2), round(1000*summary(Tm_vs_lnCt_fit)$coefficients[2,2], 2),
round(coef(Tm_vs_lnCt_fit)[1] - (310.15*coef(Tm_vs_lnCt_fit)[2]), 2), round(sqrt((summary(Tm_vs_lnCt_fit)$coefficients[1,2])^2 + (310*summary(Tm_vs_lnCt_fit)$coefficients[2,2])^2 - 2*310*summary(Tm_vs_lnCt_fit)$cov.unscaled[1,2]*(summary(Tm_vs_lnCt_fit)$sigma^2)),2))
names(Tm_vs_lnCt) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Tm_vs_lnCt)
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_2_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
if (Weight_Tm_M2){
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8,
ylim = c(min(Tm_data$invT) - max(propagated.weights), max(Tm_data$invT) + max(propagated.weights)))
arrows(x = Tm_data$lnCt, x0 = Tm_data$lnCt, length = 0.05,
y = Tm_data$invT - propagated.weights, y0 = Tm_data$invT + propagated.weights,
code = 3, lwd=2, angle = 90)
lines(x = Tm_data$lnCt, predict(Tm_vs_lnCt_fit), col = "red")
}else{
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
lines(x = Tm_data$lnCt, predict(Tm_vs_lnCt_fit), col = "red")
}
dev.off()
}
}else{
Tm.model = NA
}
}else{
Tm.model = NA
}
if (Mmodel == "Monomolecular.2State"){
if (Mmodel == "Monomolecular.3State"){
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_2_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(x = Tm_data$lnCt, y = Tm_data$invT,
xlab = "ln[ Ct (M) ]", ylab = "1/[ Temperature (K) ]",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
lines(Tm_data$lnCt, predict(lm(Tm_data$invT ~ Tm_data$lnCt)), col = "red")
dev.off()
}
}}
}
if (methods[2] == FALSE){
Tm_vs_lnCt <- list(NA, NA,
NA, NA,
NA, NA)
names(Tm_vs_lnCt) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Tm_vs_lnCt)
}
####Method 3 Global fitting####
if (methods[3] == TRUE){
b <- data.frame("Sample" = c(), "Pathlength" = c(), "Temperature" = c(),
"Absorbance" = c(), "Ct" = c(), "Ext" = c())
for (i in c(1:length(unique(no.background$Sample)))){
a <- subset(no.background, Sample == unique(no.background$Sample)[i])
a$Sample <- i
b <- rbind(b, a)
}
gfit_data <- b
gfit_start = list(H = mean(indvfits$H), S = mean(indvfits$S), mED = mED, bED = bED, mESS = mSS, bESS = bSS)
if (Mmodel == "Monomolecular.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Heteroduplex.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Mmodel == "Homoduplex.2State"){
gfit <- nls(Absorbance ~ GModel(H, S, mED, bED, mESS, bESS, Sample, Temperature, Ct),
start = gfit_start,
data = gfit_data,
nls.control(tol = 5e-04, minFactor = 1e-10, maxiter = 50, warnOnly = TRUE))
}
if (Save_results == "all"){
pdf(paste(file_path, "/", file_prefix, "_method_3_Gfit_raw_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(gfit_data$Temperature, gfit_data$Absorbance,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbance",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(unique(gfit_data$Sample)))){
a <- subset(gfit_data, Sample == unique(gfit_data$Sample)[i])
if (Mmodel == "Monomolecular.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature)))),
col = "red")
}
if (Mmodel == "Heteroduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1]),
col = "red")
}
if (Mmodel == "Homoduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1]),
col = "red")
}
}
dev.off()
pdf(paste(file_path, "/", file_prefix, "_method_3_Gfit_normalized_plot.pdf", sep = ""),
width = 3, height = 3, pointsize = 0.25)
plot(gfit_data$Temperature, gfit_data$Ext,
xlab = Temperature ~ (degree ~ C), ylab = "Absorbtivity (1/M*cm)",
cex.lab = 1.5, cex.axis = 1.25, cex = 0.8)
for (i in c(1:length(unique(gfit_data$Sample)))){
a <- subset(gfit_data, Sample == unique(gfit_data$Sample)[i])
if (Mmodel == "Monomolecular.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))))/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
if (Mmodel == "Heteroduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), 2*GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1])/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
if (Mmodel == "Homoduplex.2State"){
lines(c(floor(min(a$Temperature)):ceiling(max(a$Temperature))), GModel(H = coef(gfit)[1],
S = coef(gfit)[2],
mED = coef(gfit)[i + 2],
bED = coef(gfit)[i + 2 + length(unique(gfit_data$Sample))],
mESS = coef(gfit)[i + 2 + 2*length(unique(gfit_data$Sample))],
bESS = coef(gfit)[i + 2 + 3*length(unique(gfit_data$Sample))],
Temperature = c(floor(min(a$Temperature)):ceiling(max(a$Temperature))),
Ct = a$Ct[1])/(a$Pathlength[1]*a$Ct[1]),
col = "red")
}
}
dev.off()
}
Gfit_summary <- list(round(coef(gfit)[1], 2), round(summary(gfit)$coefficients[1,2], 2),
round(1000*coef(gfit)[2], 2), round(1000*summary(gfit)$coefficients[2,2], 2),
round(coef(gfit)[1] - (310.15)*coef(gfit)[2], 2), round(sqrt((summary(gfit)$coefficients[1,2])^2 + (310.15*summary(gfit)$coefficients[2,2])^2 - 2*310.15*summary(gfit)$cov.unscaled[1,2]*(summary(gfit)$sigma^2)), 2))
names(Gfit_summary) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Gfit_summary <- data.frame(Gfit_summary)
}
if (methods[3] == FALSE){
Gfit_summary <- list(NA, NA,
NA, NA,
NA, NA)
names(Gfit_summary) <- c("H", "SE.H", "S", "SE.S", "G", "SE.G")
Tm_vs_lnCt <- data.frame(Gfit_summary)
}
####Assemble Results####
indvfits$Ct = as.numeric(formatC(indvfits$Ct, format = "e", digits = 2))
indvfits$S = 1000*indvfits$S
colnames(indvfits) = c("Sample", "Ct", "dH", "dS", "dG", "Tm")
if (Mmodel == "Monomolecular.2State"){
comparison <- rbind(indvfits.mean, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "3 Global fit")), comparison)
row.names(comparison) <- c(1:2)
}
if (Mmodel == "Heteroduplex.2State"){
comparison <- rbind(indvfits.mean, Tm_vs_lnCt, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "2 Tm versus ln[Ct]", "3 Global fit")), comparison)
row.names(comparison) <- c(1:3)
}
if (Mmodel == "Homoduplex.2State"){
comparison <- rbind(indvfits.mean, Tm_vs_lnCt, Gfit_summary)
comparison <- cbind(data.frame("Method" = c("1 Individual fits", "2 Tm versus ln[Ct]", "3 Global fit")), comparison)
row.names(comparison) <- c(1:3)
}
####Expectation maximized Tm calculation####
Tm_at_0.1mM = c()
SE.Tm_at_0.1mM = c()
for (i in 1:nrow(comparison)){
if (i == 1){
df.covar = indvfits
udH = df.covar$dH - mean(df.covar$dH)
udS = (df.covar$dS - mean(df.covar$dS))/1000
sigma = sum(udH*udS)/nrow(df.covar)
dH = mean(df.covar$dH)
dS = mean(df.covar$dS)/1000
SE.dH = sd(df.covar$dH)
SE.dS = sd(df.covar$dS)/1000
}
if (i == 2){
if (Mmodel == "Monomolecular.2State"){
SE.dH = NA
SE.dS = NA
}else{
if (methods[2] == F){
sigma = NA
dH = NA
dS = NA
SE.dH = NA
SE.dS = NA
}else{
sigma = summary(Tm_vs_lnCt_fit)$cov.unscaled[1,2]*(summary(Tm_vs_lnCt_fit)$sigma^2)
dH = coef(Tm_vs_lnCt_fit)[1]
dS = coef(Tm_vs_lnCt_fit)[2]
SE.dH = coef(summary(Tm_vs_lnCt_fit))[1,2]
SE.dS = coef(summary(Tm_vs_lnCt_fit))[2,2]
}
}
}
if (i == 3){
if (methods[2] == F){
sigma = NA
dH = NA
dS = NA
SE.dH = NA
SE.dS = NA
}else{
sigma = summary(gfit)$cov.unscaled[1,2]*(summary(gfit)$sigma^2)
dH = coef(gfit)[1]
dS = coef(gfit)[2]
SE.dH = coef(summary(gfit))[1,2]
SE.dS = coef(summary(gfit))[2,2]
}
}
Tm_at_0.1mM[i] = round(calcTM(H = dH, S = dS,
Ct = 10^-4), digits = 2)
SE.Tm_at_0.1mM[i] = round(calcTM.SE(H = dH, S = dS,
SE.H = SE.dH, SE.S = SE.dS,
Ct = 10^-4, covar = sigma), digits = 2)
}
comparison$Tm_at_0.1mM = Tm_at_0.1mM
comparison$SE.Tm_at_0.1mM = SE.Tm_at_0.1mM
colnames(comparison) = c("Method", "dH", "SE.dH", "dS", "SE.dS", "dG", "SE.dG", "Tm_at_0.1mM", "SE.Tm_at_0.1mM")
if (Save_results != "none"){
write.table(comparison, paste(file_path, "/", file_prefix, "_summary.csv", sep = ""), sep = ",", row.names = FALSE)
}
if (Save_results != "none"){
write.table(indvfits, paste(file_path, "/", file_prefix, "_method_1_individual_fits.csv", sep = ""), sep = ",", row.names = FALSE)
}
####Assemble data for custom derivative plots####
df.deriv = first.derive[[1]]
if (length(first.derive) > 1){
for (i in 2:length(first.derive)){
#print(i)
df.deriv = rbind(df.deriv, first.derive[[i]])
}
}
####Assemble data for custom individual fit plots####
df.method.1 = subset(no.background, Sample == unique(no.background$Sample)[1])
df.method.1$Model = predict(fit[[1]])
if (length(unique(no.background$Sample)) >1){
for (i in 2:length(unique(no.background$Sample))){
df.method.1.i = subset(no.background, Sample == unique(no.background$Sample)[i])
df.method.1.i$Model = predict(fit[[i]])
df.method.1 = rbind(df.method.1, df.method.1.i)
}
}
####Assemble data for custom 1/Tm versus lnCt plots####
if (Mmodel != "Monomolecular.2State"){
if (methods[2] == TRUE){
Tm_data$Model = predict(Tm_vs_lnCt_fit)
}
}
####Assemble global fit data for custom fit plots####
if (methods[3] == TRUE){
gfit_data$Model = predict(gfit)
}
####Compile data for the BLTrimmer####
BLTrimmer = list(raw.df, #background subtracted raw data
Mmodel,
Tmodel)
####Make a list of MeltR A settings####
settings = list(data_frame, #1
blank, #2
NucAcid, #3
concT, #4
fitTs, #5
methods, #6
Mmodel, #7
Tmodel, #8
Save_results, #9
file_prefix, #10
file_path, #11
Weight_Tm_M2) #12
####Assemble final output####
if (Silent){}else{
print("Individual curves")
print("dH and dG are in kcal/mol and dS is in cal/mol/K. Tms are in deg Celsius")
print(indvfits[order(indvfits$Ct),])
print("Summary")
print("dH and dG are in kcal/mol and dS is in cal/mol/K. Tms are in deg Celsius")
print(comparison)
range <- data.frame("dH" = round(100*abs((range(comparison$dH)[1]-range(comparison$dH)[2])/mean(comparison$dH)), digits = 1),
"dS" = round(100*abs((range(comparison$dS)[1]-range(comparison$dS)[2])/mean(comparison$dS)), digits = 1),
"dG" = round(100*abs((range(comparison$dG)[1]-range(comparison$dG)[2])/mean(comparison$dG)), digits = 1))
print("Maximum %error across methods")
print(range)
}
output <- list("Summary" = comparison,
"Method.1.indvfits" = indvfits,
"Range" = range,
"Derivatives.data" = df.deriv,
"Method.1.data" = df.method.1,
"Method.1.fit" = fit,
"BLTrimmer.data" = BLTrimmer,
"meltR.A.settings" = settings)
if (Mmodel != "Monomolecular.2State"){
if (methods[2] == TRUE){
output$Method.2.data = Tm_data
output$Method.2.fit <- Tm_vs_lnCt_fit
}
}
if (methods[3] == TRUE){
output$Method.3.data = gfit_data
output$Method.3.fit <- gfit
}
output <- output
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.