################# 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) 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)) ## 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 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) # 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) 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 num.nts names(mtDA.multSeqs.bstr) names(mtDNA.multSeqs.bstr) mtDNA.multSeqs.bstr mtDNA.multSeqs.bstr mtDNA.multSeqs.bstr mtDNA.multSeqs.bstr --help print(mtDNA.multSeqs.bstr) print(mtDNA.multSeqs.bstr, -n40) print(mtDNA.multSeqs.bstr, -n 40) print(mtDNA.multSeqs.bstr, n=20) class(mtDNA.multSeqs.bstr) print(mtDNA.multSeqs.bstr) ?print() ?print table(mtDNA.multSeqs.bstr) mtDNA.multSeqs.bstr mtDNA.multSeqs.bstr$width mtDNA.multSeqs.bstr[,1]$width mtDNA.multSeqs.bstr[1,]$width mtDNA.multSeqs.bstr[1]$width mtDNA.multSeqs.bstr[1] mtDNA.multSeqs.bstr[1]$seq mtDNA.multSeqs.bstr[1]$width mtDNA.multSeqs.bstr[1]$names mtDNA.multSeqs.bstr$names # name the Biostring sequences and compute stats names(mtDNA.multSeqs.bstr) <- Names.vec # count nucs and sequence lengths mtDNA.multSeqs.bstr$names mtDNA.multSeqs.bstr$Names 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