在R中高效计算氨基酸组合与OGT的系统发育相关性

Hey there! Let's work through this problem to fix the efficiency bottleneck and properly account for phylogenetic signal. Your original approach struggles with the sheer number of amino acid combinations (2²⁰ - 1 = 1,048,575 total non-empty combinations) because it relies on slow loop-based operations. We'll switch to vectorized computations and use phylogenetically informed statistical methods to get reliable results.

Step 1: Streamline Data Prep & Tree Matching #

First, we'll clean up the data and ensure our species align with the phylogenetic tree (critical for valid phylogenetic analyses):

# Load core packages
library(dplyr)
library(ape)
library(geiger)
library(caper)

# Import tree and data
taxonomy_tree <- read.nexus("taxonomyforzeldospecies.nex")
zeldodata <- read.csv("COMPLETECOPYFORR.csv")

# Keep only necessary columns: Species, OGT, and amino acids A-Y
aa_data <- zeldodata %>%
  select(Species, OGT, A:Y)

# Match tree and data (removes species missing from either)
matched_dataset <- treedata(taxonomy_tree, aa_data, sort = TRUE)
phylo_tree <- matched_dataset$phy
cleaned_data <- matched_dataset$data

Step 2: Optimize Subset Sum Calculation #

Instead of looping through every combination for each species, we'll use vectorized matrix operations to compute all subset sums in one pass. This leverages R's optimized C-level backend for speed:

# Generate all non-empty amino acid subsets (binary mask matrix)
n_amino_acids <- 20
subset_masks <- expand.grid(rep(list(c(0, 1)), n_amino_acids)) %>%
  filter(rowSums(.) > 0)  # Exclude empty subset

# Convert masks to a matrix (each column = one amino acid combination)
subset_matrix <- t(as.matrix(subset_masks))
# Name columns with the corresponding amino acid combinations (e.g., "AB", "ACD")
colnames(subset_matrix) <- apply(subset_masks, 1, function(x) {
  paste(names(aa_data)[3:22][x == 1], collapse = "")
})

# Calculate all subset sums for every species (vectorized operation)
subset_sums <- as.data.frame(cleaned_data[, 3:22] %*% subset_matrix)

# Combine with original Species and OGT columns
full_analysis_data <- cbind(cleaned_data[, 1:2], subset_sums)

This cuts down computation time drastically—matrix multiplication is orders of magnitude faster than looping through each combination for every species.

Step 3: Phylogenetically Informed Correlation Analysis #

Standard Pearson correlation ignores shared ancestry between species, which can lead to spurious results. We'll use Phylogenetic Generalized Least Squares (PGLS) to model the relationship between each subset sum and OGT while accounting for phylogenetic signal:

Option 1: Sequential Processing (Simple) #

# Create comparative data object for caper (links tree to data)
comp_data <- comparative.data(phylo_tree, full_analysis_data, names.col = "Species", vcv = TRUE)

# Function to run PGLS and extract key stats
get_pgls_results <- function(subset_col, comp_data) {
  model <- pgls(OGT ~ ., data = comp_data$data[, c("OGT", subset_col)], comp_data = comp_data)
  tibble(
    combination = subset_col,
    correlation_coef = coef(model)[[2]],
    p_value = summary(model)$coefficients[2, 4],
    phylogenetic_lambda = model$param  # Measures phylogenetic signal strength
  )
}

# Run analysis for all combinations
correlation_results <- lapply(colnames(full_analysis_data)[3:ncol(full_analysis_data)], 
                              get_pgls_results, comp_data = comp_data) %>%
  bind_rows()

# Sort results by absolute correlation (optional)
sorted_results <- correlation_results %>%
  arrange(desc(abs(correlation_coef)))

Option 2: Parallel Processing (Faster for 1M+ Combinations) #

If you need to speed up the PGLS step, use parallel computing to run multiple models at once:

library(parallel)

# Set up parallel cluster
cl <- makeCluster(detectCores() - 1)
clusterExport(cl, c("comp_data", "get_pgls_results"))
clusterEvalQ(cl, {
  library(caper)
  library(dplyr)
})

# Run parallel analysis
parallel_results <- parLapply(cl, colnames(full_analysis_data)[3:ncol(full_analysis_data)], 
                              get_pgls_results, comp_data = comp_data) %>%
  bind_rows()

# Stop cluster
stopCluster(cl)

Key Improvements Over Your Original Code #

1. Vectorized Computations: Replaces slow loops with matrix multiplication to compute all subset sums in seconds instead of hours.
2. Phylogenetic Correction: Uses PGLS to account for shared ancestry, ensuring your correlation results are biologically meaningful.
3. Data Alignment: The treedata function ensures only species present in both the tree and dataset are included, avoiding analysis errors.

内容的提问来源于stack exchange，提问作者Biomage