Skip to content

Commit

Permalink
Correct the indentation of the R script using styler
Browse files Browse the repository at this point in the history
  • Loading branch information
@rmassei
rmassei committed Dec 16, 2024
1 parent 605b980 commit 026651d
Showing 1 changed file with 83 additions and 83 deletions.
166 changes: 83 additions & 83 deletions tools/tox_tools/dose_responses/dose_response.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,104 +2,104 @@ library(drc)
library(ggplot2)

fit_models <- function(data, concentration_col, response_col) {
models <- list(
LL.2 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = LL.2(), type = "binomial"),
W1.4 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = W1.4(), type = "binomial"),
W2.4 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = W2.4(), type = "binomial"),
BC.5 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = BC.5(), type = "binomial")
)
return(models)
models <- list(
LL.2 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = LL.2(), type = "binomial"),
W1.4 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = W1.4(), type = "binomial"),
W2.4 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = W2.4(), type = "binomial"),
BC.5 = drm(data[[response_col]] ~ data[[concentration_col]], data = data, fct = BC.5(), type = "binomial")
)
return(models)
}

select_best_model <- function(models) {
aic_values <- sapply(models, AIC)
best_model_name <- names(which.min(aic_values))
best_model <- models[[best_model_name]]
return(list(name = best_model_name, model = best_model))
aic_values <- sapply(models, AIC)
best_model_name <- names(which.min(aic_values))
best_model <- models[[best_model_name]]
return(list(name = best_model_name, model = best_model))
}

calculate_ec_values <- function(model) {
ec50 <- ED(model, 50, type = "relative")
ec25 <- ED(model, 25, type = "relative")
ec10 <- ED(model, 10, type = "relative")
return(list(EC50 = ec50, EC25 = ec25, EC10 = ec10))
ec50 <- ED(model, 50, type = "relative")
ec25 <- ED(model, 25, type = "relative")
ec10 <- ED(model, 10, type = "relative")
return(list(EC50 = ec50, EC25 = ec25, EC10 = ec10))
}

plot_dose_response <- function(model, data, ec_values, concentration_col, response_col, plot_file, compound_name, concentration_unit) {
# Generate a grid of concentration values for predictions
concentration_grid <- seq(min(data[[concentration_col]]), max(data[[concentration_col]]), length.out = 100)
prediction_data <- data.frame(concentration_grid)
colnames(prediction_data) <- concentration_col

# Compute predictions with confidence intervals
predictions <- predict(model, newdata = prediction_data, type = "response", interval = "confidence")
prediction_data$resp <- predictions[, 1]
prediction_data$lower <- predictions[, 2]
prediction_data$upper <- predictions[, 3]

print(prediction_data)

data$rep <- factor(data$rep)

# Create the plot
p <- ggplot(data, aes_string(x = concentration_col, y = response_col)) +
geom_point(aes(colour = rep)) + # Original data points
geom_line(data = prediction_data, aes_string(x = "conc", y = "resp"), color = "blue") + # Predicted curve
geom_ribbon(data = prediction_data, aes_string(x = "conc", ymin = "lower", ymax = "upper"), alpha = 0.2, fill = "blue") + # Confidence intervals
geom_vline(xintercept = ec_values$EC10[1], color = "green", linetype = "dashed") +
geom_vline(xintercept = ec_values$EC50[1], color = "red", linetype = "dashed") +
labs(
title = paste(compound_name, "- Dose-Response Curve"),
x = paste("Dose [", concentration_unit, "]"),
y = "Response %"
) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", color = NA),
plot.background = element_rect(fill = "white", color = NA)
)

# Save the plot to a file
jpeg(filename = plot_file)
print(p)
dev.off()
# Generate a grid of concentration values for predictions
concentration_grid <- seq(min(data[[concentration_col]]), max(data[[concentration_col]]), length.out = 100)
prediction_data <- data.frame(concentration_grid)
colnames(prediction_data) <- concentration_col

# Compute predictions with confidence intervals
predictions <- predict(model, newdata = prediction_data, type = "response", interval = "confidence")
prediction_data$resp <- predictions[, 1]
prediction_data$lower <- predictions[, 2]
prediction_data$upper <- predictions[, 3]

print(prediction_data)

data$rep <- factor(data$rep)

# Create the plot
p <- ggplot(data, aes_string(x = concentration_col, y = response_col)) +
geom_point(aes(colour = rep)) + # Original data points
geom_line(data = prediction_data, aes_string(x = "conc", y = "resp"), color = "blue") + # Predicted curve
geom_ribbon(data = prediction_data, aes_string(x = "conc", ymin = "lower", ymax = "upper"), alpha = 0.2, fill = "blue") + # Confidence intervals
geom_vline(xintercept = ec_values$EC10[1], color = "green", linetype = "dashed") +
geom_vline(xintercept = ec_values$EC50[1], color = "red", linetype = "dashed") +
labs(
title = paste(compound_name, "- Dose-Response Curve"),
x = paste("Dose [", concentration_unit, "]"),
y = "Response %"
) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", color = NA),
plot.background = element_rect(fill = "white", color = NA)
)

# Save the plot to a file
jpeg(filename = plot_file)
print(p)
dev.off()
}

dose_response_analysis <- function(data, concentration_col, response_col, plot_file, ec_file, compound_name, concentration_unit) {
# Ensure column names are correctly selected
concentration_col <- colnames(data)[as.integer(concentration_col)]
response_col <- colnames(data)[as.integer(response_col)]

# Fit models and select the best one
models <- fit_models(data, concentration_col, response_col)
best_model_info <- select_best_model(models)
best_model <- best_model_info$model
best_model_name <- best_model_info$name

# Calculate EC values
ec_values <- calculate_ec_values(best_model)

# Plot the dose-response curve
plot_dose_response(best_model, data, ec_values, concentration_col, response_col, plot_file, compound_name, concentration_unit)

# Get model summary and AIC value
model_summary <- summary(best_model)
model_aic <- AIC(best_model)

# Prepare EC values data frame with additional information
ec_data <- data.frame(
Metric = c("chemical_name", "EC10", "EC25", "EC50", "AIC"),
Value = c(compound_name, ec_values$EC10[1], ec_values$EC25[1], ec_values$EC50[1], model_aic)
)
# Ensure column names are correctly selected
concentration_col <- colnames(data)[as.integer(concentration_col)]
response_col <- colnames(data)[as.integer(response_col)]

# Fit models and select the best one
models <- fit_models(data, concentration_col, response_col)
best_model_info <- select_best_model(models)
best_model <- best_model_info$model
best_model_name <- best_model_info$name

# Calculate EC values
ec_values <- calculate_ec_values(best_model)

# Plot the dose-response curve
plot_dose_response(best_model, data, ec_values, concentration_col, response_col, plot_file, compound_name, concentration_unit)

# Get model summary and AIC value
model_summary <- summary(best_model)
model_aic <- AIC(best_model)

# Prepare EC values data frame with additional information
ec_data <- data.frame(
Metric = c("chemical_name", "EC10", "EC25", "EC50", "AIC"),
Value = c(compound_name, ec_values$EC10[1], ec_values$EC25[1], ec_values$EC50[1], model_aic)
)

# Write EC values, AIC, and model summary to the output file
write.table(ec_data, ec_file, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)
# Write EC values, AIC, and model summary to the output file
write.table(ec_data, ec_file, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)

# Append the model summary to the file
cat("\nModel Summary:\n", file = ec_file, append = TRUE)
capture.output(model_summary, file = ec_file, append = TRUE)
# Append the model summary to the file
cat("\nModel Summary:\n", file = ec_file, append = TRUE)
capture.output(model_summary, file = ec_file, append = TRUE)

return(list(best_model = best_model_name, ec_values = ec_values))
return(list(best_model = best_model_name, ec_values = ec_values))
}

args <- commandArgs(trailingOnly = TRUE)
Expand Down

0 comments on commit 026651d

Please sign in to comment.