---
title: "Phylogenies in R"
author: "Axel Arango"
date: "2025-10-16"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r,echo=FALSE,message=FALSE,warning=FALSE}
# Set so that long lines in R will be wrapped:
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80), tidy=TRUE)
```
## R Markdown
first we load packages
```{r}
library(muscle)
library(ape)
library(phangorn)
library(phytools)
```
Set our directory
```{r}
setwd("~/Documents/PostdocWue/Phylo_class/Code/datasets/")
```

load sequence data
```{r,}
adn<-read.dna("https://raw.githubusercontent.com/AxelArango/phyloclass/refs/heads/main/datasets/multitaxa.fna", format="fasta")
adn
```

read with strings
```{r}
dna<-readDNAStringSet("https://raw.githubusercontent.com/AxelArango/phyloclass/refs/heads/main/datasets/multitaxa.fna",format = "fasta")
dna
```

aligmnet with muscle
```{r}
dna_align<-muscle::muscle(dna)
dna_align
```
write aligned fasta
first we convert, then we write it
```{r}
dna_align<-as(dna_align,"DNAStringSet")
writeXStringSet(dna_align,file="multitaxa_a.fna")
```

now we load the aligned fasta
```{r}
dna<-read.dna("multitaxa_a.fna",format="fasta")
```

We can now calculate the genetic distance matrix
```{r}
dist_mat<-dist.dna(dna,model="K80")
head(dist_mat)
```

build a tree with upgma
```{r}
upgma_tree<-upgma(dist_mat,method="average")
plot(upgma_tree)
```

change names

```{r}
upgma_tree$tip.label<-c("Human","Chimpancee","Cow","Dog","Deer","Sheep","Falcon","Salmon","Bee","Bear","Gorilla")
plot(upgma_tree)

```

```{r}
nj_tree<-nj(dist_mat)
plot(nj_tree)
```

for visualization

```{r}
nj_tree2<-nj_tree
nj_tree2$tip.label<-upgma_tree$tip.label
plot(nj_tree2)

```

Maximum parsimony tree
first convert to phyDat
```{r}
phy_dat<-phyDat(dna,type="DNA")
phy_dat
```

lets create a start tree
```{r}
start_tree<-random.addition(phy_dat)
start_treeb<-start_tree;start_treeb$tip.label<-upgma_tree$tip.label
plot(start_treeb)
```

optimize topology
```{r}
pars_tree<-optim.parsimony(start_tree,phy_dat)
pars_tree2<-pars_tree;pars_tree2$tip.label<-upgma_tree$tip.label
plot(pars_tree2)
```
check parsimony score
```{r}
parsimony(pars_tree,phy_dat)
parsimony(start_tree,phy_dat)

```
search for more (and better trees)
```{r}
pars_tree_b<-pratchet(phy_dat,start_tree,minit = 10000,rearrangements = "SPR")
```

lets see it rooted
```{r}
rooted_pars<-root(pars_tree2,outgroup = "Bee",resolve.root = T)
plot(rooted_pars,main="rooted MP tree")
```

Maximum likelihood trees
make another starting tree
```{r}
starting_tree2<-nj(dist.dna(dna,model="K80"))

```

initial fitting likelihood model

```{r}
fit_init<-pml(starting_tree2,phy_dat)
fit_init
```

optimal ML tree

```{r}

fit_jc<-optim.pml(fit_init,model="JC", rearrangement = "NNI")
fit_jc
```

lets try a gtr model
```{r}
fit_gtr<-optim.pml(fit_init,model="GTR",rearrangement = "NNI")
AIC(fit_jc,fit_gtr)
```

lets plot best ML model
```{r}

mltree<-fit_gtr$tree
mltree$tip.label<-upgma_tree$tip.label
plot(mltree,main="best ML tree")
```

root the tree
```{r}

root_mltree<-root(mltree,outgroup = "Bee",resolve.root = T)
plot(root_mltree, main="best rooted ML tree")
```
 lets plot ML and MP together
 
```{r}
par(mfrow=c(1,2))
plot(rooted_pars,main="MP tree")
plot(root_mltree,main="ML tree")
```
 
 lets bootstrap
 first our nj tree
```{r}
bootnj<-boot.phylo(nj_tree,dna,function(x) nj(dist.dna(x, model="K80")),B=100)
bootnj
```
 
lets see the support values

```{r}
plot(nj_tree2,main="NJ Tree with bootstrap")
nodelabels(bootnj,cex=0.8,bg="skyblue")

```
Bootstrap ML tree
```{r}
boot_ml<-bootstrap.pml(fit_gtr,bs=100,optNni="TRUE",multicore = 6)
boot_ml

```
create a consensus tree
```{r}
con_tree<-midpoint(fit_gtr$tree)
con_tree$tip.label<-upgma_tree$tip.label
plot(con_tree)
```

The Maximum clade credibility tree

```{r}
mcc_tree<-maxCladeCred(boot_ml,rooted = T)
mcc_tree2<-mcc_tree;mcc_tree2$tip.label<-upgma_tree$tip.label
plot(mcc_tree2,main="Maximum Clade Credibility tree")
```

plot consensus tree with support

```{r}

plotBS(midpoint(fit_gtr$tree),boot_ml,type="phylogram",main="ML tree with Bootstrap")
```

date a tree
```{r}
ultra_tree<-chronos(con_tree, lambda=1, model="relaxed",quiet=F)
plot(ultra_tree)
add.scale.bar()
```

