Advanced Usage

# Custom objective: maximize total variance
variance_objective <- function(data) {
  sum(apply(data, 2, var))
}

# Test the function
test_data <- reference[, 1:50]
cat("Variance objective value:", variance_objective(test_data), "\n")
#> Variance objective value: 200.312

# Custom objective: maximize marker specificity
# Higher when genes are specific to individual cell types
marker_score <- function(data) {
  # Max expression / mean expression ratio
  col_max <- apply(data, 2, max)
  col_mean <- colMeans(data)
  col_mean[col_mean == 0] <- 1e-10
  sum(col_max / col_mean)
}

cat("Marker score value:", marker_score(test_data), "\n")
#> Marker score value: 107.8375

dw_custom <- darwin(reference)
dw_custom$optimize(
  ngen = 60,
  objectives = c("correlation", variance_objective),
  weights = c(-1, 1),  # Minimize correlation, maximize variance
  verbose = FALSE,
  parallel = FALSE
)

# The second objective column will show variance values
fitness <- dw_custom$get_fitness()
colnames(fitness) <- c("Correlation", "Variance")

ggplot(as.data.frame(fitness), aes(x = Correlation, y = Variance)) +
  geom_point(color = "#3498db", size = 3, alpha = 0.7) +
  geom_line(color = "gray60", alpha = 0.5) +
  labs(
    title = "Custom Objectives: Correlation vs Variance",
    subtitle = "Minimize correlation, maximize variance",
    x = "Correlation (lower is better)",
    y = "Variance (higher is better)"
  ) +
  theme_minimal(base_size = 12)

dw_three <- darwin(reference)
dw_three$optimize(
  ngen = 60,
  objectives = c("correlation", "distance", "condition"),
  weights = c(-1, 1, -1),  # Minimize corr/cond, maximize distance
  verbose = FALSE,
  parallel = FALSE
)

fitness3 <- dw_three$get_fitness()
cat("Three-objective optimization:\n")
#> Three-objective optimization:
cat("  Solutions:", nrow(fitness3), "\n")
#>   Solutions: 314
cat("  Correlation range:", round(range(fitness3[,1]), 2), "\n")
#>   Correlation range: 0.21 0.79
cat("  Distance range:", round(range(fitness3[,2]), 2), "\n")
#>   Distance range: 143.43 363.86
cat("  Condition range:", round(range(fitness3[,3]), 2), "\n")
#>   Condition range: 3.93 4.27

dw <- darwin(reference)
dw$optimize(ngen = 50, verbose = FALSE, parallel = FALSE)

# Emphasize minimizing correlation
dw$select(weights = c(-2, 1))
cat("Emphasize correlation - genes:", sum(dw$get_selection()), "\n")
#> Emphasize correlation - genes: 388

# Emphasize maximizing distance
dw$select(weights = c(-1, 2))
cat("Emphasize distance - genes:", sum(dw$get_selection()), "\n")
#> Emphasize distance - genes: 388

# Select by direct index
dw$select(index = 1)
cat("Solution 1 - genes:", sum(dw$get_selection()), "\n")
#> Solution 1 - genes: 400

# Select by objective rank
# (objective_index, rank) - rank 1 = best for that objective
dw$select(index = c(1, 1))  # Best correlation
cat("Best correlation - genes:", sum(dw$get_selection()), "\n")
#> Best correlation - genes: 326

dw$select(index = c(2, -1))  # Best distance (last rank)
cat("Best distance - genes:", sum(dw$get_selection()), "\n")
#> Best distance - genes: 400

# Select solution closest to target correlation value
fitness <- dw$get_fitness()
target_corr <- median(fitness$correlation)
dw$select(close_to = c(1, target_corr))
cat("Target correlation", round(target_corr, 2), "- genes:", sum(dw$get_selection()), "\n")
#> Target correlation 0.42 - genes: 397

# Benchmark (example, not run)
library(microbenchmark)

# Large dataset
large_ref <- matrix(rnorm(50 * 2000), nrow = 50)

microbenchmark(
  serial = {
    dw <- darwin(large_ref)
    dw$optimize(ngen = 10, parallel = FALSE, verbose = FALSE)
  },
  parallel = {
    dw <- darwin(large_ref)
    dw$optimize(ngen = 10, parallel = TRUE, verbose = FALSE)
  },
  times = 3
)

results <- list()
for (pop in c(30, 60, 100)) {
  dw_temp <- darwin(reference)
  dw_temp$optimize(
    ngen = 40,
    pop_size = pop,
    verbose = FALSE,
    parallel = FALSE
  )
  results[[as.character(pop)]] <- data.frame(
    dw_temp$get_fitness(),
    pop_size = factor(pop)
  )
}

df_pop <- do.call(rbind, results)

ggplot(df_pop, aes(x = correlation, y = distance, color = pop_size)) +
  geom_point(alpha = 0.6) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Effect of Population Size",
    x = "Correlation",
    y = "Distance",
    color = "Population\nSize"
  ) +
  theme_minimal(base_size = 12)

results_mut <- list()
for (mut in c(0.001, 0.01, 0.05)) {
  dw_temp <- darwin(reference)
  dw_temp$optimize(
    ngen = 40,
    mutation_prob = mut,
    verbose = FALSE,
    parallel = FALSE
  )
  
  pareto <- dw_temp$get_pareto()
  gene_counts <- sapply(pareto, sum)
  
  results_mut[[as.character(mut)]] <- data.frame(
    n_genes = gene_counts,
    mutation = factor(mut)
  )
}

df_mut <- do.call(rbind, results_mut)

ggplot(df_mut, aes(x = n_genes, fill = mutation)) +
  geom_density(alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  labs(
    title = "Effect of Mutation Probability on Gene Count",
    x = "Number of Selected Genes",
    y = "Density",
    fill = "Mutation\nProbability"
  ) +
  theme_minimal(base_size = 12)

library(Seurat)

# From Seurat object with annotations
dw <- darwin(
  seurat_obj,
  celltype_key = "cell_type",
  assay = "RNA",
  layer = "data",
  use_highly_variable = TRUE
)

# Optimize and select
dw$optimize(ngen = 100)
dw$select()

# Get marker genes
markers <- dw$get_genes()

# Get selected genes for use with other tools
genes <- dw$get_genes()
selection <- dw$get_selection()

# Create signature matrix for CIBERSORT
signature_matrix <- reference[, selection]
write.csv(signature_matrix, "signature_matrix.csv")

# Or as gene list for MuSiC
gene_list <- genes

# Check optimization quality
fitness <- dw$get_fitness()

cat("Diagnostic Summary:\n")
#> Diagnostic Summary:
cat("  Pareto front size:", nrow(fitness), "\n")
#>   Pareto front size: 99
cat("  Correlation range:", round(diff(range(fitness[,1])), 3), "\n")
#>   Correlation range: 0.149
cat("  Distance range:", round(diff(range(fitness[,2])), 3), "\n")
#>   Distance range: 47.774
cat("  Gene count range:", range(sapply(dw$get_pareto(), sum)), "\n")
#>   Gene count range: 326 400