513 lines
20 KiB
R
513 lines
20 KiB
R
"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.breaks <- g.hist$breaks # 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.probs[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)
|
|
alpha.LM
|
|
# 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.probs[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, lwd=5)
|
|
################# 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.probs[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, lwd=3)
|
|
################# 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.probs[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, lwd=3)
|
|
################# Poisson #################
|
|
g.pois <- dpois(g.seq, g.mean, log=F)
|
|
lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3)
|
|
################# 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.probs[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, lwd=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, lwd=3)
|
|
plot(yeast)
|
|
hist(yeast)
|
|
hist(g.vec)
|
|
g.pois
|
|
g.mean
|
|
alpha.LM
|
|
alpha.ML
|
|
degree(g)
|
|
sort(degree(g))
|
|
sort(degree(g),decreasing=FALSE)
|
|
sort(degree(g),decreasing=F)
|
|
sort(degree(g),decreasing=false)
|
|
sort(degree(g), decreasing = TRUE)
|
|
head(sort(degree(g), decreasing = TRUE))
|
|
stddev(degree(g))
|
|
sd(degree(g))
|
|
tail(sort(degree(g), decreasing = TRUE))
|
|
plot(log(g.breaks.clean), log(g.probs.clean))
|
|
# 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, lwd=3)
|
|
################# Poisson #################
|
|
g.pois <- dpois(g.seq, g.mean, log=F)
|
|
lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3)
|
|
################# 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.probs[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, lwd=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, lwd=3)
|
|
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
|
|
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)
|
|
msa(mtDNA.multSeqs.bstr,method="ClustalOmega")
|
|
msa(mtDNA.multSeqs.bstr,method="Muscle")
|
|
# 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
|
|
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)
|