# Function to plot frequency of distribution of different Wi given selac parameters
ComputeEquilibriumCodonFrequencies <- function(nuc.model="JC", base.freqs=rep(0.25, 4), nsites=1, C=4, Phi=0.5, q=4e-7, Ne=5e6, alpha=1.83, beta=0.10, gamma=0.0003990333, include.stop.codon=TRUE, numcode=1, diploid=TRUE, flee.stop.codon.rate=0.9999999) {
#To test: nuc.model="JC"; base.freqs=rep(0.25, 4); nsites=1; C=4; Phi=0.5; q=4e-7; Ne=5e6; alpha=1.83; beta=0.10; gamma=0.0003990333; include.stop.codon=TRUE; numcode=1; diploid=TRUE; flee.stop.codon.rate=0.9999999
nuc.mutation.rates <- CreateNucleotideMutationMatrix(1, model=nuc.model, base.freqs=base.freqs)
codon.index.matrix = CreateCodonMutationMatrixIndex()
codon_mutation_matrix <- matrix(nuc.mutation.rates[codon.index.matrix], dim(codon.index.matrix))
codon_mutation_matrix[is.na(codon_mutation_matrix)]=0
aa.distances <- CreateAADistanceMatrix(alpha=alpha, beta=beta, gamma=gamma, aa.properties=NULL, normalize=FALSE, poly.params=NULL, k=0)
Q_codon_array <- FastCreateAllCodonFixationProbabilityMatrices(aa.distances=aa.distances, nsites=nsites, C=C, Phi=Phi, q=q, Ne=Ne, include.stop.codon=include.stop.codon, numcode=numcode, diploid=diploid, flee.stop.codon.rate=0.9999999)
diag(codon_mutation_matrix) = 0
diag(codon_mutation_matrix) = -rowSums(codon_mutation_matrix)
.unique.aa <- c("K", "N", "T", "R", "S", "I", "M", "Q", "H", "P", "L", "E", "D", "A", "G", "V", "*", "Y", "C", "W", "F")
#Finish the Q_array codon mutation matrix multiplication here:
for(k in 1:21){
if(diploid == TRUE){
Q_codon_array[,,.unique.aa[k]] = (2 * Ne) * codon_mutation_matrix * Q_codon_array[,,.unique.aa[k]]
}else{
Q_codon_array[,,.unique.aa[k]] = Ne * codon_mutation_matrix * Q_codon_array[,,.unique.aa[k]]
}
diag(Q_codon_array[,,.unique.aa[k]]) = 0
diag(Q_codon_array[,,.unique.aa[k]]) = -rowSums(Q_codon_array[,,.unique.aa[k]])
}
starting.state <- 1
liks <- matrix(1/64, 1, 64)
#liks[,starting.state] <- 1
.nonstop.aa <- .unique.aa[-which(.unique.aa=="*")]
#eq.freq <- expm(Q_codon_array[,,optimal.aa] * 1000000, method=c("Ward77")) %*% liks[1,]
eq.freq.matrix <- matrix(nrow=64, ncol=length(.nonstop.aa))
rownames(eq.freq.matrix) <- rownames(Q_codon_array)
colnames(eq.freq.matrix) <- .nonstop.aa
for (i in sequence(length(.nonstop.aa))) {
eq.freq.matrix[,i] <- expm::expm(t(Q_codon_array[,,.nonstop.aa[i]]) * 1000000, method=c("Ward77")) %*% liks[1,]
}
return(eq.freq.matrix)
}
#' Function to plot a distribution of frequencies of codons given a 3d array of equilibrium frequency matrices
#'
#' @param eq.freq.matrices A 3d array of eq.freq.matrix returned from ComputeEquilibriumFrequencies
#' @param values The vector of labels for each matrix (i.e., different Phi values)
#' @param palette Color palette to use from RColorBrewer
#' @param lwd Line width
#' @param ... Other paramters to pass to plot()
PlotEquilbriumCodonDistribution <- function(eq.freq.matrices, values, palette="Set1", lwd=2, ...) {
colors <- RColorBrewer::brewer.pal(dim(eq.freq.matrices)[3],palette)
distributions <- list()
total.range <- c(NA)
for (i in sequence(dim(eq.freq.matrices)[3])) {
sorted <- apply(eq.freq.matrices[,,i], 2, sort, decreasing=TRUE)
distributions[[i]] <- apply(sorted, 1, mean)
total.range <- range(c(distributions[i], total.range), na.rm=TRUE)
}
plot(x=c(1,64), y=total.range, type="n", bty="n", xlab="index", ylab="Average frequency", ...)
for (i in sequence(dim(eq.freq.matrices)[3])) {
lines(sequence(length(distributions[[i]])), distributions[[i]], lwd=lwd, col=colors[i])
}
legend(x="topright", legend=values, fill=colors)
}
ComputeEquilibriumAAFitness <- function(nuc.model="JC", base.freqs=rep(0.25, 4), nsites=1, C=4, Phi=0.5, q=4e-7, Ne=5e6, alpha=1.83, beta=0.10, gamma=0.0003990333, include.stop.codon=TRUE, numcode=1, diploid=TRUE, flee.stop.codon.rate=0.9999999) {
eq.freq.matrix <- ComputeEquilibriumCodonFrequencies(nuc.model, base.freqs, nsites, C, Phi, q, Ne, alpha, beta, gamma, include.stop.codon, numcode, diploid, flee.stop.codon.rate)
codon.fitness.matrix <- eq.freq.matrix*NA
aa.distances <- CreateAADistanceMatrix(alpha=alpha, beta=beta, gamma=gamma, aa.properties=NULL, normalize=FALSE, poly.params=NULL, k=0)
aa.names <- unique(sapply(rownames(eq.freq.matrix), TranslateCodon, numcode=numcode))
aa.fitnesses <- matrix(nrow=length(aa.names), ncol=dim(eq.freq.matrix)[2])
rownames(aa.fitnesses) <- aa.names
colnames(aa.fitnesses) <- colnames(eq.freq.matrix)
for (col.index in sequence(dim(codon.fitness.matrix)[2])) {
for (row.index in sequence(dim(codon.fitness.matrix)[1])) {
codon.fitness.matrix[row.index, col.index] <- GetFitness(focal.protein=TranslateCodon(rownames(codon.fitness.matrix)[row.index], numcode), optimal.protein=colnames(codon.fitness.matrix)[col.index], aa.distances, nsites=nsites, C=1, Phi=Phi, q=q)
aa.fitnesses[TranslateCodon(rownames(codon.fitness.matrix)[row.index], numcode), col.index] <- codon.fitness.matrix[row.index, col.index]
}
}
return(list(aa.fitness.matrix=aa.fitnesses, codon.fitnesses=codon.fitness.matrix, equilibrium.codon.frequency = eq.freq.matrix))
}
# Computes the distribution of fitness differences between new and original mutations as well as the frequency with which those are attempted
# Returns a list; most of the return objects are the same as from ComputeEquilibriumAAFitness() but the new things are
# codon.mutation.matrix: instantaneous rate matrix for codon mutations
# codon.relative.rate.matrix: the above, but scaled so that each row sums to 1
# delta.fitness.array: 3d array: dimensions are starting codon, optimal aa, and finishing codon; entries are new codon fitness - original codon fitness
# frequency.array: 3d array: dimensions are starting codon, optimal aa, and finishing codon; entries are frequencies that that mutation is attempted given the optimal aa (column)
#
ComputeMutationFitnesses <- function(nuc.model="JC", base.freqs=rep(0.25, 4), nsites=1, C=4, Phi=0.5, q=4e-7, Ne=5e6, alpha=1.83, beta=0.10, gamma=0.0003990333, include.stop.codon=TRUE, numcode=1, diploid=TRUE, flee.stop.codon.rate=0.9999999) {
equilibrium.values <- ComputeEquilibriumAAFitness(nuc.model=nuc.model, base.freqs=base.freqs, nsites=nsites, C=C, Phi=Phi, q=q, Ne=Ne, alpha=alpha, beta=beta, gamma=gamma, include.stop.codon=include.stop.codon, numcode=numcode, diploid=diploid, flee.stop.codon.rate=flee.stop.codon.rate)
codon.fitness.matrix <- equilibrium.values$codon.fitnesses
equilibrium.codon.frequency <- equilibrium.values$equilibrium.codon.frequency
nuc.mutation.rates <- CreateNucleotideMutationMatrix(1, model=nuc.model, base.freqs=base.freqs)
codon.index.matrix = CreateCodonMutationMatrixIndex()
codon_mutation_matrix <- matrix(nuc.mutation.rates[codon.index.matrix], dim(codon.index.matrix))
codon_mutation_matrix[is.na(codon_mutation_matrix)]=0
diag(codon_mutation_matrix) = 0
scaling.factors <- rowSums(codon_mutation_matrix)
diag(codon_mutation_matrix) = -rowSums(codon_mutation_matrix)
codon.relative.rate.matrix <- codon_mutation_matrix
diag(codon.relative.rate.matrix) <- 0
codon.relative.rate.matrix <- codon.relative.rate.matrix / scaling.factors #so we know the freq of each move
delta.fitness.array <- array(dim=c(64,20,64)) #row is which codon coming from, col is which aa is optimal, depth is which codon going to
dimnames(delta.fitness.array) <- list(rownames(equilibrium.codon.frequency), colnames(equilibrium.codon.frequency), rownames(equilibrium.codon.frequency))
frequency.array <- delta.fitness.array
for (optimal.aa.index in sequence(20)) {
for (from.codon.index in sequence(64)) {
for (to.codon.index in sequence(64)) {
frequency.array[from.codon.index, optimal.aa.index, to.codon.index] <- codon.relative.rate.matrix[from.codon.index, to.codon.index] * equilibrium.codon.frequency[from.codon.index, optimal.aa.index]
delta.fitness.array[from.codon.index, optimal.aa.index, to.codon.index] <- codon.fitness.matrix[to.codon.index, optimal.aa.index] - codon.fitness.matrix[from.codon.index, optimal.aa.index] #new codon - old codon
}
}
}
return(list(aa.fitness.matrix=equilibrium.values$aa.fitnesses, codon.fitnesses=equilibrium.values$codon.fitnesses, equilibrium.codon.frequency = equilibrium.values$equilibrium.codon.frequency, codon.mutation.matrix=codon_mutation_matrix, codon.relative.rate.matrix=codon.relative.rate.matrix, delta.fitness.array=delta.fitness.array, frequency.array=frequency.array))
}
#' Function to plot a distribution of fitnesses W or selection coefficients S for a given optimal aa and other terms.
#'
#' @param aa.fitness.matrices, A 3d array of aa.fitness.matrix returned from ComputeEquilibriumAAFitness (first element in return)
#' @param values The vector of labels for each matrix (i.e., different Phi values)
#' @param optimal.aa Single letter code for the optimal aa. If NULL, integrates across aa.
#' @param palette Color palette to use from RColorBrewer
#' @param lwd Line width
#' @param include.stop.codon Include stop codons
#' @param type If "histogram", do a histogram plot; if "density", do a density plot
#' @param fitness If TRUE, plot fitness W; if FALSE, plot selection coefficient S (= W- 1)
#' @param scale.x.axis.by.Ne if TRUE, x axis is transformed from S to S*Ne; if FALSE no scaling is done
#' @param legend.title Sets the title of the figure legend.
#' @param Ne used to scale x axis when scale.x.axis.by.Ne is TRUE
#' @param ... Other paramters to pass to plot()
PlotPerAAFitness <- function(aa.fitness.matrices, values, optimal.aa=NULL, palette="Set1", lwd=2, include.stop.codon=FALSE, type="histogram", fitness=TRUE, scale.x.axis.by.Ne=FALSE, legend.title=NULL, Ne=10^6, ...) {
colors <- RColorBrewer::brewer.pal(dim(aa.fitness.matrices)[3],palette)
distributions <- list()
y.range <- c()
##scale values based on 'fitness' and set legend label if necessary
if(fitness) {
if(is.null(legend.title)) {
legend.title="W"
}
}else{
aa.fitness.matrices <- aa.fitness.matrices -1 #S_i = W_i - W_*
if(is.null(legend.title)) {
legend.title="s"
}
if(scale.x.axis.by.Ne){
aa.fitness.matrices <- aa.fitness.matrices * Ne
if(is.null(legend.title)) {
legend.title="s Ne"
}
}
}
x.range <- c(NA)
for (i in sequence(dim(aa.fitness.matrices)[3])) {
distribution <- NA
local.matrix <- aa.fitness.matrices[,,i]
if(!include.stop.codon) {
local.matrix <- local.matrix[which(rownames(local.matrix) != "*"),]
}
input.values <- NA
if (is.null(optimal.aa)) {
input.values <- c(local.matrix)
} else {
input.values <- c(local.matrix[,optimal.aa])
}
distribution <- list()
if(type=="density") {
distribution <- stats::density(input.values, to=1)
distribution$y <- distribution$y / sum(distribution$y)
distributions[[i]] <- distribution
}
if(type=="histogram") {
# table.of.dist <- table(input.values)
distribution <- graphics::hist(input.values, plot=FALSE)
distribution$x <- distribution$mids
distribution$y <- distribution$counts / length(input.values)
if(length(distribution$x)==1) {
distribution$x <- median(input.values)
distribution$y <- 1
}
# distribution$x <- as.numeric(names(table.of.dist))
# distribution$y <- unname(table.of.dist)/length(input.values)
distributions[[i]] <- distribution
}
y.range <- range(c(y.range, distribution$y), na.rm=TRUE)
x.range <- range(c(x.range, distribution$x), na.rm=TRUE)
}
plot(x=x.range, y=y.range, type="n", bty="n", xlab=ifelse(fitness, "W", ifelse(scale.x.axis.by.Ne, "s Ne", "s")), ylab="Frequency", ...)
if(type=="histogram") {
for (i in sequence(length(distributions))) {
points(distributions[[i]]$x, distributions[[i]]$y, col=colors[i], pch=20)
for (j in sequence(length(distributions[[i]]$x))) {
lines(rep(distributions[[i]]$x[j],2), c(0, distributions[[i]]$y[j]), col=add.alpha(colors[i],0.2), lwd=lwd)
}
}
} else {
lines(distributions[[i]]$x, distributions[[i]]$y, col=colors[i], lwd=lwd)
}
legend(x="topleft", legend=values, fill=colors)
}
# from http://www.magesblog.com/2013/04/how-to-change-alpha-value-of-colours-in.html
add.alpha <- function(col, alpha=1){
if(missing(col))
stop("Please provide a vector of colours.")
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x[1], x[2], x[3], alpha=alpha))
}
#' Function to plot a distribution of fitnesses based on codon equilibrium freqs
#'
#' @param codon.fitnesses.matrices A 3d array of aa.fitness.matrix returned from ComputeEquilibriumAAFitness (first element in return)
#' @param codon.eq.matrices A 3d array of codon equilibrium frequencies
#' @param values The vector of labels for each matrix (i.e., different Phi values)
#' @param optimal.aa Single letter code for the optimal aa. If NULL, integrates across aa.
#' @param palette Color palette to use from RColorBrewer
#' @param lwd Line width
#' @param include.stop.codon Include stop codons
#' @param type If "histogram", do a histogram plot; if "density", do a density plot
#' @param fitness If TRUE, plot W; if FALSE, plot S (= 1 - W)
#' @param numcode The genetic code
#' @param ... Other paramters to pass to plot()
PlotExpectedFitness <- function(codon.fitnesses.matrices, codon.eq.matrices, values, optimal.aa=NULL, palette="Set1", lwd=2, include.stop.codon=FALSE, type="histogram", fitness=TRUE, numcode=1, ...) {
colors <- RColorBrewer::brewer.pal(dim(codon.fitnesses.matrices)[3],palette)
distributions <- list()
y.range <- c()
if(!fitness) {
codon.fitnesses.matrices <- 1 - codon.fitnesses.matrices
}
x.range <- c(NA)
for (i in sequence(dim(codon.fitnesses.matrices)[3])) {
distribution <- NA
local.codon.fitness.matrix <- codon.fitnesses.matrices[,,i]
local.codon.equilibrium.matrix <- codon.eq.matrices[,,i]
if(!include.stop.codon) {
aa.vector <- unname(sapply(rownames(local.codon.fitness.matrix), TranslateCodon, numcode=numcode))
local.codon.fitness.matrix <- local.codon.fitness.matrix[which(aa.vector != "*"),]
local.codon.equilibrium.matrix <- local.codon.equilibrium.matrix[which(aa.vector != "*"),]
}
input.values <- NA
if (is.null(optimal.aa)) {
input.values <- rep(c(local.codon.fitness.matrix), round(10000*c(local.codon.equilibrium.matrix)))
} else {
input.values <- rep(c(local.codon.fitness.matrix[,optimal.aa]), round(10000*c(local.codon.equilibrium.matrix[,optimal.aa])))
}
distribution <- list()
if(type=="density") {
distribution <- stats::density(input.values, to=1)
distribution$y <- distribution$y / sum(distribution$y)
distributions[[i]] <- distribution
}
if(type=="histogram") {
# table.of.dist <- table(input.values)
distribution <- graphics::hist(input.values, plot=FALSE)
distribution$x <- distribution$mids
distribution$y <- distribution$counts / length(input.values)
if(length(distribution$x)==1) {
distribution$x <- median(input.values)
distribution$y <- 1
}
# distribution$x <- as.numeric(names(table.of.dist))
# distribution$y <- unname(table.of.dist)/length(input.values)
distributions[[i]] <- distribution
}
y.range <- range(c(y.range, distribution$y), na.rm=TRUE)
x.range <- range(c(x.range, distribution$x), na.rm=TRUE)
}
plot(x=x.range, y=y.range, type="n", bty="n", xlab=ifelse(fitness, "W", "S"), ylab="Frequency", ...)
if(type=="histogram") {
for (i in sequence(length(distributions))) {
points(distributions[[i]]$x, distributions[[i]]$y, col=colors[i], pch=20)
for (j in sequence(length(distributions[[i]]$x))) {
lines(rep(distributions[[i]]$x[j],2), c(0, distributions[[i]]$y[j]), col=add.alpha(colors[i],0.2), lwd=lwd)
}
}
} else {
lines(distributions[[i]]$x, distributions[[i]]$y, col=colors[i], lwd=lwd)
}
legend(x="topleft", legend=values, fill=colors)
}
# Generates a line for mutation fitness spectra
LineMutationFitnessSpectra <- function(mutation.fitness.object, optimal.aa=NULL) {
delta.fitness.array <- mutation.fitness.object$delta.fitness.array
frequency.array <- mutation.fitness.object$frequency.array
if(!is.null(optimal.aa)) {
delta.fitness.array <- delta.fitness.array[,which(dimnames(delta.fitness.array)[[2]]==optimal.aa),]
frequency.array <- frequency.array[,which(dimnames(frequency.array)[[2]]==optimal.aa),]
}
result <- stats::density(x=c(delta.fitness.array), weights=c(frequency.array)/sum(frequency.array))
return(result)
}
#' Plot fitness of mutations, weighted by frequency of those mutations
#'
#' @param mutation.fitness.object.list List that contains multiple objects from ComputeMutationFitnesses() calls
#' @param values The vector of labels for each matrix (i.e., different Phi values)
#' @param optimal.aa Single letter code for the optimal aa. If NULL, integrates across aa.
#' @param palette Color palette to use from RColorBrewer
#' @param lwd Line width
#' @param ... other arguments to pass to plot()
#'
PlotMutationFitnessSpectra <- function(mutation.fitness.object.list, values, optimal.aa=NULL, palette="Set1", lwd=2, ...) {
colors <- add.alpha(RColorBrewer::brewer.pal(length(mutation.fitness.object.list),palette),0.5)
results.to.plot <- lapply(mutation.fitness.object.list, LineMutationFitnessSpectra, optimal.aa=optimal.aa)
x.range <- c(NA)
y.range <- c(NA)
for (i in sequence(length(results.to.plot))) {
x.range <- range(results.to.plot[[i]]$x, na.rm=TRUE)
y.range <- range(results.to.plot[[i]]$y, na.rm=TRUE)
}
plot(x=x.range, y=y.range, bty="n", xlab="W", ylab="density", type="n",...)
for (i in sequence(length(results.to.plot))) {
lines(results.to.plot[[i]], lwd=lwd,col=colors[i])
}
legend(x="topleft", legend=values, fill=colors)
}
#' Function to plot info by site in a gene
#'
#' @param all.info The output of GetGeneSiteInfo
#' @param aa.properties The aa.properties you want to use; if NULL, uses Grantham
#' @param mean.width Sliding window width
PlotGeneSiteInfo <- function(all.info, aa.properties=NULL, mean.width=10) {
if(is.null(aa.properties)) {
aa.properties <- structure(c(0, 2.75, 1.38, 0.92, 0, 0.74, 0.58, 0, 0.33, 0, 0,
1.33, 0.39, 0.89, 0.65, 1.42, 0.71, 0, 0.13, 0.2, 8.1, 5.5, 13,
12.3, 5.2, 9, 10.4, 5.2, 11.3, 4.9, 5.7, 11.6, 8, 10.5, 10.5,
9.2, 8.6, 5.9, 5.4, 6.2, 31, 55, 54, 83, 132, 3, 96, 111, 119,
111, 105, 56, 32.5, 85, 124, 32, 61, 84, 170, 136), .Dim = c(20L,
3L), .Dimnames = list(c("A", "C", "D", "E", "F", "G",
"H", "I", "K", "L", "M", "N", "P", "Q", "R",
"S", "T", "V", "W", "Y"), c("c", "p", "v"))) #properties from Grantham paper
}
aa.properties.reordered <- aa.properties[order(match(rownames(aa.properties), .unique.aa)),]
info.by.site <- all.info$site.aa.information
dimnames(info.by.site)[[1]] <- .unique.aa
info.by.site <- info.by.site[which(.unique.aa!="*"),,]
AIC.site.information <- -2*info.by.site
get.delta <- function(x) {
return(x-min(x))
}
weights.integrating.phi <- matrix(nrow=dim(AIC.site.information)[1], ncol=dim(AIC.site.information)[2])
for (site.number in sequence(dim(AIC.site.information)[2])) { #fine, I give up on the apply across three dimensions. BCO.
local.matrix <- exp(-0.5*get.delta(AIC.site.information[,site.number,]))
local.matrix <- local.matrix / sum(local.matrix)
weights.integrating.phi[,site.number] <- rowSums(local.matrix)
}
rownames(weights.integrating.phi) <- dimnames(AIC.site.information)[[1]]
reverse.weighted.mean <- function(w, x) {
return(stats::weighted.mean(x, w))
}
average.c <- apply(weights.integrating.phi, 2, reverse.weighted.mean, x=aa.properties.reordered[,"c"])
average.p <- apply(weights.integrating.phi, 2, reverse.weighted.mean, x=aa.properties.reordered[,"p"])
average.v <- apply(weights.integrating.phi, 2, reverse.weighted.mean, x=aa.properties.reordered[,"v"])
sliding.c <- zoo::rollmean(average.c, k=mean.width)
sliding.p <- zoo::rollmean(average.p, k=mean.width)
sliding.v <- zoo::rollmean(average.v, k=mean.width)
#average.phi <- apply()
par(mfcol=c(1,4))
plot(x=sequence(length(average.c)), y=average.c, main="Composition", xlab="Site", pch=20,ylab="", bty="n", col=rgb(0,0,0,.5))
lines(x=sequence(length(sliding.c)), y=sliding.c, lwd=2)
plot(x=sequence(length(average.p)), y=average.p, main="Polarity", xlab="Site", pch=20, ylab="", bty="n", col=rgb(0,0,0,.5))
lines(x=sequence(length(sliding.p)), y=sliding.p, lwd=2)
plot(x=sequence(length(average.v)), y=average.v, main="Molecular volume", xlab="Site", pch=20, ylab="", bty="n", col=rgb(0,0,0,.5))
lines(x=sequence(length(sliding.v)), y=sliding.v, lwd=2)
}
ComputeMutationFitnessesUnderGammaRates <- function(nuc.model="JC", base.freqs=rep(0.25, 4), nsites=1, C=4, Phi=0.5, q=4e-7, Ne=5e6, alpha=1.83, beta=0.10, gamma=0.0003990333, include.stop.codon=TRUE, numcode=1, diploid=TRUE, flee.stop.codon.rate=0.9999999, shape.gamma=1, n.pulls=1000) {
Phi.vector <- Phi*rgamma(n.pulls, shape=shape.gamma, rate=shape.gamma)
ComputeMutationFitnessesPhiFirst <- function(Phi=0.5, nuc.model="JC", base.freqs=rep(0.25, 4), nsites=1, C=4, q=4e-7, Ne=5e6, alpha=1.83, beta=0.10, gamma=0.0003990333, include.stop.codon=TRUE, numcode=1, diploid=TRUE, flee.stop.codon.rate=0.9999999) {
return(ComputeMutationFitnesses(nuc.model=nuc.model, base.freqs=base.freqs, nsites=nsites, C=C, Phi=Phi, q=q, Ne=Ne, alpha=alpha, beta=beta, gamma=gamma, include.stop.codon=include.stop.codon, numcode=numcode, diploid=diploid, flee.stop.codon.rate=flee.stop.codon.rate))
}
results <- lapply(Phi.vector, ComputeMutationFitnessesPhiFirst, nuc.model=nuc.model, base.freqs=base.freqs, nsites=nsites, C=C, q=q, Ne=Ne, alpha=alpha, beta=beta, gamma=gamma, include.stop.codon=include.stop.codon, numcode=numcode, diploid=diploid, flee.stop.codon.rate=flee.stop.codon.rate)
return(results)
}
# PlotCDFOfMutations <- function(mutation.fitness.object.list, values, optimal.aa=NULL, palette="Set1", lwd=2, ...) {
# colors <- add.alpha(RColorBrewer::brewer.pal(length(mutation.fitness.object.list),palette),0.5)
# results.to.plot <- lapply(mutation.fitness.object.list, LineMutationFitnessSpectra, optimal.aa=optimal.aa)
# }
# TODO
#
# Y axis: fixation probability relative to neutral
#
# Do for different amino acids
#
# Phi*g
#
# X axis is log(W)*Ne
#
# Lines are cdfs for diff amino acids
#
# Could have theta curve also on the plot
#
# When have gamma, pull from distribution, not just the category midpoints
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.