adrimaz / flan Goto Github PK

Dear all,

I was having a look at using flan for estimating mutation rates, from my data. I usually have my data in a dataframe format, and would like to estimate the mutation rate for a number of clones. So my data will look something like

genotype; mutants; final_counts; plating_efficiency
A; 0; 1e9; 1
A; 10; 3e9; 1
A; 5; 2e9; 1
B; 1; 1e9; 0.3
B; 1; 2e9; 0.3
B; 50; 3e9; 0.3

I can easily calculate the mutation rate for each clone with:

library(data.table) as.data.table(data)[, mutestim(mutants,final_counts,method="GF"), by = genotype]

However, as you can see from the mock data, sometimes I will have variation (per genotype) for the plating efficiency. Such that for genotype A, I will have plated the entire culture, but for genotype B, I will have plated just 1/3.

So the problem is when I run:

as.data.table(data)[, mutestim(mutants,final_counts,plateeff=plating_efficiency,method="GF"), by = genotype]

I get the error:

Error in mutestim(res, final, method = "GF", plateff = plating_efficiency) : 'plateff' must be a single positive and <= 1 number. In addition: Warning message: In if (plateff > 1 | length(plateff) > 1) { : the condition has length > 1 and only the first element will be used

So is there a way around this? or would it be possible for a future version of flan to allow plateeff to be chosen from a data frame. I guess this would be useful for a lot of people.

Very happy to discover flan recently and I've been trying to use it to work out what effect the issue of variance in Nt (raised in the 2014 PLOSone paper Ycart B, Veziris N (2014) Unbiased estimation of mutation rates under fluctuating final counts. PLoS One 9: e101434.) might have on our results.

However, there's some behaviour I don't understand that looks problematic (though it could just be my not having understood it correctly) – if you give mutestim() a vector of counts, a CV for the Nt and nothing else, it comes back with an estimate of the mutation probability, rather than an estimate of the number of mutations, which is what I would have expected when it hasn't been told the actual value of Nt. So:

library(flan)
x <- c(0, 0, 0, 2, 2, 1, 1, 0, 7, 0, 1, 0, 1, 2, 0, 5)
mutestim(x, cvfn = 0.2)

gives
$mutprob
[1] 8.319557e-10

$sd.mutprob
[1] 2.802634e-10

$fitness
[1] 2.34619

$sd.fitness
[1] 1.626095

Which, oddly, looks very similar to the result of:
mutestim(x, fn = rep(9.909846e8, length(x)), cvfn = 0.2)

It seems unlikely that the algorithm has silently estimated an Nt value of 9.9e8 (even though this is a reasonably realistic value), so I might have expected instead a result with $mutations as the first element of the result and perhaps a warning to say the cvfn value was being ignored.

A supplementary question related to this (and possibly clarifying my (mis-)understanding):
– the sort of data we have has Nt estimated separately from the cultures giving the mutant counts, but the estimate comes from at least 3 independent cultures (usually a single count from each one), so it is possible to calculate a CV for the Nt estimate. Is it therefore valid to use the sort of call shown above, repeating the single estimated Nt value (say 9.909846e8) from those counts in the fn argument to mutestim(), as if each one had been estimated separately, and using the estimated CV (say 0.2) in the cvfn argument, and so use the estimates given by mutestim(x, fn = rep(9.909846e8, length(x)), cvfn = 0.2)?

Thanks!
Chris