diagnostics_and_testing.Rmd

---
title: "Experiments with information theoretic tree generation"
output:
  github_document:
    pandoc_args: --webtex
---

```{r setup, include=FALSE}
knitr::opts_knit$set(root.dir = '/Users/shashanksule/Documents/info_theoretic_phylo/')
library("adegenet")
library("ape")
library("apTreeshape")
library("BoSSA")
library("diversitree")
library("pegas")
library("phangorn")
library("phylobase")
#library("phyloch")
library("seqinr")
library("readr")
source("utilities.R")
library("progress")
library("ggplot2")
library("phyclust")
library("TreeDist")
library("TreeTools")
library("dplyr")
library("parallel")
library("adephylo")
library("frequency")
library("plot3D")
```

This is the place where we'll be testing and exploring our divisive and agglomerative algorithms that make trees from aligned nucleotide sequences. Look for more extensive documentation on the computing process [here](https://github.com/ShashankSule/info_theoretic_phylo/blob/main/Computing_docs.md). 

# Simulated Sequences 

We generate 100 blocks of sequence data for 10 species where each block corresponds to a tree randomly generated in the TN93 model. Don't run this unless you want to generate new data! 

```{r echo=TRUE, results='hide', eval = FALSE}
#setwd("/Users/shashanksule/Documents/info_theoretic_phylo/")
n <- 1
trees <- list(n)
sequences <- list(n)
seeds <- list(n)
for(i in 1:n){
  
  #Initialize the .dat file and store the seed
  seeds[[i]] <- who_dat(seqs = 16,
                sites = 1000, 
                model = 6, 
                parameters = "5 5", 
                gamma = "0.5 4", 
                mutation = 5.5,
                equilibrium = "0.25 0.28 0.34 0.13",
                spit_seed = TRUE)
  
  #Simulate the trees and the sequences
  system("paml4.8/src/evolverRandomTree 5 paml4.8/MCbaseRTree.dat") 
  trees[[i]] <- write.tree(read.nexus("mctrees.nex"))
  sequences[[i]] <- ReadCharacters("mc.nex")
}
```

## Cummings et al 2003 parameters with TN93, 16 species, and 1000 sites

```{r echo=TRUE, results='hide', eval = FALSE}
num_sequence <- 100
trees <- list(num_sequence)
divisive <- list(num_sequence)
sequences <- list(num_sequence)
# agg_trees_nuc <- list(num_sequence)
# nj_trees_nuc <- list(num_sequence)

#read the original trees
for(i in c(1:num_sequence)){
  filename <- paste("./TN93_Jul21/Trees/tree", i, ".txt", sep = "")
  trees[[i]] <- readChar(filename, file.info(filename)$size)
}

#read the divisive trees
for(i in c(1:num_sequence)){
  filename <- paste("./TN93_Jul21/DivisiveTrees/tree", i, ".txt", sep = "")
  divisive[[i]] <- readChar(filename, file.info(filename)$size)
}
divisive <- lapply(divisive,  function(x) paste(x,";", sep = ""))

for(i in c(1:num_sequence)){
  filename <- paste("./TN93_Jul21/Sequences/sequence", i, ".nex", sep = "")
  sequences[[i]] <- ReadCharacters(filename)
}
```

For the tree that works, the min edge to mutation ratio is 0.01454545 under TN93 assumptions with mutation 5.5; the minimax is ge. 

For a tree that doesn't work, the min edge to mutation ratio is 0.002945818 under TN93 assumptions with mutation 5.5; the minimax is 0.003957794

Let's compute some diagnostics! 

```{r, echo = TRUE, results= 'hide'}
nj_trees_nuc <- sequences %>%
            lapply(as.DNAbin) %>%
            lapply(dist.dna, model = "TN93", base.freq = c(0.34, 0.28, 0.13, 0.25)) %>%
            lapply(nj) %>%
            lapply(write.tree)
```

1. Robinson-Foulds distances: 

```{r include=FALSE}
data <- data.frame(as.character(trees), as.character(divisive), as.character(nj_trees_nuc))
div_distance  <- mapply(RobinsonFoulds,
       lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
       lapply(data$as.character.divisive., function(x) read.tree(text = as.character(x))))
# agg_distance <- mapply(RobinsonFoulds,
#        lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
#        lapply(data$as.character.agglomerative., function(x) read.tree(text = as.character(x))))
nj_distance <- mapply(RobinsonFoulds,
       lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
       lapply(data$as.character.nj_trees_nuc., function(x) read.tree(text = as.character(x))))
# diagnostics <- data.frame(as.character(div_distance),
#                           as.character(agg_distance),
#                           as.character(nj_distance))
```

```{r, echo = FALSE}


# div_mean <- diagnostics$as.character.div_distance. %>%
#             as.character() %>%
#             as.numeric() %>%
#             mean()
# agg_mean <- diagnostics$as.character.agg_distance. %>% 
#             as.character() %>% 
#             as.numeric() %>% 
#             mean()
# nj_mean <-  diagnostics$as.character.nj_distance. %>% 
#             as.character() %>% 
#             as.numeric() %>% 
#             mean()

print(div_distance)
print(paste("Mean distance of divisive trees from original: ", mean(div_distance), sep = ""))
print(paste("The proportion of divisive trees has 0 RF distance from the original tree is: ", sum(div_distance == 0) / num_sequence, sep = ""))
hist(div_distance)

print(nj_distance)
print(paste("Mean distance of nj trees from original: ", mean(nj_distance), sep = ""))
print(paste("The proportion of nj trees has 0 RF distance from the original tree is: ", sum(nj_distance == 0) / num_sequence, sep = ""))
hist(nj_distance)
summary(nj_distance)
var(nj_distance)
# print(paste("Mean distance of divisive trees from original: ", div_mean, sep = ""))
# print(paste("Mean distance of agglomerative trees from original: ", agg_mean, sep = ""))
# print(paste("Mean distance of NJ trees from original: ", nj_mean, sep = ""))
```

2. Check the rootings using all.equals

```{r, echo = FALSE}
div_root  <- mapply(all.equal,
       lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
       lapply(data$as.character.divisive., function(x) read.tree(text = as.character(x))),
       use.edge.length = FALSE)
# agg_root <- mapply(all.equal,
#        lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
#        lapply(data$as.character.agglomerative., function(x) read.tree(text = as.character(x))),
#        use.edge.length = FALSE)
nj_root <- mapply(all.equal,
       lapply(data$as.character.trees., function(x) read.tree(text = as.character(x))),
       lapply(data$as.character.nj_trees_nuc., function(x) read.tree(text = as.character(x))),
       use.edge.length = FALSE)
```

```{r, echo = FALSE}
print(div_root)
# print(agg_root)
print(nj_root)

summary(div_root)
# summary(agg_root)
summary(nj_root)
```


3. Low tolerance ultrametricity: 

```{r include=FALSE}
metricity <- data[,c(1,2,3)]
metricity[] <- metricity %>%
               lapply(function(x) is.ultrametric(read.tree(text = as.character(x))))

```

```{r}
summary(metricity)
```

4. High tolerance ultrametricity: 

```{r include=FALSE}
low_metricity <- data[,c(1,2,3)]
low_metricity[] <- low_metricity %>%
                   lapply(function(x) is.ultrametric(read.tree(text = as.character(x)), tol = 1e-6))
```

```{r}
summary(low_metricity)
print("Tolerance used: 1e-6")
```

5. Check additivity: 

```{r include=FALSE}
# original_met <- logical(length = num_sequence)
original_add <- logical(length = num_sequence)
# div_met <- logical(length = num_sequence)
div_add <- logical(length = num_sequence)
# agg_met <- logical(length = num_sequence)
# agg_add <- logical(length = num_sequence)
# # nj_met <- logical(length = num_sequence)
nj_add <- logical(length = num_sequence)

for(i in c(1:num_sequence)){
  # original_met[i] <- is.ultrametric(read.tree(text = as.character(trees[i])))
  #original_add[i] <- is_additive(read.tree(text = as.character(trees[i])), 1e-15)
  original_add[i] <- is_additive(read.tree(text = as.character(trees[i])))
}

for(i in c(1:num_sequence)){
  # div_met[i] <- is.ultrametric(read.tree(text = as.character(divisive[i])))
  #div_add[i] <- is_additive(read.tree(text = as.character(divisive[i])), 1e-15)
  div_add[i] <- is_additive(read.tree(text = as.character(divisive[i])))
}

# for(i in c(1:num_sequence)){
#   # agg_met[i] <- is.ultrametric(read.tree(text = as.character(agglomerative[i])))
#   #agg_add[i] <- is_additive(read.tree(text = as.character(agglomerative[i])), 1e-15)
#   agg_add[i] <- is_additive(read.tree(text = as.character(agglomerative[i])))
# }
# 
for(i in c(1:num_sequence)){
  # nj_met[i] <- is.ultrametric(read.tree(text = as.character(nj_trees_nuc[i])))
  #nj_add[i] <- is_additive(read.tree(text = as.character(nj_trees_nuc[i])), 1e-15)
  nj_add[i] <- is_additive(read.tree(text = as.character(nj_trees_nuc[i])))
}
```

```{r}
# print(original_met)
# print(div_met)
# print(agg_met)
# print(nj_met)

# print(original_add)
# print(div_add)
# print(agg_add)
# print(nj_add)
summary(original_add)
summary(div_add)
summary(nj_add)
```


6. Relative distance between originals: 

```{r}
rel_distances <- matrix(rep(0,1e4), nrow = 100)
for(i in 1:99){
  for(j in (i+1):100)
 rel_distances[i,j] <- RobinsonFoulds(read.tree(text = trees[[i]]), read.tree(text = trees[[j]]))
}

```

```{r}
par(mar=c(5.1, 4.1, 4.1, 4.1))
plot(rel_distances[1:50,1:50])
```


Now we'll filter those that have zero robinson-foulds distances 

```{r}
index <- which(diagnostics$as.numeric.div_distance. + diagnostics$as.numeric.agg_distance. == 0)
```


Finally, we'll make a single plot with the branch lengths compared against the neighbour-joining ones: 

```{r}
# i <- index[1]
# nj_lengths <- read.tree(text = nj_trees_nuc[[i]])$edge.length
#     d_lengths <- read.tree(text = divisive[[i]])$edge.length
#     plot(nj_lengths, d_lengths)

cl <- rainbow(length(index))

for (i in index){
    nj_lengths <- read.tree(text = nj_trees_nuc[[i]])$edge.length
    d_lengths <- read.tree(text = divisive[[i]])$edge.length
    plot(nj_lengths[order(nj_lengths)], d_lengths[order(nj_lengths)])
}
```

## JC69 with 12 species and 2000 sites 

Reading the data and making the trees: 

```{r}
num_sequence <- 105
trees <- list(num_sequence)
agglomerative <- list(num_sequence)
divisive <- list(num_sequence)
sequences <- list(num_sequence)
# agg_trees_nuc <- list(num_sequence)
# nj_trees_nuc <- list(num_sequence)

#read the original trees
for(i in c(1:num_sequence)){
  filename <- paste("./JC69_Jul28/Trees/tree", i, ".txt", sep = "")
  trees[[i]] <- readChar(filename, file.info(filename)$size)
}

#read the divisive trees 
for(i in c(1:num_sequence)){
  filename <- paste("./JC69_Jul28/divisive_trees/tree", i, ".txt", sep = "")
  divisive[[i]] <- readChar(filename, file.info(filename)$size)
}
divisive <- lapply(divisive,  function(x) paste(x,";", sep = ""))

#read the agglomerative trees
for(i in c(1:num_sequence)){
  filename <- paste("./JC69_Jul28/agglomerative_trees/tree", i, ".txt", sep = "")
  agglomerative[[i]] <- readChar(filename, file.info(filename)$size)
}
agglomerative <- lapply(agglomerative,  function(x) paste(x,";", sep = ""))

for(i in c(1:num_sequence)){
  filename <- paste("./JC69_Jul28/Sequences/sequence", i, ".nex", sep = "")
  sequences[[i]] <- ReadCharacters(filename)
}

 #make the neighbour joining trees 

nj_trees_nuc <- sequences %>%
            lapply(as.DNAbin) %>%
            lapply(dist.dna, model = "JC69") %>%
            lapply(nj) %>%
            lapply(write.tree)

```

1. Robinson-Foulds Data: 

```{r}
data <- data.frame(Original = as.character(trees), 
                   Agglomerative = as.character(agglomerative), 
                   Divisive = as.character(divisive), 
                   NJ = as.character(nj_trees_nuc))
div_distance <- mapply(RobinsonFoulds,
       lapply(data$Divisive, function(x) read.tree(text = as.character(x))),
       lapply(data$Original, function(x) read.tree(text = as.character(x))))
agg_distance <- mapply(RobinsonFoulds,
       lapply(data$Original, function(x) read.tree(text = as.character(x))),
       lapply(data$Agglomerative, function(x) read.tree(text = as.character(x))))
nj_distance <- mapply(RobinsonFoulds,
       lapply(data$Original, function(x) read.tree(text = as.character(x))),
       lapply(data$NJ, function(x) read.tree(text = as.character(x))))

```

```{r}
layout(matrix(c(1,2,3,3), nrow = 2, ncol = 2))
hist(div_distance)
hist(agg_distance)
hist(nj_distance)
```

```{r}
dist_data <- data.frame(Divisive = div_distance, Agglomerative = agg_distance, NJ = nj_distance)
summary(dist_data)
```

```{r}
acc_div <- div_distance == 0
acc_agg <- agg_distance == 0 
acc_nj <- nj_distance == 0 
acc_data <- data.frame(Divisive = div_distance, Agglomerative = agg_distance, NJ = nj_distance)
table(acc_div)
table(acc_agg)
table(acc_nj)
```

```{r}
plot(dist_data$NJ, dist_data$Divisive)
```

```{r}
##  Create cuts:
x_c <- cut(dist_data$NJ, 2)
y_c <- cut(dist_data$Divisive, 2)

##  Calculate joint counts at cut levels:
z <- table(x_c, y_c)

##  Plot as a 3D histogram:
hist3D(z=z, border="black")
```


# Cummings et al 1995

```{r}
names <- c("atpase6",
"atpase8",
"coii",
"coiii",
"cytb",
"nd1",
"nd2",
"nd3",
"nd4",
"nd4l",
"nd5",
"nd6")

d_file_paths_nuc <- paste("./1995_trees/Nucleotide/Divisive/", names, ".txt", sep = "")
a_file_paths_nuc <- paste("./1995_trees/Nucleotide/Agglomerative/", names, ".txt", sep = "")
d_file_paths_cod <- paste("./1995_trees/Codon/Divisive/", names, ".txt", sep = "")

d_trees_nuc <- d_file_paths_nuc %>%
                  mclapply(readLines) %>%
                  mclapply(function(x) paste(x,";", sep = "")) %>%
                  mclapply(function(x) read.tree(text = x)) %>%
                  mclapply(ladderize)
a_trees_nuc <- a_file_paths_nuc %>%
                  mclapply(readLines) %>%
                  mclapply(function(x) paste(x,";", sep = "")) %>%
                  mclapply(function(x) read.tree(text = x)) %>%
                  mclapply(ladderize)
d_trees_cod <- d_file_paths_cod %>%
                  mclapply(readLines) %>%
                  mclapply(function(x) paste(x,";", sep = "")) %>%
                  mclapply(function(x) read.tree(text = x)) %>%
                  mclapply(ladderize)

```

## Nucleotide Data

1. Divisive

```{r}

gen_tree <- read.tree(text = "((Carp, Loach), (Xenopus, (Chicken, ((Mouse, Rat),(Human, (Seal, (Whale, Bovine)))))));")
plot(gen_tree)
title("Original tree")

#d_trees_nuc <- d_trees_nuc[-3]
for(i in c(1:12)){
  
  plot(d_trees_nuc[[i]])
  title(names[[i]])
  axisPhylo()
}
```

```{r}
d_distances <- as.numeric(mclapply(d_trees_nuc, function(x) RobinsonFoulds(gen_tree, x)))
d_metricity <- as.numeric(mclapply(d_trees_nuc, is.ultrametric))
d_additivity <- as.numeric(mclapply(d_trees_nuc, is_additive))
hist(d_distances, xlab = "Robinson Foulds Distance", ylab = "Frequency", main = 
       "Robinson Foulds distances of divisive trees \ncomputed from C. et al '95 datasets to original")

text(7.5,1.5, paste("Proportion of \nidentical trees:", sum(as.numeric(d_distances == 0))/13))
text(7.5,2.5, paste("Proportion of \nultrametric trees:", sum(as.numeric(d_metricity == 1))/13))
text(7.5,3.5, paste("Proportion of \nadditive trees:", sum(as.numeric(d_additivity == 1))/13))
```

2. Agglomerative Trees 

```{r}
plot(gen_tree)
title("Original tree")

#d_trees_nuc <- d_trees_nuc[-3]
for(i in c(1:12)){
  
  plot(a_trees_nuc[[i]])
  title(names[[i]])
  axisPhylo()
}
```






```{r}
a_distances <- as.numeric(mclapply(a_trees_nuc, function(x) RobinsonFoulds(gen_tree, x)))
a_metricity <- as.numeric(mclapply(a_trees_nuc, is.ultrametric))
a_additivity <- as.numeric(mclapply(a_trees_nuc, is_additive))
hist(a_distances, xlab = "Robinson Foulds Distance", ylab = "Frequency", main = 
       "Robinson Foulds distances of agglomerative trees \ncomputed from C. et al '95 datasets to original")

text(9.5,1.5, paste("Proportion of \nidentical trees:", sum(as.numeric(a_distances == 0))/12))
text(9.5,2.5, paste("Proportion of \nultrametric trees:", sum(as.numeric(a_metricity == 1))/12))
text(9.5,3.5, paste("Proportion of \nadditive trees:", sum(as.numeric(a_additivity == 1))/12))
```


## Codon 

1. Divisive Trees

```{r}
plot(gen_tree)
title("Original tree")

#d_trees_nuc <- d_trees_nuc[-3]
for(i in c(1:12)){
  
  plot(d_trees_cod[[i]])
  title(names[[i]])
  axisPhylo()
}
```

```{r}
d_distances <- as.numeric(mclapply(d_trees_cod, function(x) RobinsonFoulds(gen_tree, x)))
d_metricity <- as.numeric(mclapply(d_trees_cod, is.ultrametric))
d_additivity <- as.numeric(mclapply(d_trees_cod, is_additive))
hist(d_distances, xlab = "Robinson Foulds Distance", ylab = "Frequency", main = 
       "Robinson Foulds distances of divisive trees \ncomputed from C. et al '95 datasets to original")

text(7.5,1.5, paste("Proportion of \nidentical to genome trees:", sum(as.numeric(d_distances == 0))/13))
text(7.5,2.5, paste("Proportion of \nultrametric trees:", sum(as.numeric(d_metricity == 1))/13))
text(7.5,3.5, paste("Proportion of \nadditive trees:", sum(as.numeric(d_additivity == 1))/13))
```