Nothing
# 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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.