library(phangorn) library(devtools) library(NELSI) library(phytools) library(phylobase) library(ggplot2) ######READ_IN_TREES##### ######################## n_genes <- 434 species_tree <- read.tree("###TREE_FILE") node_vector <- seq(length(species_tree$tip.label)+1, length(species_tree$tip.label)+length(branching.times(species_tree)), 1) ######GET_SPECIES_PAIRS##### ############################ species_pairs <- vector("list", length(node_vector)) for (i in 1:length(node_vector)){ if (length(extract.clade(species_tree, node_vector[[i]])$tip.label) == 2){ species_pairs[[i]] <- extract.clade(species_tree, node_vector[[i]])$tip.label } } ######PROCESS_SPECIES_PAIRS################# ############################################ removal <- vector(mode="numeric", length = 0) for (i in 1:length(species_pairs)){ if (length(species_pairs[[i]]) > 0){ if (species_pairs[[i]][[1]] == "O_pteripes_MM_2409"){ removal <- append(removal, i) } if (species_pairs[[i]][[1]] == "O_sp_JAR_4375"){ removal <- append(removal, i) } } if (length(species_pairs[[i]]) == 0){ removal <- append(removal, i) } } species_pairs_final <- species_pairs[-removal] ######EXTRACT_LINEAGE_EFFECTS_TAXON_PAIRS##### ############################################## concatenated_species_pairs_lengths <- vector("list", length(species_pairs_final)) for (i in 1:length(concatenated_species_pairs_lengths)){ concatenated_species_pairs_lengths[[i]] <- extract.clade(species_tree, findMRCA(species_tree, species_pairs_final[[i]], type = c("node")))$edge.length } lineage_effects <- vector("list", length(species_pairs_final)) for (i in 1:length(lineage_effects)){ if (concatenated_species_pairs_lengths[[i]][[1]] > concatenated_species_pairs_lengths[[i]][[2]]){ lineage_effects[[i]] <- concatenated_species_pairs_lengths[[i]][[1]]/concatenated_species_pairs_lengths[[i]][[2]] } else if (concatenated_species_pairs_lengths[[i]][[1]] <= concatenated_species_pairs_lengths[[i]][[2]]){ lineage_effects[[i]] <- concatenated_species_pairs_lengths[[i]][[2]]/concatenated_species_pairs_lengths[[i]][[1]] } } ######ANALYSIS_OF_LOGNORMAL_DISTRIBUTION##### ############################################# pairwise_branch_lengths <- vector("list", 100000) for (i in 1:length(pairwise_branch_lengths)){ pairwise_branch_lengths[[i]] <- rlnorm(2, -0.01125, 0.587405/4) } simulated_lineage_effects <- vector("list", length(pairwise_branch_lengths)) for (i in 1:length(pairwise_branch_lengths)){ if (pairwise_branch_lengths[[i]][[1]] > pairwise_branch_lengths[[i]][[2]]){ simulated_lineage_effects[[i]] <- pairwise_branch_lengths[[i]][[1]]/pairwise_branch_lengths[[i]][[2]] } else if (pairwise_branch_lengths[[i]][[1]] < pairwise_branch_lengths[[i]][[2]]){ simulated_lineage_effects[[i]] <- pairwise_branch_lengths[[i]][[2]]/pairwise_branch_lengths[[i]][[1]] } }