### library(phangorn) library(devtools) library(NELSI) library(phytools) library(phylobase) rep_length <- 1 ### lineage_specific_rate <- 0.587405/4 gene_specific_rate <- 0.587405/4 residual_rate <- 0.587405/4 ### lineage_specific_rate_mean <- -0.01125 gene_specific_rate_mean <- -0.01125 residual_rate_mean <- -0.01125 ### lineage_specific_rates <- lineage_specific_rate residual_rates <- residual_rate gene_specific_rates <- gene_specific_rate lineage_specific_rates_means <- lineage_specific_rate_mean residual_rates_means <- residual_rate_mean gene_specific_rates_means <- gene_specific_rate_mean ### locus_number <- 400 sequence_length <- 800 for (a in 1:rep_length){ tree <- read.tree("three_taxon_tree.tre") setwd(directories[[a]]) file_names <- paste(seq(1, locus_number, 1),".nexus", sep = "") tree_file_names <- c(paste(seq(1, locus_number, 1),".tre", sep = "")) concatenated_file_names <- paste(c(1,5,20,400),"concatenated.nexus", sep = "") ###LINEAGE_EFFECT for (i in 1:length(tree[[2]])){ tree[[2]][[i]] <- tree[[2]][[i]] * rlnorm(1, meanlog = lineage_specific_rates_means[[a]], sdlog = lineage_specific_rates[[a]]) } ###GENERATE_GENE_TREES gene_trees <- vector("list", length(file_names)) for (i in 1:length(gene_trees)){ gene_trees[[i]] <- tree } ###GENE_EFFECT for (i in 1:length(gene_trees)){ gene_trees[[i]][[2]] <- gene_trees[[i]][[2]] * rlnorm(1, meanlog = gene_specific_rates_means[[a]], sdlog = gene_specific_rates[[a]]) } ###RESIDUAL_EFFECT for (j in 1:length(gene_trees)){ for (i in 1:length(gene_trees[[j]][[2]])){ gene_trees[[j]][[2]][[i]] <- gene_trees[[j]][[2]][[i]] * rlnorm(1, meanlog = residual_rates_means[[a]], sdlog = residual_rates[[a]]) } } ###GENERATE_SEQUENCES locus_sequences <- vector("list", length(file_names)) for (i in 1:length(file_names)){ locus_sequences[[i]] <- simSeq(gene_trees[[i]], sequence_length, Q=NULL, bf=NULL, rootseq=NULL, type="DNA", rate = 0.05) } ###WRITE_SEQUENCES_TO_FILE for (i in 1:length(file_names)){ write.phyDat(locus_sequences[[i]], file = file_names[[i]], format = "nexus") } ###WRITE_TREES_TO_FILE for (i in 1:length(tree_file_names)){ write.tree(gene_trees[[i]], file = tree_file_names[[i]]) } ###CONCATENATE_SEQEUNCES_FOR_FASTER_ANALYSIS sequences <- vector("list", locus_number) for (i in 1:locus_number){ sequences[[i]] <- read.nexus.data(file = file_names[[i]]) } concatenated_one <- vector("list", 3) concatenated_five <- vector("list", 3) concatenated_twenty <- vector("list", 3) concatenated_four_hundered <- vector("list", 3) for (i in 1:3){ concatenated_one[[i]] <- sequences[[1]][[i]] } concatenated_one <- phyDat(list(A = concatenated_one[[1]], B = concatenated_one[[2]], C = concatenated_one[[3]])) for (j in 1:3){ for (i in 1:5){ concatenated_five[[j]] <- append(concatenated_five[[j]], sequences[[i]][[j]]) } } concatenated_five <- phyDat(list(A = concatenated_five[[1]], B = concatenated_five[[2]], C = concatenated_five[[3]])) for (j in 1:3){ for (i in 1:20){ concatenated_twenty[[j]] <- append(concatenated_twenty[[j]], sequences[[i]][[j]]) } } concatenated_twenty <- phyDat(list(A = concatenated_twenty[[1]], B = concatenated_twenty[[2]], C = concatenated_twenty[[3]])) for (j in 1:3){ for (i in 1:400){ concatenated_four_hundered[[j]] <- append(concatenated_four_hundered[[j]], sequences[[i]][[j]]) } } concatenated_four_hundered <- phyDat(list(A = concatenated_four_hundered[[1]], B = concatenated_four_hundered[[2]], C = concatenated_four_hundered[[3]])) write.phyDat(concatenated_one, file = concatenated_file_names[[1]], format = "nexus") write.phyDat(concatenated_five, file = concatenated_file_names[[2]], format = "nexus") write.phyDat(concatenated_twenty, file = concatenated_file_names[[3]], format = "nexus") write.phyDat(concatenated_four_hundered, file = concatenated_file_names[[4]], format = "nexus") #### setwd("###LOCATION_OF_EXPERIMENT") remove("tree", "file_names", "gene_trees", "locus_sequences", "sequences", "concatenated_one", "concatenated_five", "concatenated_twenty", "concatenated_four_hundered") }