03_make_data.R

# 00 - setup ---------------------------------------------------------------

## set workdir to repo
setwd("D:\\git-repos\\data-repos\\nymphalid-phenomics")

## "pacman" does automatic package loading
if (!require("pacman")) install.packages("pacman")

pkgs <- c(
  # Data manipulation
  "data.table", "lubridate", "stringi", "stringr",
  "googlesheets4","fs",
  
  # Visualization
  "ggplot2", "cowplot", "ggridges", "ggnewscale", "pals",
  "plotly", "ggimage","gridExtra", "egg", "corrplot",
  "scales", "ks", "ggtree",
  
  # Statistical modeling and multivariate analysis
  "mgcv", "gratia", "FactoMineR", "factoextra",
  "AICcmodavg", "caret", "sda", "fields",
  "umap", "Rtsne", "mclust", "MLmetrics",
  
  # load tree
  "phytools",
  
  ## for utils
  "terra","MASS"
  
)

## load custom tools
source("scripts\\utils\\utils_r.R")

pacman::p_load(char = pkgs, install = F)
rm(pkgs)

theme_set(theme_cowplot())
theme_update(
  panel.grid.major = element_line(colour = "grey90"),
  panel.border = element_rect(colour = "black", fill=NA, linewidth=1))

# 01 - load data and reformat ---------------------------------------------------------------

## load files
data_emb = fread("data/data_primary/embeddings.csv", stringsAsFactors = TRUE, na.strings = "")
data_emb_moths = fread("data/data_primary/embeddings_moths.csv", stringsAsFactors = TRUE, na.strings = "")
data_feat = fread("data/data_primary/features.csv", stringsAsFactors = TRUE, na.strings = "")
data_meta = fread("data/data_primary/meta.csv", stringsAsFactors = TRUE, na.strings = "")
data_meta_moths = fread("data/data_primary/meta_moths.csv", stringsAsFactors = TRUE, na.strings = "")
data_labels = fread("data/data_primary/labels.csv", stringsAsFactors = TRUE, na.strings = "")
feature_key = fread("data/data_primary/feature_key.csv", stringsAsFactors = TRUE, na.strings = "")
data_tree_chaz2021_all = read.tree("data/data_primary/tree.tree")

## align embeddings and meta-data
data_meta = data_meta[order(subfamily, species, uuid, class_dv)]
data_emb = data_emb[match(data_meta$mask_name, data_emb$mask_name)]
data_meta_moths = data_meta_moths[order(species, mask_name)]
data_emb_moths = data_emb_moths[match(data_meta_moths$mask_name, data_emb_moths$mask_name)]

# 02 - variables  --------------------------------------------------------------

## get colnames
cols_feat = names(data_feat)[substr(names(data_feat),1,4) == "feat"]
cols_emb = names(data_emb)[substr(names(data_emb),1,3) == "emb"]

## colors
colors_unpala = c("unknown"="gray80", "low"="#1b9e77","medium"="#7570b3","high"="#d95f02")
labels_unpala = c("unknown"="NA", "low"="Low","medium"="Medium","high"="High")
colors_moths = c("high"="firebrick1", "low"="darkolivegreen")
labels_moths = c("high"="Aposematism", "low"="Camouflage")

# fam_cols =  c("black", glasbey(15)[5:15])
fam_cols = c(`0` = "black", heliconiinae = "#FF00B6", satyrinae = "#005300", 
             charaxinae = "#FFD300", nymphalinae = "#009FFF", danainae = "#9A4D42", 
             pseudergolinae = "#00FFBE", apaturinae = "#783FC1", biblidinae = "#1F9698", 
             limenitidinae = "#FFACFD", cyrestinae = "#B1CC71", libytheinae = "#F1085C"
)
fam_labs = c("Root", tools::toTitleCase(names(fam_cols)[2:12])) 
names(fam_labs) = names(fam_cols)

# > labels - phylo-imp ---------------------------------------------------

## check
ev_lev = c("sp", "gn")
data_labels[,table(unpalatability)]

## get tree subset (uses all 366 species in the phylogeny)
tree <- drop.tip(data_tree_chaz2021_all, setdiff(
  data_tree_chaz2021_all$tip.label,
  data_labels$species))

## trait probability matrix
data_labels_wide = dcast(
  data_labels[species %in% tree$tip.label],
  species ~ unpalatability,
  fun.aggregate = length,
  value.var = "unpalatability",
  fill = 0
)
data_labels_wide[
  , (c("high","medium","low")) := lapply(.SD, as.numeric), 
  .SDcols = c("high","medium","low")]
data_labels_wide[
  species %in% data_labels[!evidence_level %in% ev_lev]$species,
  c("high","low") := .(1/2, 1/2)
]
data_labels_wide[
  species %in% data_labels[evidence_level %in% ev_lev & unpalatability=="medium"]$species,
  c("high","low") := .(2/3, 1/3)
]
data_labels_wide[
  species %in% data_labels[evidence_level == "gn" & unpalatability=="high"]$species,
  c("high","low") := .(2/3, 1/3)
]
data_labels_wide[
  species %in% data_labels[evidence_level == "gn" & unpalatability=="low"]$species,
  c("high","low") := .(1/3, 2/3)
]

data_labels_wide[, medium := NULL]
labels_unpala_wide = as.matrix(data_labels_wide[, c("low","high"), with=F])
rownames(labels_unpala_wide) = data_labels_wide$species

table(data_labels_wide$high)
table(data_labels_wide$low)

## sim map - all
set.seed(42)
mod3 <- make.simmap(tree, labels_unpala_wide, model="ARD", nsim=1000, pi="fitzjohn")
tips <- mod3[[1]]$tip.label
Ntip <- length(tips)
tip_by_sim <- rbindlist(lapply(seq_along(mod3), function(s) {
  tr <- mod3[[s]]
  tip_edge_idx <- match(1:Ntip, tr$edge[,2])
  tip_state <- vapply(tip_edge_idx, function(i) tail(names(tr$maps[[i]]), 1), character(1))
  data.table(sim = s, species = tips, tip_state = tip_state)
}))
tip_post <- tip_by_sim[, .N, by=.(species, tip_state)]
tip_post[, prob := N / sum(N), by=species]
tip_post_wide <- dcast(tip_post, species ~ tip_state, value.var="prob", fill=0)
tip_post_wide[, unpalatability_est :=  c("low","high")[max.col(.SD, ties.method="first")], .SDcols =  c("low","high")]
tip_post_wide[, prob_est  := do.call(pmax, c(.SD, list(na.rm=TRUE))), .SDcols =  c("low","high")]

## save new label estimates
data_labels_new = merge(data_labels, tip_post_wide[, c("species","unpalatability_est", "prob_est")], all.x=T)
setnames(data_labels_new, "unpalatability", "unpalatability_orig")
fwrite(data_labels_new, "data/data_secondary/labels_imputed.csv")

# > labels - format meta data  --------------------------------------------------

## imputed labels
data_labels_new = fread("data/data_secondary/labels_imputed.csv", na.strings = "", stringsAsFactors = F)
data_labels_new[, unpalatability := unpalatability_est]

## format per-img metadata
data_meta_all = merge(data_meta, data_labels_new[, c("species", "subfamily", "tribe","unpalatability")], all.x=T)
data_meta_all[is.na(unpalatability), unpalatability := "unknown"]
data_meta_all[is.na(sex) | sex=="", sex := "unknown"]
data_meta_all[, unpalatability := factor(unpalatability, levels=c("unknown","low","medium","high"))]
data_meta_all[, species_label := ifelse(species %in% data_tree_chaz2021_all$tip.label, "label", "no-label")]

## format per-uuid metadata
data_meta_uuid = unique(data_meta_all[, c("uuid", "species", "subfamily", "tribe", "unpalatability", "sex", "species_label")])
data_meta_uuid = data_meta_uuid[order(species_label, species, uuid)]
data_meta_uuid[, row := 1:.N,]

## format aggregated metadata
data_meta_agg = unique(data_meta_uuid[, c("species", "species_label","subfamily", "tribe", "unpalatability")])
data_meta_agg = merge(data_meta_agg, data_meta_uuid[, .N, by=.(species)])

## moths order
data_meta_moths[, unpalatability := factor(unpalatability, levels=c("low","high"))]


# > labels - plot ------------------------------------------------------

## optional: plot imputed labels for chemical defense

# ## prep
# colors_unpala_tree <- c(
#   a = "#1b9e77",
#   b = "#7570b3",
#   c = "#d95f02")
# labels_unpala_tree <- c(
#   a = "Low",
#   b = "Medium",
#   c = "High")
# tree <- drop.tip(data_tree_chaz2021_all, setdiff(
#   data_tree_chaz2021_all$tip.label,
#   intersect(data_meta_all$species, data_labels$species)))
# ladderize = F
# 
# ## tree heatmap
# data_labels_new_sub = data_labels_new[species %in% tree$tip.label,]
# data_labels_new_sub[!is.na(evidence_level), unpalatability_certain := unpalatability_orig]
# unpalatability_labs <- data.frame(data_labels_new_sub[, c("unpalatability_certain","unpalatability_est")])
# unpalatability_labs[] = lapply(
#   unpalatability_labs,
#   factor,
#   levels = c("low","medium","high", "low_uncertain"),
#   labels = c("a","b","c","d")
# )
# rownames(unpalatability_labs) <- data_labels_new_sub$species
# 
# ## add subfamily info
# tree_phylo <- as.phylo(tree)
# species_by_family <- data_meta_agg[species %in% tree$tip.label, .(species_list = list(species)), by = "subfamily"]
# species_by_family <- setNames(species_by_family$species_list, species_by_family$subfamily)
# tree_phylo = groupOTU(tree_phylo, species_by_family, group_name = "subfamily")
# 
# ## format tree data
# data_tree = data.table(ggtree(tree_phylo , ladderize = ladderize)$data)[isTip == TRUE,]
# data_tree = data_tree[order(y)]
# setnames(data_tree, "label", "species")
# data_tree[, subfamily := factor(subfamily, levels=c(
#   "nymphalinae", "cyrestinae", "biblidinae", "apaturinae", "pseudergolinae",
#   "heliconiinae",  "limenitidinae", "satyrinae", "charaxinae", "libytheinae", "danainae"))]
# data_tree = merge(data_tree, data_labels_new_sub)
# 
# p1 =
#   ggtree(tree_phylo, ladderize = ladderize) +
#   geom_tree(aes(colour = subfamily), size = 1) +
#   coord_cartesian(clip="off") +
#   scale_color_manual(values=fam_cols,guide="none") +
#   
#   theme_void() + theme(
#     legend.position = "bottom",
#     plot.margin = margin(t=120,b=10,l=15,r=100, unit="pt")
#   )
# 
# p1 <- p1 + new_scale_fill()
# p1 = gheatmap(p1, unpalatability_labs[c("unpalatability_certain", "unpalatability_est")], 
#               offset=-1.75,  width=0.2, colnames = T, colnames_position = "top",
#               font.size = 4, family="bold", colnames_angle=90, hjust=0, colnames_offset_y = 1,
#               custom_column_labels=c("Labels (evidence)", "Labels (imputed)")
# ) +
#   scale_fill_manual(values=colors_unpala_tree, guide="none", na.value = colors_unpala["unknown"]) +
#   geom_text(data=data_tree, aes(x=x+20, y=y, label=species),
#             size=2.5, color="black", hjust=0) +
#   
#   geom_text(data=data_tree,
#             aes(x=x+14, y=y, label=round(prob_est,2), alpha=-prob_est),
#             size=1.5, hjust=0, show.legend = F,  fontface = "bold") +
#   scale_alpha_continuous(
#     range = c(0.2, 1), limits = c(-0.8, -0.5), oob = scales::squish
#   )
# 
# ## legend
# p_leg1 =
#   ggplot(tree_phylo) +
#   
#   geom_line(
#     aes(x=1,y=1,color=subfamily)) +
#   scale_color_manual(values=fam_cols, "Subfamily", guide=guide_legend(
#     order=1, override.aes = list(linewidth=1), keywidth = unit(15, "pt"),
#     title.position = "top", ncol=4, byrow=T),
#     labels=fam_labs) +
#   
#   new_scale_fill() +
#   geom_bar(
#     data=unpalatability_labs,
#     aes(x=1,y=1,fill=unpalatability_certain), stat="identity") +
#   scale_fill_manual(
#     "Chemical defense",
#     values=colors_unpala_tree,
#     labels=labels_unpala_tree,
#     guide=guide_legend(order=2, title.position = "top")) +
#   
#   theme(
#     legend.box = "vertical",
#     legend.box.just = "left",
#     legend.position = "bottom",
#     legend.box.margin = margin(b = 50, l=50, unit="pt"),
#   )
# 
# leg1 <- get_plot_component(p_leg1, 'guide-box', return_all = TRUE)[[3]]
# pg_final = cowplot::plot_grid(p1,leg1, ncol=1, rel_heights = c(0.95,0.05))
# 
# width=2000
# height=10000
# ggsave(paste0("figures/labels_impute_check.png"), pg_final, bg="white",
#        width=width, height=height, units = c("px"))

# > PCA -----------------------------------------------------------------

## prep
data_pca = merge(data_meta_all[
  ,c("mask_name", "uuid","species", "subfamily", "unpalatability", "class_dv", 
     "species_label", "sex")], data_emb, by="mask_name")
data_pca = data_pca[order(species, uuid, class_dv)]

## run
pca1 <- PCA(data_pca[, cols_emb, with=F], ncp=Inf, graph=FALSE, scale.unit=F)

## select components that expl. x% variance
pc_cutoff = which(pca1$eig[,3] >= 95)[1]
cols_pc_dims = str_remove(colnames(pca1$ind$coord[,1:pc_cutoff]), "\\.")
pca1_eigenvalues1 = data.table(pca1$eig)
setnames(pca1_eigenvalues1, c("eigenvalue","var_explained", "cum_var_explained"))

## dim reduction butterflies
pca_coords_butterflies = data.table(pca1$ind$coord[,1:pc_cutoff])
setnames(pca_coords_butterflies, cols_pc_dims)
pca_coords_butterflies = cbind(
  data_pca[, -cols_emb, with=F],
  pca_coords_butterflies)

## dim reduction moths
pca_coords_moths = predict(pca1, data_emb_moths[, cols_emb, with=FALSE])
pca_coords_moths = data.table(pca_coords_moths$coord[,1:pc_cutoff])
setnames(pca_coords_moths, cols_pc_dims)
pca_coords_moths = cbind(data_meta_moths, pca_coords_moths)

## aggregation
pca_coords_butterflies_agg = pca_coords_butterflies[
  ,lapply(.SD, mean), 
  .SDcols = cols_pc_dims, 
  by = c("species", "unpalatability")]
pca_coords_butterflies_agg_dv <- pca_coords_butterflies[
  , lapply(.SD, mean),
  .SDcols = cols_pc_dims, 
  by = c("species", "unpalatability", "class_dv")]
pca_coords_moths_agg = pca_coords_moths[
  ,as.list(unlist(lapply(.SD, mean))),
  .SDcols=cols_pc_dims,
  by=c("species", "unpalatability")]
pca_coords_moths_agg[,class_dv := "dorsal"]

# > LDA (w. CV) -------------------------------------------------------------

## select labels
species_incl = data_meta_agg[unpalatability %in% c("low","high"),]$species

## in-distribution (train + test)
pca_coords_butterflies_agg_incl_dv <- droplevels(pca_coords_butterflies_agg_dv[species %in% species_incl])

## create folds based on species identity
set.seed(42)
k_folds = 5
data_folds <- unique(droplevels(data_meta_agg[species %in% species_incl, .(species, unpalatability)]))
data_folds[, fold := createFolds(unpalatability, k=k_folds, list=FALSE)]
print(data_folds[, .N, by=.(fold, unpalatability)][order(fold)])

## run lda
pc_range <- 2:50
results_metrics <- data.table()
results_ld_scores <- data.table()

for (pc in pc_range) {
  message("pcs = ", pc)
  
  for (k in 1:k_folds) {
    train_species <- data_folds[fold != k, species]
    test_species  <- data_folds[fold == k, species]
    
    train_data <- pca_coords_butterflies_agg_incl_dv[species %in% train_species]
    test_data  <- pca_coords_butterflies_agg_incl_dv[species %in% test_species]
    datasets <- c(list("train"=train_data, 
                       "test"=test_data, 
                       "test_moths"=pca_coords_moths_agg))
    
    lda_mod <- lda(
      as.matrix(train_data[, cols_pc_dims[1:pc], with=FALSE]),
      as.factor(train_data$unpalatability)
    )
    
    for (dataset_name in names(datasets)) {
      ev <- lda_evaluate_cv(lda_mod, datasets[[dataset_name]], pc, cols_pc_dims)
      
      if (nrow(ev$metrics)) {
        results_metrics <- rbind(
          results_metrics,
          cbind(data.table(pcs=pc, fold=k, dataset=dataset_name), ev$metrics),
          fill=TRUE
        )
      }
      
      if (nrow(ev$scores)) {
        results_ld_scores <- rbind(
          results_ld_scores,
          cbind(data.table(pcs=pc, fold=k, dataset=dataset_name), ev$scores),
          fill=TRUE
        )
      }
    }
  }
}

# results_metrics: per (pcs, fold, dataset, class_dv) metrics
# results_ld_scores: per-species predictions + LD scores (x)
data_lda_results <- results_metrics[
  , .(
    accuracy  = mean(accuracy,  na.rm=TRUE),
    precision = mean(precision, na.rm=TRUE),
    recall    = mean(recall,    na.rm=TRUE),
    F1_score  = mean(F1_score,  na.rm=TRUE),
    accuracy_sd  = sd(accuracy,  na.rm=TRUE),
    F1_score_sd  = sd(F1_score,  na.rm=TRUE)
  ),
  by = .(pcs, dataset, class_dv)
]
fwrite(data_lda_results, "data/data_secondary/lda_results.csv")

## final components
final_components = 8
results_metrics[, dataset := factor(dataset, levels=c("train", "test", "test_moths"))]
ggplot(results_metrics[pcs<=50]) +
  facet_wrap(~dataset, ncol=4) + 
  geom_line(aes(x=pcs, y=accuracy,
                group=interaction(fold,class_dv),
                color=factor(fold),
                linetype=class_dv), linewidth=1) 

## confusion matrices 
species_level <- results_ld_scores[pcs == final_components, .(
  unpalatability = unique(unpalatability)[1],
  LD1_mean = mean(LD1, na.rm=TRUE),
  LD1_sd   = sd(LD1,   na.rm=TRUE),
  pred_vote = names(which.max(table(pred)))
), by=.(dataset, class_dv, species)]
species_level[, unpalatability := factor(unpalatability, levels=c("low","high"))]
species_level[, pred_vote := factor(pred_vote, levels=c("low","high"))]
conf_mats_clean <- species_level[
  !is.na(unpalatability),
  .(cm = list(table(
    Predicted = pred_vote,
    Actual    = unpalatability
  ))),
  by=.(dataset, class_dv)
]
for (i in seq_len(nrow(conf_mats_clean))) {
  cat("\n", conf_mats_clean$dataset[i], "|", conf_mats_clean$class_dv[i], "\n")
  print(conf_mats_clean$cm[[i]])
}

lda_folds <- lapply(1:k_folds, function(k){
  tr_sp <- data_folds[fold != k, species]
  tr <- droplevels(pca_coords_butterflies_agg[species %in% tr_sp])
  lda(as.matrix(tr[, cols_pc_dims[1:final_components], with=FALSE]),
      as.factor(tr$unpalatability))
})

save(lda_folds, file="data/analyses_primary/lda_ensemble_model.RData")

# > predict LD (ensemble) ----------------------------------------------

load(file="data/analyses_primary/lda_ensemble_model.RData")

ens_b <- lda_ensemble_predict(lda_folds, pca_coords_butterflies, cols_pc_dims, final_components)
data_ld_scores_butterflies <- data.table(
  pca_coords_butterflies[, -cols_pc_dims, with=FALSE],
  LD1 = ens_b$LD1_mean,
  LD1_sd = ens_b$LD1_sd,
  unpalatability_pred = ens_b$pred_class,
  prob_high = ens_b$post_mean[, "high"]
)

ens_m <- lda_ensemble_predict(lda_folds, pca_coords_moths, cols_pc_dims, final_components)
data_ld_scores_moths <- data.table(
  pca_coords_moths[, c("species","activity","diet","strategy","unpalatability")],
  LD1 = ens_m$LD1_mean,
  LD1_sd = ens_m$LD1_sd,
  unpalatability_pred = ens_m$pred_class,
  prob_high = ens_m$post_mean[, "high"]
)

## LD cutoff estimate
data_scores <- copy(data_ld_scores_butterflies)[unpalatability_pred %chin% c("low","high")]
setorder(data_scores, LD1)
data_scores[, prev := data.table::shift(unpalatability_pred)]
data_scores[, prev_LD1 := data.table::shift(LD1)]
data_scores[!is.na(prev) & unpalatability_pred != prev, mean(c(prev_LD1, LD1))]

# > table S2 (confusion mat.) ----------------------------------------------------------------

## fold-level confusion matrices for final_components
cms_by_fold <- results_ld_scores[
  pcs == final_components & !is.na(unpalatability) & !is.na(pred),
  {
    y_true <- factor(unpalatability, levels=c("low","high"))
    y_pred <- factor(pred, levels=c("low","high"))
    
    tab <- table(Predicted=y_pred, Actual=y_true)
    
    .(
      low_low   = tab["low",  "low"],
      high_low  = tab["high", "low"],
      low_high  = tab["low",  "high"],
      high_high = tab["high", "high"]
    )
  },
  by = .(pcs, fold, dataset, class_dv)
]

cms_mean <- cms_by_fold[
  ,
  .(
    low_low   = mean(low_low,   na.rm=TRUE),
    high_low  = mean(high_low,  na.rm=TRUE),
    low_high  = mean(low_high,  na.rm=TRUE),
    high_high = mean(high_high, na.rm=TRUE),
    
    low_low_sd   = sd(low_low,   na.rm=TRUE),
    high_low_sd  = sd(high_low,  na.rm=TRUE),
    low_high_sd  = sd(low_high,  na.rm=TRUE),
    high_high_sd = sd(high_high, na.rm=TRUE)
  ),
  by = .(pcs, dataset, class_dv)
]

for (i in seq_len(nrow(cms_mean))) {
  
  cat("\n", cms_mean$dataset[i], "|", cms_mean$class_dv[i], "\n")
  
  tab <- matrix(
    c(
      cms_mean$low_low[i],
      cms_mean$high_low[i],
      cms_mean$low_high[i],
      cms_mean$high_high[i]
    ),
    nrow = 2,
    byrow = FALSE,
    dimnames = list(
      Predicted = c("low","high"),
      Actual = c("low","high")
    )
  )
  
  print(round(tab, 0))
}

# > PCA + supp ---------------------------------------------------------

# prep
data_feat_scaled = copy(data_feat)
data_feat_scaled[, (cols_feat) := lapply(.SD, scale), .SDcols=cols_feat]
data_pca2 = merge(data_emb, data_feat_scaled)
data_pca2 = merge(data_pca2, data_ld_scores_butterflies[,c("mask_name", "LD1")])

## run
pca2 <- PCA(data_pca2[, c(cols_emb, cols_feat, "LD1"), with=F], ncp=Inf, graph=FALSE, scale.unit=F, quanti.sup=769:925)

## pca cols
cols_pc_dims_all = paste0("Dim", 1:pc_cutoff)
cols_pc_dims_apo = paste0("Dim", 1:final_components)

## eigenvalues
pca_eigenvalues = data.table(pca2$eig)
setnames(pca_eigenvalues, c("eigenvalue","var_explained", "cum_var_explained"))

## coordinates
pca_coords_butterflies = data.table(pca2$ind$coord[,1:pc_cutoff])
setnames(pca_coords_butterflies, cols_pc_dims_all)
pca_coords_butterflies = cbind(mask_name=data_pca2$mask_name, pca_coords_butterflies)
pca_coords_butterflies = merge(data_meta_all, pca_coords_butterflies, by="mask_name")

## suporting variables
data_pca_quanti_supp = data.table(
  feature_string=rownames(pca2$quanti.sup$coord),
  Dim1=pca2$quanti.sup$coord[,1],
  Dim2=pca2$quanti.sup$coord[,2])
data_pca_quanti_supp = merge(feature_key, data_pca_quanti_supp, all=T)

# > diff metrics  ------------------------------------------------

pca_coords_butterflies_agg = pca_coords_butterflies[
  ,lapply(.SD, mean), 
  .SDcols = cols_pc_dims_all, 
  by = c("species")]
pca_coords_butterflies_agg_dv = pca_coords_butterflies[
  ,lapply(.SD, mean), 
  .SDcols = cols_pc_dims_all, 
  by = c("species", "class_dv")]
pca_coords_butterflies_agg_sex <- pca_coords_butterflies[
  !sex=="unknown",
  lapply(.SD, mean),
  .SDcols = cols_pc_dims_all, 
  by = c("species", "sex")]
pca_coords_butterflies_agg_dv_sex <- pca_coords_butterflies[
  !sex=="unknown",
  lapply(.SD, mean),
  .SDcols = cols_pc_dims_all, 
  by = c("species", "class_dv", "sex")]
pca_coords_butterflies_agg_dv_sex = pca_coords_butterflies_agg_dv_sex[
  species %in% pca_coords_butterflies_agg_dv_sex[, .N, by="species"][N==4]$species]

# dv diff
pca_coords_butterflies_agg_dv_wide = dcast(pca_coords_butterflies_agg_dv, species~class_dv, value.var = cols_pc_dims_all)
pca_coords_butterflies_agg_dv_wide[, dv_diff_apo := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_dorsal"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_ventral"), with = FALSE])
), by = c("species")]
pca_coords_butterflies_agg_dv_wide[, dv_diff_all := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_dorsal"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_ventral"), with = FALSE])
), by = c("species")]
data_dv_diff_agg = pca_coords_butterflies_agg_dv_wide[,c("species", "dv_diff_apo", "dv_diff_all")]

## dv diff by sex
pca_coords_butterflies_agg_dv_sex_wide = dcast(pca_coords_butterflies_agg_dv_sex, species+sex~class_dv, value.var = cols_pc_dims_all)
pca_coords_butterflies_agg_dv_sex_wide[, dv_diff_apo := euclid_dist( 
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_dorsal"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_ventral"), with = FALSE])
), by = c("species", "sex")]
pca_coords_butterflies_agg_dv_sex_wide[, dv_diff_all := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_dorsal"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_ventral"), with = FALSE])
), by = c("species", "sex")]
data_dv_diff_agg_sex = pca_coords_butterflies_agg_dv_sex_wide[,c("species", "sex","dv_diff_apo", "dv_diff_all")]
data_dv_diff_agg_sex <- data_dv_diff_agg_sex[
  species %in% data_dv_diff_agg_sex[, .N, by = species][N == 2, species]
]

## sex diff
pca_coords_butterflies_agg_sex_wide = dcast(pca_coords_butterflies_agg_sex, species~sex, value.var = cols_pc_dims_all)
pca_coords_butterflies_agg_sex_wide[, sex_diff_apo := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_male"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_female"), with = FALSE])
), by = c("species")]
pca_coords_butterflies_agg_sex_wide[, sex_diff_all := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_male"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_female"), with = FALSE])
), by = c("species")]
data_sex_diff_agg = pca_coords_butterflies_agg_sex_wide[,c("species", "sex_diff_apo", "sex_diff_all")]
data_sex_diff_agg = data_sex_diff_agg[!is.na(sex_diff_apo)]

## sex diff by dv
pca_coords_butterflies_agg_sex_dv_wide = dcast(pca_coords_butterflies_agg_dv_sex, species+class_dv~sex, value.var = cols_pc_dims_all)
pca_coords_butterflies_agg_sex_dv_wide[, sex_diff_apo := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_male"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:final_components, "_female"), with = FALSE])
), by = c("species", "class_dv")]
pca_coords_butterflies_agg_sex_dv_wide[, sex_diff_all := euclid_dist(
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_male"), with = FALSE]),
  as.numeric(.SD[1, paste0("Dim", 1:pc_cutoff, "_female"), with = FALSE])
), by = c("species", "class_dv")]
data_sex_diff_agg_dv = pca_coords_butterflies_agg_sex_dv_wide[,c("species", "class_dv", "sex_diff_apo", "sex_diff_all")]
data_sex_diff_agg_dv <- data_sex_diff_agg_dv[
  species %in% data_sex_diff_agg_dv[!is.na(sex_diff_apo), .N, by = species][N == 2, species]
]

## mean diff
cols_pc_dims_apo = paste0("Dim", 1:final_components)
pca_coords_butterflies_sub = pca_coords_butterflies[!sex=="unknown" ]
pca_coords_butterflies_sub_n = pca_coords_butterflies_sub[, .N, by=c("species", "sex", "class_dv")]
pca_coords_butterflies_sub_n = pca_coords_butterflies_sub_n[N>=3]
pca_coords_butterflies_sub = pca_coords_butterflies_sub[
  species %in% pca_coords_butterflies_sub_n[, .N, by="species"][N==4]$species]
data_mean_diff_agg_dv_sex <- pca_coords_butterflies_sub[, {
  mat_7   <- as.matrix(.SD[,  paste0("Dim", 1:final_components), with=F])
  mat_177 <- as.matrix(.SD[,  paste0("Dim", 1:pc_cutoff), with=F])
  .(mean_diff_apo = mean_euclid_dist(mat_7), mean_diff_all = mean_euclid_dist(mat_177))
}, by = c("species", "class_dv", "sex")]


# EXPORT ------------------------------------------------------------------

## export all processed data for the next step 

fwrite(pca_coords_butterflies, "data/data_secondary/pca2_coords_butterflies.csv")
fwrite(pca_eigenvalues, "data/data_secondary/pca2_eigenvalues.csv")
fwrite(data_pca_quanti_supp, "data/data_secondary/pca2_quanti_supp.csv")
fwrite(pca_coords_moths_agg, "data/data_secondary/pca1_coords_moths_agg.csv")

fwrite(data_ld_scores_butterflies, "data/data_secondary/ld_scores_butterflies.csv")
fwrite(data_ld_scores_moths, "data/data_secondary/ld_scores_moths.csv")

fwrite(data_dv_diff_agg, "data/data_secondary/dv_diff_agg.csv")
fwrite(data_dv_diff_agg_sex, "data/data_secondary/dv_diff_agg_sex.csv")
fwrite(data_sex_diff_agg, "data/data_secondary/sex_diff_agg.csv")
fwrite(data_sex_diff_agg_dv, "data/data_secondary/sex_diff_agg_dv.csv")
fwrite(data_mean_diff_agg_dv_sex, "data/data_secondary/mean_diff_agg_dv_sex.csv")

# ROIs --------------------------------------------------------------------

## this pulls the butterfly images closest to the respective species centroids in PCA space
## needs to be generated before making the figures in the next script

## export all 
pca_coords_butterflies[, roi_path := paste0( "data_raw//segmentation_masks//butterflies_all_clean", "//" , species, "//", mask_name)]

## class dv
pca_coords_butterflies_summ = pca_coords_butterflies[, lapply(.SD, mean),.SDcols = cols_pc_dims_all, by=c("species", "class_dv")]
setnames(pca_coords_butterflies_summ, cols_pc_dims_all, paste0(cols_pc_dims_all, "_c"))
pca_coords_butterflies_all = merge(pca_coords_butterflies, pca_coords_butterflies_summ, by=c("species", "class_dv"), all.x=T)
pca_coords_butterflies_all[, dist := sqrt(rowSums((.SD - pca_coords_butterflies_all[
  , paste0(cols_pc_dims_all, "_c"), with=FALSE])^2)), .SDcols = cols_pc_dims_all]
data_sel <- pca_coords_butterflies_all[pca_coords_butterflies_all[, .I[which.min(dist)], by=c("species", "class_dv")]$V1]
data_sel = data_sel[,-c(cols_pc_dims_all, paste0(cols_pc_dims_all, "_c")), with=F]
data_sel[, new_path := file.path(paste0("data_raw//segmentation_masks//butterflies_centroids_dv//", class_dv, "//", species, ".png"))]
data_sel[, roi_path_exists := file.exists(roi_path)]
data_sel[roi_path_exists == TRUE, {
  mapply(function(from, to) {
    dir.create(dirname(to), recursive = TRUE, showWarnings = FALSE)
    fs::file_copy(from, to, overwrite = TRUE)
  }, roi_path, new_path)
}]

## sex
pca_coords_butterflies_summ = pca_coords_butterflies[
  !sex=="unknown" & class_dv=="dorsal", 
  lapply(.SD, mean),.SDcols = cols_pc_dims_all, by=c("species", "sex")]
setnames(pca_coords_butterflies_summ, cols_pc_dims_all, paste0(cols_pc_dims_all, "_c"))
pca_coords_butterflies_all = merge(pca_coords_butterflies, pca_coords_butterflies_summ, by=c("species", "sex"), all.x=T)
pca_coords_butterflies_all[, dist := sqrt(rowSums((.SD - pca_coords_butterflies_all[
  , paste0(cols_pc_dims_all, "_c"), with=FALSE])^2)), .SDcols = cols_pc_dims_all]
data_sel <- pca_coords_butterflies_all[pca_coords_butterflies_all[, .I[which.min(dist)], by=c("species", "sex")]$V1]
data_sel = data_sel[,-c(cols_pc_dims_all, paste0(cols_pc_dims_all, "_c")), with=F]
data_sel[, new_path := file.path(paste0("data_raw//segmentation_masks//butterflies_centroids_sex//", sex, "//", species, ".png"))]
data_sel[, roi_path_exists := file.exists(roi_path)]
data_sel[roi_path_exists == TRUE, {
  mapply(function(from, to) {
    dir.create(dirname(to), recursive = TRUE, showWarnings = FALSE)
    fs::file_copy(from, to, overwrite = TRUE)
  }, roi_path, new_path)
}]

## class dv x sex
pca_coords_butterflies_summ = pca_coords_butterflies[!sex=="unknown", lapply(.SD, mean),.SDcols = cols_pc_dims_all, by=c("species", "class_dv", "sex")]
setnames(pca_coords_butterflies_summ, cols_pc_dims_all, paste0(cols_pc_dims_all, "_c"))
pca_coords_butterflies_all = merge(pca_coords_butterflies, pca_coords_butterflies_summ, by=c("species", "class_dv", "sex"), all.x=T)
pca_coords_butterflies_all[, dist := sqrt(rowSums((.SD - pca_coords_butterflies_all[
  , paste0(cols_pc_dims_all, "_c"), with=FALSE])^2)), .SDcols = cols_pc_dims_all]
data_sel <- pca_coords_butterflies_all[pca_coords_butterflies_all[, .I[which.min(dist)], by=c("species", "class_dv", "sex")]$V1]
data_sel = data_sel[,-c(cols_emb, paste0(cols_emb, "_c")), with=F]
data_sel[, new_path := file.path(paste0("data_raw//segmentation_masks//butterflies_centroids_dv_sex//",class_dv, "_", sex, "//", species, ".png"))]
data_sel[, roi_path_exists := file.exists(roi_path)]
data_sel[roi_path_exists == TRUE, {
  mapply(function(from, to) {
    dir.create(dirname(to), recursive = TRUE, showWarnings = FALSE)
    fs::file_copy(from, to, overwrite = TRUE)
  }, roi_path, new_path)
}]