# Setup: Load essential R packages for chemical analysis

# Avoid interactive CRAN prompts
options(repos = c(CRAN = "https://cloud.r-project.org"))

required_packages <- c("ggplot2", "dplyr", "plotly", "corrplot", "reshape2", "RColorBrewer")

# Install missing packages quietly
missing <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing)) install.packages(missing, quiet = TRUE)

# Attach packages without printing list results or masking chatter
suppressPackageStartupMessages({
  for (pkg in required_packages) {
    suppressWarnings(library(pkg, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))
  }
})
# Alternatively:
# invisible(lapply(required_packages, function(pkg)
#   suppressWarnings(library(pkg, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))
# ))

# Theme fix: use linewidth instead of size to avoid ggplot2 deprecation warning
chemistry_theme <- ggplot2::theme_minimal() +
  ggplot2::theme(
    plot.title = ggplot2::element_text(size = 16, face = "bold", color = "#2E86AB"),
    plot.subtitle = ggplot2::element_text(size = 12, color = "#A23B72"),
    axis.title = ggplot2::element_text(size = 12, face = "bold"),
    axis.text = ggplot2::element_text(size = 10),
    legend.title = ggplot2::element_text(size = 11, face = "bold"),
    panel.grid.minor = ggplot2::element_blank(),
    panel.border = ggplot2::element_rect(color = "gray80", fill = NA, linewidth = 0.5)
  )

cat("Chemical Analysis Toolkit Loaded Successfully!\n")
cat("Ready to explore the molecular world with R\n")

# Molecular descriptor analysis and QSAR modeling
# Create synthetic dataset representing typical pharmaceutical compounds
set.seed(42)  # For reproducibility

# Generate molecular descriptors for drug-like compounds
n_compounds <- 100
molecular_descriptors <- data.frame(
  compound_id = paste0("COMP_", sprintf("%03d", 1:n_compounds)),
  molecular_weight = runif(n_compounds, 150, 800),
  logP = runif(n_compounds, -2, 6),  # Lipophilicity
  h_bond_donors = rpois(n_compounds, 2),
  h_bond_acceptors = rpois(n_compounds, 4),
  rotatable_bonds = rpois(n_compounds, 5),
  polar_surface_area = runif(n_compounds, 20, 200),
  aromatic_rings = rpois(n_compounds, 2)
)

# Calculate synthetic bioactivity based on realistic relationships
molecular_descriptors$bioactivity <- with(molecular_descriptors, {
  # Lipinski's Rule of Five influences
  lipinski_penalty <- ifelse(molecular_weight > 500, -0.5, 0) +
                     ifelse(logP > 5, -0.3, 0) +
                     ifelse(h_bond_donors > 5, -0.2, 0) +
                     ifelse(h_bond_acceptors > 10, -0.2, 0)
  
  # Optimal activity model
  base_activity <- 2 + 
                   0.3 * aromatic_rings + 
                   0.1 * h_bond_acceptors - 
                   0.05 * rotatable_bonds +
                   0.2 * logP - 0.02 * logP^2 +  # Optimal logP around 5
                   lipinski_penalty +
                   rnorm(n_compounds, 0, 0.3)  # Random variation
  
  pmax(0, base_activity)  # Ensure non-negative
})

# Create Lipinski's Rule of Five compliance
molecular_descriptors$lipinski_compliant <- with(molecular_descriptors, 
  (molecular_weight <= 500) & (logP <= 5) & 
  (h_bond_donors <= 5) & (h_bond_acceptors <= 10)
)

# Build predictive QSAR model
qsar_model <- lm(bioactivity ~ molecular_weight + logP + I(logP^2) + 
                 h_bond_donors + h_bond_acceptors + rotatable_bonds + 
                 polar_surface_area + aromatic_rings, 
                 data = molecular_descriptors)

# Add predictions to dataset
molecular_descriptors$predicted_activity <- predict(qsar_model)
molecular_descriptors$residuals <- residuals(qsar_model)

# Create comprehensive QSAR visualization
p_qsar <- ggplot(molecular_descriptors, aes(x = predicted_activity, y = bioactivity)) +
  geom_point(aes(color = lipinski_compliant, size = molecular_weight), alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "red", alpha = 0.3) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("FALSE" = "orange", "TRUE" = "green"), 
                     name = "Lipinski\nCompliant") +
  scale_size_continuous(range = c(2, 6), name = "Molecular\nWeight") +
  labs(
    title = "QSAR Model: Predicted vs Observed Bioactivity",
    subtitle = "Machine learning prediction of drug-like properties from molecular descriptors",
    x = "Predicted Bioactivity",
    y = "Observed Bioactivity",
    caption = "Lipinski's Rule of Five compliance affects drug-like properties"
  ) +
  chemistry_theme

print(p_qsar)

# Model performance metrics
r_squared <- summary(qsar_model)$r.squared
rmse <- sqrt(mean(molecular_descriptors$residuals^2))
mae <- mean(abs(molecular_descriptors$residuals))

cat("\n QSAR Model Performance:\n")
cat("==========================\n")
cat(sprintf("R² = %.3f\n", r_squared))
cat(sprintf("RMSE = %.3f\n", rmse))
cat(sprintf("MAE = %.3f\n", mae))

# Lipinski compliance analysis
lipinski_analysis <- molecular_descriptors %>%
  group_by(lipinski_compliant) %>%
  summarise(
    count = n(),
    avg_bioactivity = mean(bioactivity),
    avg_mw = mean(molecular_weight),
    avg_logP = mean(logP),
    .groups = 'drop'
  )

cat("\n Lipinski Rule Analysis:\n")
cat("==========================\n")
print(lipinski_analysis)

cat("\n Key QSAR Insights:\n")
cat("• Lipinski-compliant compounds show better drug-like properties\n")
cat("• Optimal logP around 2-3 for bioactivity\n")
cat("• Molecular weight <500 Da preferred for oral drugs\n")
cat("• R integration enables rapid cheminformatics analysis!\n")