# Lab 10 for the University of Tulsa's CS-6643 Bioinformatics Course
# Phylogenetic Analysis
# Professor: Dr. McKinney, Fall 2022
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#### Part A: GenBank sequences and a multiple fasta file
if (!require("ape")) install.packages("ape")
library(ape) # needed for read.GenBank

# fetch the mtDNA sequences
mtDNA.MultiSeqs.list<-read.GenBank(c("AF011222","AF254446","X90314","AF089820",
                                     "AF176766","AF451972", "AY079510",
                                     "AF050738","AF176722","AF315498", 
                                     "AF176731","AF451964"), as.character=TRUE)
# look at species names
mtDNA.Species<-attr(mtDNA.MultiSeqs.list,"species")
# use species as name instead of genbank id
names(mtDNA.MultiSeqs.list)<-mtDNA.Species
# need to fix some names
names(mtDNA.MultiSeqs.list)[1] <- paste("German_Neanderthal",sep="")
names(mtDNA.MultiSeqs.list)[2] <- paste("Russian_Neanderthal",sep="")
names(mtDNA.MultiSeqs.list)[3] <- paste("Human")
names(mtDNA.MultiSeqs.list)[6] <- paste("Puti_Orangutan",sep="")
names(mtDNA.MultiSeqs.list)[12] <- paste("Jari_Orangutan",sep="")

length(mtDNA.MultiSeqs.list$Human)

# look at one of the sequences using $
mtDNA.MultiSeqs.list$Human

## Convert to Biostrings object for the sequences
if (!require("BiocManager")) install.packages("BiocManager")
library(BiocManager)
if (!require("Biostrings")) BiocManager::install("Biostrings")
library(Biostrings)
# loop through the list to create vector of strings for Biostrings input
Names.vec <- c() # initialize speices names string vector
Seqs.vec <- c()  # initialize sequence string vector
for (mtDNA.name in names(mtDNA.MultiSeqs.list))
{
  Names.vec <- c(Names.vec,mtDNA.name) # concatenate vector
  Seqs.vec <-c(Seqs.vec,paste(mtDNA.MultiSeqs.list[[mtDNA.name]],collapse=""))
}
mtDNA.multSeqs.bstr <- DNAStringSet(Seqs.vec) # convert to Biostring

# name the Biostring sequences and compute stats
names(mtDNA.multSeqs.bstr) <- Names.vec # count nucs and sequence lengths
num.nts <- alphabetFrequency(mtDNA.multSeqs.bstr)[,1:4]
mtDNA.lengths <- rowSums(num.nts)
proportion.nts <- num.nts/mtDNA.lengths

# Obtain name and length of species with longest sequence
nlengthnames <- cbind(mtDNA.lengths, Names.vec)
idx <- which.max(nlengthnames[,1])
nlengthnames[idx,]


#### Part B: Multiple Sequence Alignment
if (!require("BiocManager")) install.packages("BiocManager")
library(BiocManager)
if (!require("msa")) BiocManager::install("msa")
library(msa)
mtDNA.msa <- msa(mtDNA.multSeqs.bstr,method="ClustalOmega")
msaPrettyPrint(mtDNA.msa, file="mtDNA.pdf", output="pdf", showNames="left",
               showLogo="none", askForOverwrite=FALSE, verbose=TRUE )
## loop to make results data frame
num_seqs <- length(Names.vec)
# initialize data frame
align.stats.df <- data.frame(species=rep(NA,num_seqs), seqlen=rep(0,num_seqs), 
                             numgaps=rep(0,num_seqs), nt_a=rep(NA,num_seqs), 
                             nt_c=rep(NA,num_seqs), nt_g=rep(NA,num_seqs), 
                             nt_t=rep(NA,num_seqs))
# DNAbin type required for dist.dna and helpful for other calculations
mtDNA.msa.DNAbin <- as.DNAbin(mtDNA.msa) 
for (i in 1:num_seqs){
  seq_name <- Names.vec[i]
  seq.vec <- as.character(mtDNA.msa.DNAbin[i,])
  num.gaps <- sum(seq.vec=="-")
  prop.nt.i <- proportion.nts[i,]
  align.stats.df[i,] <- c(seq_name, mtDNA.lengths[i], num.gaps,
                          round(prop.nt.i[1],digits=2), round(prop.nt.i[2],digits=2), 
                          round(prop.nt.i[3],digits=2), round(prop.nt.i[4],digits=2))
}

# write to file
write.table(align.stats.df,file="alignstats.tab",sep = "\t", row.names=FALSE, quote=FALSE)

# you can use $ operator to grab a named column from a data.frame 
# similar to grabbing a named variable from a list 
align.stats.df$species
align.stats.df$nt_a             # strings by default
as.numeric(align.stats.df$nt_a) # convert to numeric


#### Part C: DNA distance matrices and phylogenetic trees
# Compute Distances
mtDNA.dist <- dist.dna(mtDNA.msa.DNAbin,model="K80")
# manually find closest species
mtDNA.dist.mat <-as.matrix(mtDNA.dist)
diag(mtDNA.dist.mat)<-1  # force diagonal to be 1, not 0
which(mtDNA.dist.mat == min(mtDNA.dist.mat), arr.ind = TRUE)

## Make tree from distance matrix
hc<- hclust(as.dist(mtDNA.dist.mat)) # transform to dist object first
plot(hc,xlab="species",ylab="distance")  

## UPGMA
if (!require("phangorn")) install.packages("phangorn")
library(phangorn)
mtDNA.tree.nj <- NJ(mtDNA.dist) # phangorn function
plot(mtDNA.tree.nj, main="Neighbor Joining Tree (rooted) for primates")
mtDNA.tree.upgma  <- upgma(mtDNA.dist)
plot(mtDNA.tree.upgma, show.node.label = TRUE, main="UPGMA Tree for Primates")

source("msaUtils.R") # load msaConvert function into memory
mtDNA.msa.phangorn <-msaConvert(mtDNA.msa,type="phangorn::phyDat")
parsimony(mtDNA.tree.nj, mtDNA.msa.phangorn)
# bootstrap to show support for tree edges
# creates trees from bootstrap samples and checks how often
# that edge appears. Show consistency of tree edge. 
bs.trees <- bootstrap.phyDat(mtDNA.msa.phangorn, FUN=function(x)NJ(dist.dna(as.DNAbin(x),model="K80")), bs=100)
plotBS(mtDNA.tree.nj, bs.trees, "phylogram", main="Neighbor Joining")

parsimony(mtDNA.tree.upgma, mtDNA.msa.phangorn)

bs.upgma.trees <- bootstrap.phyDat(mtDNA.msa.phangorn, FUN=function(x)upgma(dist.dna(as.DNAbin(x),model="K80")), bs=100)
plotBS(mtDNA.tree.upgma, bs.upgma.trees, "phylogram", main="UPGMA")


#### Part D: Multidimensional Scaling
# 2d MDS viz
locs<-cmdscale(as.dist(mtDNA.dist))
x<-locs[,1]
y<-locs[,2]
plot(x,y,main="Multi-dimensional Scaling",xlab="MDS dimension-1",ylab="MDS dimension-2", xlim=c(-.3,.35))
text(x,y,rownames(locs),cex=0.5)
     
library(rgl)
locs<-cmdscale(as.dist(mtDNA.dist),k=3)
x<-locs[,1]
y<-locs[,2]
z<-locs[,3]
plot3d(x,y,z)
text3d(x=x,y=y,z=z,texts=rownames(locs),cex=1.5) 
play3d(spin3d(axis=c(0,1,1), rpm=3), duration=30)