Finalizing multidimensional scaling

.Rhistory

@@ -1,107 +1,3 @@
-################# Max-Log-Likelihood #################
-n <- length(g.breaks.clean)
-kmin <- g.breaks.clean[1]
-alpha.ML <- 1 + n/sum(log(g.breaks.clean)/kmin)
-alpha.ML
-lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4)
-# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course
-# Degree Distribution
-# Professor: Dr. McKinney, Spring 2022
-# Noah Schrick - 1492657
-library(igraph)
-library(igraphdata)
-data(yeast)
-g <- yeast
-g.netname <- "Yeast"
-################# Set up Work #################
-g.vec <- degree(g)
-g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname,
-" Network"))
-legend("topright", c("Guess", "Poisson", "Least-Squares Fit",
-"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6",
-"#006CD1", "#E66100", "#D35FB7"))
-g.mean <- mean(g.vec)
-g.seq <- 0:max(g.vec) # x-axis
-################# Guessing Alpha #################
-alpha.guess <- 1.5
-lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1)
-################# Poisson #################
-g.pois <- dpois(g.seq, g.mean, log=F)
-lines(g.seq, g.pois, col="#006CD1", lty=2)
-################# Linear model: Least-Squares Fit #################
-g.breaks <- g.hist$breaks[-c(1,2,3)] # remove 0
-g.probs <- g.hist$density[-1] # make lengths match
-# Need to clean up probabilities that are 0
-nz.probs.mask <- g.probs!=0
-g.breaks.clean <- g.breaks[nz.probs.mask]
-g.probs.clean <- g.breaks[nz.probs.mask]
-#plot(log(g.breaks.clean), log(g.probs.clean))
-g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean))
-summary(g.fit)
-alpha.LM <- coef(g.fit)[2]
-lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3)
-################# Max-Log-Likelihood #################
-n <- length(g.breaks.clean)
-kmin <- g.breaks.clean[1]
-alpha.ML <- 1 + n/sum(log(g.breaks.clean)/kmin)
-alpha.ML
-lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4)
-# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course
-# Degree Distribution
-# Professor: Dr. McKinney, Spring 2022
-# Noah Schrick - 1492657
-library(igraph)
-library(igraphdata)
-data(yeast)
-g <- yeast
-g.netname <- "Yeast"
-################# Set up Work #################
-g.vec <- degree(g)
-g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname,
-" Network"))
-legend("topright", c("Guess", "Poisson", "Least-Squares Fit",
-"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6",
-"#006CD1", "#E66100", "#D35FB7"))
-g.mean <- mean(g.vec)
-g.seq <- 0:max(g.vec) # x-axis
-################# Guessing Alpha #################
-alpha.guess <- 1.5
-lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1)
-################# Poisson #################
-g.pois <- dpois(g.seq, g.mean, log=F)
-lines(g.seq, g.pois, col="#006CD1", lty=2)
-################# Linear model: Least-Squares Fit #################
-g.breaks <- g.hist$breaks[-c(1)] # remove 0
-g.probs <- g.hist$density[-1] # make lengths match
-# Need to clean up probabilities that are 0
-nz.probs.mask <- g.probs!=0
-g.breaks.clean <- g.breaks[nz.probs.mask]
-g.probs.clean <- g.breaks[nz.probs.mask]
-#plot(log(g.breaks.clean), log(g.probs.clean))
-g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean))
-summary(g.fit)
-alpha.LM <- coef(g.fit)[2]
-lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3)
-################# Max-Log-Likelihood #################
-n <- length(g.breaks.clean)
-kmin <- g.breaks.clean[1]
-alpha.ML <- 1 + n/sum(log(g.breaks.clean)/kmin)
-alpha.ML
-lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4)
-# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course
-# Degree Distribution
-# Professor: Dr. McKinney, Spring 2022
-# Noah Schrick - 1492657
-library(igraph)
-library(igraphdata)
-data(yeast)
-g <- yeast
-g.netname <- "Yeast"
-################# Set up Work #################
-g.vec <- degree(g)
-g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname,
-" Network"))
-legend("topright", c("Guess", "Poisson", "Least-Squares Fit",
 "Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6",
 "#006CD1", "#E66100", "#D35FB7"))
 g.mean <- mean(g.vec)
@@ -406,6 +302,13 @@ plot(log(g.breaks.clean), log(g.probs.clean))
 g.breaks.clean <- g.breaks[nz.probs.mask]
 g.probs.clean <- g.probs[nz.probs.mask]
 plot(log(g.breaks.clean), log(g.probs.clean))
+BiocManager::install()
+BiocManager::install()
+BiocManager::install("stringi")
+# 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
@@ -420,17 +323,15 @@ mtDNA.MultiSeqs.list<-read.GenBank(c("AF011222","AF254446","X90314","AF089820",
 mtDNA.Species<-attr(mtDNA.MultiSeqs.list,"species")
 # use species as name instead of genbank id
 names(mtDNA.MultiSeqs.list)<-mtDNA.Species
-names(mtDNA.MultiSeqs.list)
 # 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="")
-names(mtDNA.MultiSeqs.list)
+length(mtDNA.MultiSeqs.list$Human)
 # look at one of the sequences using $
 mtDNA.MultiSeqs.list$Human
-length(mtDNA.MultiSeqs.list$Human)
 ## Convert to Biostrings object for the sequences
 if (!require("BiocManager")) install.packages("BiocManager")
 library(BiocManager)
@@ -450,63 +351,162 @@ 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
-num.nts
+# Obtain name and length of species with longest sequence
-names(mtDA.multSeqs.bstr)
+nlengthnames <- cbind(mtDNA.lengths, Names.vec)
-names(mtDNA.multSeqs.bstr)
+idx <- which.max(nlengthnames[,1])
-mtDNA.multSeqs.bstr
+nlengthnames[idx,]
-mtDNA.multSeqs.bstr
+#### Part B: Multiple Sequence Alignment
-mtDNA.multSeqs.bstr
+if (!require("BiocManager")) install.packages("BiocManager")
-mtDNA.multSeqs.bstr --help
+library(BiocManager)
-print(mtDNA.multSeqs.bstr)
+if (!require("msa")) BiocManager::install("msa")
-print(mtDNA.multSeqs.bstr, -n40)
+library(msa)
-print(mtDNA.multSeqs.bstr, -n 40)
+msa(mtDNA.multSeqs.bstr,method="ClustalOmega")
-print(mtDNA.multSeqs.bstr, n=20)
+msa(mtDNA.multSeqs.bstr,method="Muscle")
-class(mtDNA.multSeqs.bstr)
+# Lab 10 for the University of Tulsa's CS-6643 Bioinformatics Course
-print(mtDNA.multSeqs.bstr)
+# Phylogenetic Analysis
-?print()
+# Professor: Dr. McKinney, Fall 2022
-?print
+# Noah L. Schrick - 1492657
-table(mtDNA.multSeqs.bstr)
+## Set Working Directory to file directory - RStudio approach
-mtDNA.multSeqs.bstr
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
-mtDNA.multSeqs.bstr$width
+#### Part A: GenBank sequences and a multiple fasta file
-mtDNA.multSeqs.bstr[,1]$width
+if (!require("ape")) install.packages("ape")
-mtDNA.multSeqs.bstr[1,]$width
+library(ape) # needed for read.GenBank
-mtDNA.multSeqs.bstr[1]$width
+# fetch the mtDNA sequences
-mtDNA.multSeqs.bstr[1]
+mtDNA.MultiSeqs.list<-read.GenBank(c("AF011222","AF254446","X90314","AF089820",
-mtDNA.multSeqs.bstr[1]$seq
+"AF176766","AF451972", "AY079510",
-mtDNA.multSeqs.bstr[1]$width
+"AF050738","AF176722","AF315498",
-mtDNA.multSeqs.bstr[1]$names
+"AF176731","AF451964"), as.character=TRUE)
-mtDNA.multSeqs.bstr$names
+# 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
-mtDNA.multSeqs.bstr$names
+num.nts <- alphabetFrequency(mtDNA.multSeqs.bstr)[,1:4]
-mtDNA.multSeqs.bstr$Names
+mtDNA.lengths <- rowSums(num.nts)
-mtDNA.lengths
-table(mtDNA.lengths, Names.vec)
-cbind(mtDNA.lengths, Names.vec)
-sort(cbind(mtDNA.lengths, Names.vec))
-cbind(mtDNA.lengths, Names.vec)
-cbind(mtDNA.lengths, Names.vec)
-table(cbind(mtDNA.lengths, Names.vec))
-rbind(cbind(mtDNA.lengths, Names.vec))
-sort(rbind(cbind(mtDNA.lengths, Names.vec)))
-rbind(mtDNA.lengths, Names.vec)
-cbind(mtDNA.lengths, Names.vec)
-max(cbind(mtDNA.lengths, Names.vec))
-max(cbind(mtDNA.lengths, Names.vec))[,1]
-max(cbind(mtDNA.lengths, Names.vec)[,1])
-max(cbind(mtDNA.lengths, Names.vec)[1,])
-max(cbind(mtDNA.lengths, Names.vec)[,1])
-max(cbind(mtDNA.lengths, Names.vec))
-cbind(mtDNA.lengths, Names.vec)
-nlengthnames <- cbind(mtDNA.lengths, Names.vec)
-max(nlengthnames[,1])
-nlengthnames <- cbind(mtDNA.lengths, Names.vec)
-nlengthnames[which.max(nlengthnames[,1])]
-idx <- which.max(nlengthnames[,1])
-idx
-nlengthnames[idx, idx]
-nlengthnames[idx]
-nlengthnames
-nlengthnames[idx,]
 proportion.nts <- num.nts/mtDNA.lengths
-proportion.nts
+# 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
+align.stats.df[1]
+align.stats.df[1,]
+#### 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)
+min(mtDNA.dist.mat)
+## 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(myDist))
+#### 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=1.5)
+?text()
+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)
+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)
+plot3d(x,y,z)
+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)
+q
+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)
+plot3d(x,y,z)
+text3d(x=x,y=y,z=z,texts=rownames(locs),cex=1.5)

Schrick-Noah_CS-6643_Lab-10.R

@@ -93,3 +93,57 @@ write.table(align.stats.df,file="alignstats.tab",sep = "\t", row.names=FALSE, qu
 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)
+