inst/doc/anacoda.R

In AnaCoDa: Analysis of Codon Data under Stationarity using a Bayesian Framework

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  genome <- initializeGenomeObject(file = "genome.fasta")
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)))
#  model <- initializeModelObject(parameter = parameter, model = "ROC")
#  mcmc <- initializeMCMCObject(samples = 5000, thinning = 10, adaptive.width=50)
#  runMCMC(mcmc = mcmc, genome = genome, model = model)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  parameter <- initializeParameterObject(genome = genome, sphi = c(0.5, 2), num.mixtures = 2, geneAssignment = sample.int(2, length(genome), replace = T))

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  genome <- initializeGenomeObject(file = "genome.fasta", observed.expression.file = "synthesis_values.csv")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  # One case of observed data
#  sepsilon <- 0.1
#  # Two cases of observed data
#  sepsilon <- c(0.1, 0.5)
#  # ...
#  # Five cases of observed data
#  sepsilon <- c(0.1, 0.5, 1, 0.8, 3)
#  
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)), init.sepsilon = sepsilon)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  model <- initializeModelObject(parameter = parameter, model = "ROC",  fix.observation.noise = TRUE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  genome <- initializeGenomeObject(file = "genome.fasta")
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)))
#  model <- initializeModelObject(parameter = parameter, model = "ROC")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  mcmc <- initializeMCMCObject(samples, thinning=1, adaptive.width=100, est.expression=FALSE, est.csp=TRUE, est.hyper=TRUE, est.mix=TRUE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  mcmc <- initializeMCMCObject(samples, thinning=1, adaptive.width=100, est.expression=TRUE, est.csp=FALSE, est.hyper=TRUE, est.mix=TRUE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  mcmc <- initializeMCMCObject(samples, thinning=1, adaptive.width=100, est.expression=TRUE, est.csp=TRUE, est.hyper=FALSE, est.mix=TRUE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  mcmc <- initializeMCMCObject(samples, thinning=1, adaptive.width=100, est.expression=TRUE, est.csp=TRUE, est.hyper=TRUE, est.mix=FALSE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  genome <- initializeGenomeObject(file = "genome.fasta")
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)))
#  parameter$fixDM()
#  parameter$fixDEta()

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  genome <- initializeGenomeObject(file = "genome.fasta")
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)))
#  
#  fix_dm <- TRUE
#  fix_deta <- FALSE
#  
#  parameter$initMutationCategories(dM.file,1,fix_dm)
#  parameter$initSelectionCategories(dEta.file,1,fix_deta)
#  
#  parameter <- initializeParameterObject(genome = genome, sphi = c(1,1), num.mixtures = 2, geneAssignment = sample.int(2, length(genome), replace = T), mixture.definition = "mutationShared")
#  
#  fix_dm <- TRUE
#  fix_deta <- FALSE
#  
#  parameter$initMutationCategories(dM.file,1,fix_dm)
#  parameter$initSelectionCategories(dEta.file,2,fix_deta)
#  

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  parameter <- initializeParameterObject(genome = genome, sphi = c(0.5, 2), num.mixtures = 2
#                                         , geneAssignment = sample.int(2, length(genome), replace = T),
#                                         mixture.definition = "allUnique")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  parameter <- initializeParameterObject(genome = genome, sphi = c(0.5, 2), num.mixtures = 2
#                                         , geneAssignment = sample.int(2, length(genome), replace = T),
#                                         mixture.definition = "mutationShared")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  parameter <- initializeParameterObject(genome = genome, sphi = c(0.5, 2), num.mixtures = 2
#                                         , geneAssignment = sample.int(2, length(genome), replace = T),
#                                         mixture.definition = "selectionShared")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  #     [,1] [,2]
#  #[1,]    1    1
#  #[2,]    1    2
#  #[3,]    1    3
#  def.matrix <- matrix(c(1,1,1,1,2,3), ncol=2)
#  parameter <- initializeParameterObject(genome = genome, sphi = c(0.5, 2, 1), num.mixtures = 3,
#                                         geneAssignment = sample.int(3, length(genome), replace = T),
#                                         mixture.definition.matrix = def.matrix)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  #     [,1] [,2]
#  #[1,]    1    1
#  #[2,]    2    2
#  #[3,]    3    3
#  def.matrix <- matrix(c(1,2,3,1,2,3), ncol=2)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  #     [,1] [,2]
#  #[1,]    1    1
#  #[2,]    2    1
#  #[3,]    1    2
#  def.matrix <- matrix(c(1,2,1,1,1,2), ncol=2)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  # writing a restart file every 1000 samples
#  setRestartSettings(mcmc, "restart_file", 1000, write.multiple=TRUE)
#  # writing a restart file every 1000 samples but overwriting it every time
#  setRestartSettings(mcmc, "restart_file", 1000, write.multiple=FALSE)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  initializeParameterObject(init.with.restart.file = "restart_file.rst")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  #save objects after a run
#  runMCMC(mcmc = mcmc, genome = genome, model = model)
#  writeParameterObject(parameter = parameter, file = "parameter_out.Rda")
#  writeMCMCObject(mcmc = mcmc, file = "mcmc_out.Rda")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  #save objects after a run
#  parameter <- loadParameterObject(file = "parameter_out.Rda")
#  mcmc <- loadMCMCObject(file = "mcmc_out.Rda")

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  parameter <- initializeParameterObject(genome = genome, sphi = 1, num.mixtures = 1, geneAssignment = rep(1, length(genome)),model="FONSE",init.initiation.cost=a1)
#  
#  model <- initializeModelObject(parameter,"FONSE")
#  

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  
#  trace <- parameter$getTraceObject()
#  
#  plot(trace,what="InitiationCost")
#  
#  trace_a1 <- trace$getInitiationCostTrace()
#  mean_a1 <- mean(trace_a1)
#  sd_a1 <- sd(trace_a1)
#  ci_a1 <- quantile(trace_a1,probs = c(0.025,0.975))
#  

## ----echo=TRUE,eval=FALSE-----------------------------------------------------
#  
#  genome.pa <- initializeGenomeObject("rfp.csv",fasta=FALSE,positional=FALSE)
#  genome.panse <- initializeGenomeObject("rfp.csv",fasta=FALSE,positional=TRUE)
#  

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  csp_mat <- getCSPEstimates(parameter = parameter, CSP="Mutation", mixture = 1, samples = 1000)
#  head(csp_mat)
#  #  AA Codon  Posterior     0.025%     0.975%
#  #1  A   GCA -0.2435340 -0.2720696 -0.2165220
#  #2  A   GCC  0.4235546  0.4049132  0.4420680
#  #3  A   GCG  0.7004484  0.6648690  0.7351707
#  #4  C   TGC  0.2016298  0.1679025  0.2387024
#  #5  D   GAC  0.5775052  0.5618199  0.5936979
#  #6  E   GAA -0.4524295 -0.4688044 -0.4356677
#  
#  getCSPEstimates(parameter = parameter, filename = "mutation.csv", CSP="Mutation", mixture = 1, samples = 1000)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  phi_mat <- getExpressionEstimates(parameter = parameter,
#                                    gene.index = 1:length(genome),
#                                    samples = 1000)
#  head(phi_mat)
#  #           Mean  Mean.log10      Std.Dev   log10.Std.Dev      0.025     0.975 log10.0.025 log10.0.975
#  #[1,] 0.2729446 -0.6188447 0.0001261525    2.362358e-04 0.07331819 0.5455295 -1.13478830 -0.26319141
#  #[2,] 1.4221716  0.1498953 0.0001669425    5.194123e-05 1.09593642 1.7562065  0.03978491  0.24457557
#  #[3,] 0.7459888 -0.1512764 0.0002313539    1.529267e-04 0.31559618 1.2198282 -0.50086958  0.08629407
#  #[4,] 0.6573082 -0.2030291 0.0001935466    1.400333e-04 0.31591233 1.0699855 -0.50043989  0.02937787
#  #[5,] 1.6316901  0.2098120 0.0001846631    4.986347e-05 1.28410352 2.0035207  0.10860000  0.30179215
#  #[6,] 0.6179711 -0.2286806 0.0001744928    1.374863e-04 0.28478950 0.9683327 -0.54550116 -0.01397541

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  # sampling 100 genes at random
#  phi_mat <- getExpressionEstimates(parameter = parameter,
#                                    gene.index = sample(1:length(genome), 100),
#                                    samples = 1000)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  selection.coefficients <- getSelectionCoefficients(genome = genome,
#                                                     parameter = parameter,
#                                                     samples = 1000)
#  head(selection.coefficients)
#  #                    GCA          GCC        GCG GCT        TGC TGT GAC         GAT ...
#  #SAKL0A00132g -0.1630284 -0.008695144 -0.2097771   0 -0.1014373   0   0 -0.05092397 ...
#  #SAKL0A00154g -0.8494558 -0.045305847 -1.0930388   0 -0.5285367   0   0 -0.26533820 ...
#  #SAKL0A00176g -0.4455753 -0.023764823 -0.5733448   0 -0.2772397   0   0 -0.13918105 ...
#  #SAKL0A00198g -0.3926068 -0.020939740 -0.5051875   0 -0.2442824   0   0 -0.12263567 ...
#  #SAKL0A00220g -0.9746002 -0.051980440 -1.2540685   0 -0.6064022   0   0 -0.30442861 ...
#  #SAKL0A00242g -0.3691110 -0.019686586 -0.4749542   0 -0.2296631   0   0 -0.11529644 ...

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  cai.weights <- getCAIweights(referenceGenome = genome)
#  head(cai.weights)
#  #      GCA       GCC       GCG       GCT       TGC       TGT
#  #0.7251276 0.6282192 0.2497737 1.0000000 0.6222628 1.0000000
#  
#  nc.per.aa <- getNcAA(genome = genome)
#  head(nc.per.aa)
#  #                    A        C        D        E        F        G ...
#  #SAKL0A00132g 3.611111 1.000000 2.200000 2.142857 1.792453 4.109589 ...
#  #SAKL0A00154g 1.843866 2.500000 2.035782 1.942505 1.986595 2.752660 ...
#  #SAKL0A00176g 5.142857       NA 1.857143 1.652174 1.551724 3.122449 ...
#  #SAKL0A00198g 3.800000       NA 1.924779 1.913043 2.129032 4.136364 ...
#  #SAKL0A00220g 3.198529 1.666667 1.741573 1.756757 2.000000 1.371638 ...
#  #SAKL0A00242g 4.500000       NA 2.095890 2.000000 1.408163 3.734043 ...
#  

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  selection.coefficients <- getSelectionCoefficients(genome = genome,
#                                                     parameter = parameter,
#                                                     samples = 1000)
#  s <- exp(selection.coefficients)
#  cai.weights <- getCAIweights(referenceGenome = ref.genome)
#  
#  codon.names <- colnames(s)
#  h <- hist(s[, 1], plot = F)
#  plot(NULL, NULL, axes = F, xlim = c(0,1), ylim = range(c(0,h$counts)),
#       xlab = "s", ylab = "Frequency", main = codon.names[1], cex.lab = 1.2)
#  lines(x = h$breaks, y = c(0,h$counts), type = "S", lwd=2)
#  abline(v = cai.weights[1], lwd=2, lty=2)
#  axis(1, lwd = 3, cex.axis = 1.2)
#  axis(2, lwd = 3, cex.axis = 1.2)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  trace <- getTrace(parameter)
#  plot(x = trace, what = "Mutation", mixture = 1)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  trace <- parameter$getTraceObject()
#  plot(x = trace, what = "Expression", mixture = 1, geneIndex = 669)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  plot(mcmc, what = "LogPosterior", zoom.window = c(9000, 10000))

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  # use the last 500 samples from mixture 1 for posterior estimate.
#  plot(x = model, genome = genome, samples = 500, mixture = 1)

## ---- echo = TRUE, eval = FALSE-----------------------------------------------
#  plot(parameter, what = "Selection", samples = 500)