petrelharp / context Goto Github PK

View Code? Open in Web Editor NEW

0.0 0.0 0.0 41.92 MB

Context-dependent mutation rate inference machinery.

R 27.00% Shell 4.19% Makefile 0.44% Python 1.13% TeX 19.88% HTML 47.34% Dockerfile 0.02%

context's People

Contributors

Watchers

context's Issues

Error: Must specify a long clade name.

stoat context/json-cpg ‹master*› » ./workflow.sh 
Loading required package: optparse
Loading required package: Matrix
Loading required package: methods
Loading required package: expm

Attaching package: ‘expm’

The following object is masked from ‘package:Matrix’:

    expm

Loading required package: stringdist
Loading required package: parallel
Loading required package: jsonlite

Attaching package: ‘jsonlite’

The following object is masked from ‘package:utils’:

    View

Creating a generic function for ‘coef’ from package ‘stats’ in the global environment
Creating a generic function for ‘fitted’ from package ‘stats’ in the global environment
Creating a generic function for ‘residuals’ from package ‘stats’ in the global environment
Loading required package: optparse
Loading required package: Matrix

Attaching package: ‘Matrix’

The following object is masked from ‘package:IRanges’:

    expand

Loading required package: expm

Attaching package: ‘expm’

The following object is masked from ‘package:Matrix’:

    expm

Loading required package: stringdist
Loading required package: jsonlite

Attaching package: ‘jsonlite’

The following object is masked from ‘package:utils’:

    View

Creating a new generic function for ‘counts’ in the global environment
Creating a generic function for ‘coef’ from package ‘stats’ in the global environment
Creating a generic function for ‘fitted’ from package ‘stats’ in the global environment
Error: Must specify a long clade name.
Execution halted

Set up Travis for testing our R code

Using https://github.com/craigcitro/r-travis

@petrelharp I'm curious how you feel about this sort of thing. I know that you are a bit more on the side of hand-testing rather than Testing, but I think it would be nice to know when something we've done has broken the repo in some way.

My suggestion is to actually do the feature-branch thing and have those be the units of testing. This works great (though it requires a little discipline.) See e.g. https://speakerdeck.com/holman/how-github-uses-github-to-build-github.

make unit test for Rscript interface

source("testing-code/testing-ising.R") tests the underlying functions, but doesn't test the Rscript -- interface. Write a script that does a similar test, but using e.g. Rscript flex-genmat.R and Rscript fit-model.R, etcetera.

translate examples to new style and clean up unused code

would be nice to keep cpg, ising, tasep, but translate them to new style

Rscript-callable script that parses output and gives summaries

This is @petrelharp 's idea for sure, and it's a good one. I thought that we could discuss what would happen here.

Not to be all hipster-coder, but I think it would be fun to have it spit out an R-markdown document and/or its html counterpart with the inferences and MCMC diagnostics...

object 'mutrates' not found

cd nestly
mkdir -p _simple/cpg-model.json/3_1_1/
Rscript /home/matsen/re/context/sim-seq.R --seed 0 --jobid 0 -c models/cpg-model.json -t 0.1 -s 10000 -o _simple/cpg-model.json/3_1_1/sim.RData
Rscript /home/matsen/re/context/count-seq.R -i _simple/cpg-model.json/3_1_1/sim.RData -w 3 -s 1 -l 1 -o _simple/cpg-model.json/3_1_1/counts.tsv
Rscript /home/matsen/re/context/make-genmat.R -c models/cpg-model.json -w 3 -o _simple/cpg-model.json/3_1_1/genmat.RData
Rscript /home/matsen/re/context/fit-model.R --seed 0 --jobid 0 -i _simple/cpg-model.json/3_1_1/counts.tsv -m _simple/cpg-model.json/3_1_1/genmat.RData -o _simple/cpg-model.json/3_1_1/fit.RData

....

Read 1 item
Error in unlist(c(mutrates * opt$tlen, selcoef, fixfn.params)) : 
  object 'mutrates' not found
Calls: with -> with.default -> eval -> eval -> unlist
Error in " /usr/lib/R/bin/exec/R --slave --no-restore --file=/home/matsen/re/context/fit-model.R --args --seed 0 --jobid 0 -i _simple/cpg-model.json/3_1_1/counts.tsv -m _simple/cpg-model.json/3_1_1/genmat.RData -o _simple/cpg-model.json/3_1_1/fit.RData ": dumping frames.

make an R package

http://r-pkgs.had.co.nz/intro.html

write script to get cheap estimates of initial values and scaling for parameters

Since bad scaling, or initial values, can prevent convergence, these have to be specified in the config file now. Take code commented out of fit-model.R, and add some code to explore the likelihood surface (around the truth in simulated data?) to get cheap estimates of these things.

add in tree information to `context` class

Strictify terminology

Here's the current state of the notation. I've signified terms that appear in the tex with a * and those that appear in the code (including model-desc.md) with a ^.

pattern*^: a general word for a given subsequence of some length
k-mer: a pattern appearing in an input sequence
mutation pattern or mutpat^: a pair (a,b) of patterns where a mutates to b (do a and b have to be the same length?)
motif^: a pattern that appears in a model as an a in an (a,b) mutpat for some b

put all testing scripts in testing-code

and make them runnable: either just the quick ones, or all of them.

script to run mcmc

following in the footsteps of fit-model.R

compare full likelihood for TASEP to analytic equation

The likelihood for N particles is found in Schütz 1997: the probability for particles at A_n to end up at B_n after time t is:

script to obtain/view residuals

following in the footsteps of fit-model.R and the like.

unable to do inference under TASEP

After running:

Rscript sim-tasep.R -s 50 -t 20

and then

Rscript tasep-inference.R -i tasep-sims/selsims-2014-08-03-10-25-0249890.RData -o logs

I get this in the .Rout file (it seems like it should be ending up in logs)?

Loading required package: Matrix
Loading required package: methods
Loading required package: expm

Attaching package: ‘expm’

The following object is masked from ‘package:Matrix’:

    expm

Loading required package: stringdist
Loading required package: parallel
Loading required package: Biostrings
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following object is masked from ‘package:Matrix’:

    as.vector

The following object is masked from ‘package:stats’:

    xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,
    duplicated, eval, evalq, Filter, Find, get, intersect, is.unsorted,
    lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int,
    pmin, pmin.int, Position, rank, rbind, Reduce, rep.int, rownames,
    sapply, setdiff, sort, table, tapply, union, unique, unlist

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:Matrix’:

    expand

Loading required package: XVector
Loading required package: mcmc
Error in as.vector(x, "character") : 
  cannot coerce type 'closure' to vector of type 'character'
Calls: paste ... eval -> .Method -> as.character -> as.character.default
No traceback available 
Error in .Method(..., sep = sep, collapse = collapse) : 
  cannot coerce type 'closure' to vector of type 'character'
Calls: paste -> standardGeneric -> eval -> eval -> eval -> .Method
8: as.character.default(function (x) 
   {
       standardGeneric("lwin")
   })
7: as.character(function (x) 
   {
       standardGeneric("lwin")
   })
6: .Method(..., sep = sep, collapse = collapse)
5: eval(expr, envir, enclos)
4: eval(.dotsCall, env)
3: eval(.dotsCall, env)
2: standardGeneric("paste")
1: paste(basedir, "/win-", lwin, "-", win, "-", rwin, sep = "")
Error in .Method(..., sep = sep, collapse = collapse) : 
  cannot coerce type 'closure' to vector of type 'character'
Calls: paste -> standardGeneric -> eval -> eval -> eval -> .Method
6: .Method(..., sep = sep, collapse = collapse)
5: eval(expr, envir, enclos)
4: eval(.dotsCall, env)
3: eval(.dotsCall, env)
2: standardGeneric("paste")
1: paste(basename, "-results.RData", sep = "")
Error in .Method(..., sep = sep, collapse = collapse) : 
  cannot coerce type 'closure' to vector of type 'character'
Calls: paste -> standardGeneric -> eval -> eval -> eval -> .Method
6: .Method(..., sep = sep, collapse = collapse)
5: eval(expr, envir, enclos)
4: eval(.dotsCall, env)
3: eval(.dotsCall, env)
2: standardGeneric("paste")
1: paste(basename, "-results.tsv", sep = "")
Error in lwin + win : non-numeric argument to binary operator
6: .Method(..., sep = sep, collapse = collapse)
5: eval(expr, envir, enclos)
4: eval(.dotsCall, env)
3: eval(.dotsCall, env)
2: standardGeneric("paste")
1: paste(basename, "-plot", sep = "")
Error in patlen + lwin : non-numeric argument to binary operator
Calls: meangenmatrix
No traceback available 
1: meangenmatrix(lwin = 1, rwin = 1, patlen = winlen, mutpats = mutpats, 
       selpats = selpats, mutrates = mutrates * tlen, selcoef = numeric(0), 
       boundary = "none")
No traceback available 
Error in nchar(ipatterns[1]) : 
  error in evaluating the argument 'x' in selecting a method for function 'nchar': Error in rownames(genmatrix) : 
  error in evaluating the argument 'x' in selecting a method for function 'rownames': Error: object 'genmatrix' not found
Calls: collapsepatmatrix -> nchar
1: rownames(genmatrix)
2: nchar(ipatterns[1])
1: collapsepatmatrix(ipatterns = rownames(genmatrix), lwin = lwin, 
       rwin = rwin)
Error in rowSums(genmatrix) : 
  error in evaluating the argument 'x' in selecting a method for function 'rowSums': Error: object 'genmatrix' not found
Calls: computetransmatrix -> rowSums
No traceback available 
2: rowSums(genmatrix)
1: computetransmatrix(genmatrix, projmatrix, names = TRUE)
No traceback available 
Error in nchar(ipatterns[1]) : 
  error in evaluating the argument 'x' in selecting a method for function 'nchar': Error in rownames(projmatrix) : 
  error in evaluating the argument 'x' in selecting a method for function 'rownames': Error: object 'projmatrix' not found
Calls: counttrans -> nchar
1: rownames(projmatrix)
2: nchar(ipatterns[1])
1: counttrans(rownames(projmatrix), colnames(projmatrix), simseqs[[1]]$initseq, 
       simseqs[[1]]$finalseq, lwin = lwin)
Error in rowSums(counts[[1]]) : 
  error in evaluating the argument 'x' in selecting a method for function 'rowSums': Error in counts[[1]] : object of type 'closure' is not subsettable
No traceback available 
1: rowSums(counts[[1]])
Error in rownames(counts[[1]]) : 
  error in evaluating the argument 'x' in selecting a method for function 'rownames': Error in counts[[1]] : object of type 'closure' is not subsettable
Calls: leftchanged -> rownames
No traceback available 
Error in counts[[1]] : object of type 'closure' is not subsettable
2: rownames(counts[[1]])
1: leftchanged(rownames(counts[[1]]), colnames(counts[[1]]), lwin = lwin, 
       win = win)
No traceback available 
Error in as.vector(G@muttrans %*% mutrates) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error: object 'genmatrix' not found
Calls: optimize -> <Anonymous> -> f -> update -> as.vector
No traceback available 
Error in eval(expr, envir, enclos) : object 'ans' not found
Calls: data.frame -> rbind -> standardGeneric -> eval -> eval -> eval
5: as.vector(G@muttrans %*% mutrates)
4: update(genmatrix, mutrates = mutrates, selcoef = numeric(0))
3: f(arg, ...)
2: (function (arg) 
   f(arg, ...))(3.81966011868139)
1: optimize(likfun, lower = 1e-08, upper = 10)
Error in colnames(estimates) <- "muttime" : object 'estimates' not found
6: eval(expr, envir, enclos)
5: eval(quote(list(...)), env)
4: eval(quote(list(...)), env)
3: standardGeneric("rbind")
2: rbind(ans = ans$minimum, truth = truth)
1: data.frame(rbind(ans = ans$minimum, truth = truth))
Error in apply(estimates, 1, likfun) : object 'estimates' not found
No traceback available 
Error in is.data.frame(x) : object 'estimates' not found
Calls: write.table -> is.data.frame
1: apply(estimates, 1, likfun)
Error in system.time(out <- .Call("metrop", func1, initial, nbatch, blen,  : 
  object 'estimates' not found
Calls: metrop -> metrop.function -> system.time
2: is.data.frame(x)
1: write.table(estimates, file = resultsfile, quote = FALSE, sep = "\t")
Timing stopped at: 0.001 0 0 
3: system.time(out <- .Call("metrop", func1, initial, nbatch, blen, 
       nspac, scale, func2, debug, env1, env2, PACKAGE = "mcmc"))
2: metrop.function(lud, initial = estimates["ans", "muttime"], nbatch = nbatches, 
       blen = blen, scale = 0.01)
1: metrop(lud, initial = estimates["ans", "muttime"], nbatch = nbatches, 
       blen = blen, scale = 0.01)
No traceback available 
Error in lapply(1:nrow(estimates), function(k) { : 
  error in evaluating the argument 'X' in selecting a method for function 'lapply': Error in nrow(estimates) : 
  error in evaluating the argument 'x' in selecting a method for function 'nrow': Error: object 'estimates' not found
1: nrow(estimates)
1: lapply(1:nrow(estimates), function(k) {
       x <- unlist(estimates[k, ])
       list(predictcounts(win, lwin, rwin, initcounts = rowSums(counts[[1]]), 
           mutrates = x[1], selcoef = numeric(0), genmatrix = genmatrix, 
           projmatrix = projmatrix))
   })
Error in rownames(estimates) : 
  error in evaluating the argument 'x' in selecting a method for function 'rownames': Error: object 'estimates' not found
No traceback available 
Error in lwin - lcountwin : non-numeric argument to binary operator
Calls: projectcounts
1: rownames(estimates)
1: projectcounts(lwin = lwin, countwin = cwin, lcountwin = 0, rcountwin = 0, 
       counts = counts[[1]])
Error in lapply(all.expected, function(x) { : 
  error in evaluating the argument 'X' in selecting a method for function 'lapply': Error: object 'all.expected' not found
No traceback available 
Error in save(counts, genmatrix, projmatrix, subtransmatrix, lud, likfun,  : 
  object 'genmatrix' not found
1: lapply(all.expected, function(x) {
       list(projectcounts(lwin = lwin, countwin = cwin, lcountwin = 0, 
           rcountwin = 0, counts = x[[1]]))
   })
Error in eval(expr, envir, enclos) : object 'plotfile' not found
Calls: pdf ... gsub -> paste -> standardGeneric -> eval -> eval -> eval
1: save(counts, genmatrix, projmatrix, subtransmatrix, lud, likfun, 
       truth, ans, estimates, nonoverlapping, nov.counts, mmean, 
       all.expected, cwin, subcounts, all.subexpected, mrun, win, 
       lwin, rwin, file = datafile)
Error in plot(mrun$batch, type = "l", xlab = "mcmc gens", ylab = "tlen*jumprate") : 
  object 'mrun' not found
8: eval(expr, envir, enclos)
7: eval(quote(list(...)), env)
6: eval(quote(list(...)), env)
5: standardGeneric("paste")
4: paste(plotfile, "-mcmc.pdf", sep = "")
3: gsub("%%", "", s)
2: checkIntFormat(file)
1: pdf(file = paste(plotfile, "-mcmc.pdf", sep = ""), width = 6, 
       height = 4, pointsize = 10)
Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : 
  plot.new has not been called yet
Calls: abline -> int_abline
1: plot(mrun$batch, type = "l", xlab = "mcmc gens", ylab = "tlen*jumprate")
Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : 
  object 'estimates' not found
Calls: abline -> int_abline
2: int_abline(a = a, b = b, h = h, v = v, untf = untf, ...)
1: abline(h = truth, lwd = 2)
Error in strwidth(legend, units = "user", cex = cex, font = text.font) : 
  plot.new has not been called yet
Calls: legend -> strwidth
2: int_abline(a = a, b = b, h = h, v = v, untf = untf, ...)
1: abline(h = estimates["ans", ], lty = 2)
null device 
          1 
Error in eval(expr, envir, enclos) : object 'plotfile' not found
Calls: pdf ... gsub -> paste -> standardGeneric -> eval -> eval -> eval
2: strwidth(legend, units = "user", cex = cex, font = text.font)
1: legend("topright", lty = 1:2, legend = c("truth", "point estimate"))
8: eval(expr, envir, enclos)
7: eval(quote(list(...)), env)
6: eval(quote(list(...)), env)
5: standardGeneric("paste")
4: paste(plotfile, "-1.pdf", sep = "")
3: gsub("%%", "", s)
2: checkIntFormat(file)
1: pdf(file = paste(plotfile, "-1.pdf", sep = ""), width = 6, height = 4, 
       pointsize = 10)
Error in ncol(counts[[1]]) : 
  error in evaluating the argument 'x' in selecting a method for function 'ncol': Error in counts[[1]] : object of type 'closure' is not subsettable
No traceback available 
null device 
          1 
Error in eval(expr, envir, enclos) : object 'plotfile' not found
Calls: pdf ... gsub -> paste -> standardGeneric -> eval -> eval -> eval
1: ncol(counts[[1]])
Error in rainbow(2 + length(all.expected)) : 
  object 'all.expected' not found
8: eval(expr, envir, enclos)
7: eval(quote(list(...)), env)
6: eval(quote(list(...)), env)
5: standardGeneric("paste")
4: paste(plotfile, "-2.pdf", sep = "")
3: gsub("%%", "", s)
2: checkIntFormat(file)
1: pdf(file = paste(plotfile, "-2.pdf", sep = ""), width = 6, height = 4, 
       pointsize = 10)
Error in eval(expr, envir, enclos) : object 'all.subexpected' not found
Calls: order -> standardGeneric -> eval -> eval -> eval
1: rainbow(2 + length(all.expected))
null device 
          1 
Error in eval(expr, envir, enclos) : object 'plotfile' not found
Calls: pdf ... gsub -> paste -> standardGeneric -> eval -> eval -> eval
5: eval(expr, envir, enclos)
4: eval(quote(list(...)), env)
3: eval(quote(list(...)), env)
2: standardGeneric("order")
1: order(all.subexpected[["truth"]][[1]][, k])
8: eval(expr, envir, enclos)
7: eval(quote(list(...)), env)
6: eval(quote(list(...)), env)
5: standardGeneric("paste")
4: paste(plotfile, "-3.pdf", sep = "")
3: gsub("%%", "", s)
2: checkIntFormat(file)
1: pdf(file = paste(plotfile, "-3.pdf", sep = ""), width = 6, height = 4, 
       pointsize = 10)
Error in as.vector(all.expected[["truth"]][[1]]) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error: object 'all.expected' not found
Calls: plot -> as.vector
No traceback available 
Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...) : 
  plot.new has not been called yet
Calls: abline -> int_abline
2: as.vector(all.expected[["truth"]][[1]])
1: plot(as.vector(all.expected[["truth"]][[1]]), as.vector(counts[[1]]), 
       log = "xy", xlab = "true expected counts", ylab = "counts")
2: int_abline(a = a, b = b, h = h, v = v, untf = untf, ...)
1: abline(0, 1)
Error in as.vector(all.expected[["truth"]][[1]]) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error: object 'all.expected' not found
Calls: points -> as.vector
No traceback available 
Error in strwidth(legend, units = "user", cex = cex, font = text.font) : 
  plot.new has not been called yet
Calls: legend -> strwidth
2: as.vector(all.expected[["truth"]][[1]])
1: points(as.vector(all.expected[["truth"]][[1]]), as.vector(all.expected[["ans"]][[1]]), 
       col = "green", pch = 20, cex = 0.5)
null device 
          1 
[1] "2014-08-03-10-27"

switch sign so that 'likfun' is actually log likelihood

How to avoid parameter sprawl?

An explicit penalty (e.g. L^1) as part of the likelihood
Shrinkage, Bayesian style (e.g. this post)
Adding parameters one at a time and doing sequential model selection

We won't be able to do number 2 fully, but perhaps it's a framework to start from?

I'd like to understand if you are thinking of all mutations as coming from these context-dependent processes, or if they are sitting on top of some base mutation process, like in your CpG example.

test json-tree-cpg

And, I'll take you up on your previous offer: I would like to test
json-tree-cpg in a similar way to json-cpg; could you add that in? I
would, but it isn't as straightforward as adding another file to
nestly/models/, as some of the calls and options are different:
fit-tree-model.R rather than fit-model.R, and make-genmat takes no -o
flag, as the location of the genmats is specified in the config file.

make `boundary` flag in JSON to use meangenmatrix

specialize to the case of an observed X(0)

Maybe random X(0) later.

investigate multiple peaks in likelihood surface

Want to investigate the multiple peaks thing? Simulate one set of
data, start a bunch of optimizations from different starting points?
You'll have to make sure the scaling parameters are reasonable, first.

@petrelharp , do you have some systematic exploration of starting points in mind, or were you thinking of just randomly perturbing the starting points used for simulation, say, multiplicatively?

WRT scale, it looks like the code will fill in something if I don't provide parameters. Were you saying that I should make sure that these were reasonable?

make sure peel-transmat gets columns in same order as counts file

... and refactor to precompute the tree traversal.

add a 'context.mcmc' or somesuch object carrying around relevant Bayesian whatnot

get better guessing or input of tuning parameters

namely, the parscale parameter for optim, and the stepsize parameter for mcmc.

write predictcounts and projectcounts for tree counts

remove (or, fixup?) `nestly/`?

The build_and_test.sh script does some scons + nestly things that fail:

$ cd nestly && scons simple && scons -j 6 seed
bash: cd: nestly: No such file or directory
peter@skua:~/projects/context/nestly$ scons simple
scons: Reading SConscript files ...
TypeError: expected str, bytes or os.PathLike object, not dict:
  File "/home/peter/projects/context/nestly/SConstruct", line 205:
    nest = cl_targets[target_str]()      # Initialize the nest using one of the functions above
  File "/home/peter/projects/context/nestly/SConstruct", line 67:
    nest.add('model', model_configs,
  File "/home/peter/.local/lib/python3.9/site-packages/nestly/scons.py", line 94:
    return self.nest.add(name, nestable, **kw)
  File "/home/peter/.local/lib/python3.9/site-packages/nestly/core.py", line 194:
    new_outdir = os.path.join(outdir, label_func(to_label))
  File "/usr/lib/python3.9/posixpath.py", line 142:
    p = os.fspath(p)

We should remove this, I assume, unless there is a very easy fix?

systematic checking of model inference with simulated counts

My first thought would be to do this through the command line, via nestly.

count-seq.R expects counttrans to return a matrix

workflow.sh runs count-seq.R like so:

Rscript ../count-seq.R -i sim-cpg-123456.RData -w 4 -s 2 -l 1

which has these lines:

# this returns a matrix
counts <- if (revcounts) {
    counttrans( longpats, shortpats, simseqs[[1]]$initseq, simseqs[[1]]$finalseq, lwin=countseq.opt$lwin )
} else {
    counttrans( longpats, shortpats, simseqs[[1]]$finalseq, simseqs[[1]]$initseq, lwin=countseq.opt$lwin )
}

and indeed it seems like the later code expects it to be a matrix or data frame:

countframe <- data.frame( reference=rownames(counts)[row(counts)], 
                         derived=colnames(counts)[col(counts)],
                         count=as.vector(counts)
                         )

However, counttrans says

    counts <- lapply( counts, function (x) new("tuplecounts",counts=x,lwin=lwin) )

and R doesn't want to convert from that object:

Error in as.vector(counts) : 
  no method for coercing this S4 class to a vector
Calls: data.frame -> as.vector -> as.vector
Execution halted

write down curvature of likelihood surface

... or something, to work out how confounded different parameterizations are.

with small data, adhoc guess at rate parameters may be NaN

add an error check

Add flag for setting random seed in `sim-seq.R` from the command line.

So we can have reproducible sims.

Also include in the SConstruct initial dictionary.

padding, wings, fringe?

Shall we decide on a term?

prove something about approximation for full likelihood?

This is harder than the case where we don't condition on the final sequence, as demonstrated by the following two examples:

TASEP: let x be a string of 0s with a 1 at the left-hand end, and y be a string of 0s with a 1 at the right-hand end. If both these strings are of length n, then p(y[n] | y[k:n-1], x[k:n]) is zero (or, small, if entering from the boundary is allowed), while p(y[n] | y[1:n-1], x[1:n]) is in fact very large, since conditioning on there being a particle at position 1 at the start, and there being no particles at positions y[1:n-1] implies that the particle must be at y[n] (or further to the right). So - we get totally the wrong answer by conditioning on a shorter sequence... but such situations are also very unlikely.
See attached photo. The two sequences (AGC -> CTT) might have happened through two scenarios (suppose), shown in top and bottom. Rates are on the right. This is an example where a mutation (the A->C in the leftmost site) affects the probability of an earlier one (the G->C in the top scenario). This means that information is propagating backwards in time, as well as forwards, and the conditional probability of the rightmost site could be very different depending on whether the outcome at the leftmost site is included in the conditioning or not.

add more stuff to json

e.g. tlen for sim-seq so that everything we need to check the answers is in there

Writeup: do we mind a collision of ell's?

\ell is used both for motif length and for left margin.

We would like the likfun function to automatically have access to the stuff in the context object

In context-inference-fns.R...

# We would like the likfun function to automatically have access to the stuff in the context object
#  ... where is 'self'?!?
# EVIL:
setGeneric("likfun", function (x) { standardGeneric("likfun") } )
setMethod("likfun", signature=c(x="context"), definition=function(x) {
          f <- x@likfun
          environment(f) <- list2env( list(genmatrix=x@genmatrix,projmatrix=x@projmatrix,counts=x@data), parent=globalenv())
          return(f) } )

self is there, no? (Looking at http://practicalcomputing.org/node/80). It seems that the issue comes about if we try to isolate the likfun from its context object. Why is it important that you do that?

remove -g options

since Rscript -g doesn't work.

`getmutpats` errors out

Specifically, it seems to have an issue when going through the loop at k=1. Here's what I tried evaluating to isolate that case.

> mutpats <- list()
> patterns <- getpatterns(patlen,bases)
> k <- 1
> kmers <- getpatterns(k,bases)
> mutpats <- c( mutpats,
+         apply(combn(kmers,2),2,list), # make lists of rows (2 = apply over colment lists of kmers 
+         apply(combn(kmers,2)[2:1,],2,list) # and the reverse of the 2-element lists
+         )
Error in apply(combn(kmers, 2)[2:1, ], 2, list) : 
  dim(X) must have a positive length
>

which is the same result as

> bases <- c("X","O")
> getmutpats(2,bases)
Error in apply(combn(kmers, 2)[2:1, ], 2, list) : 
  dim(X) must have a positive length

It appears to me that this comes from R's implicit casting rule for lists of length 1.

mutation patterns
window sizes
selection patterns

rather than editing one of those gigantic R scripts.

petrelharp / context Goto Github PK

context's People

Contributors

Watchers

context's Issues

Recommend Projects

Recommend Topics

Recommend Org