Dear all,

I am developing a special model to allow polymorphisms in a phylogenetic analysis.
Until now I assumed that observed allele frequencies are the "true" ones, but now I want to account for sampling variance ("uncertainty", or "ambiguity" one could say). That is, observed states at leaves will not be really observed any more, and tips of the tree will be similar to inner nodes.

The perfect solution for me would be to specify the input file not the usual way: ">species1 ATTGCA...", but explicitly declaring at each site the probability of each character: ">species1 (A:0.9,C:0.05,G:0.05)TTG(A:0.05,C:0.9,G:0.05)A..." or something similar (in fact, instead of nucleotides I have codons).
I guess that to obtain this result I would have to modify the functions "ReadDataFile", "CreateFilter", probably also "LikelihoodFunction"? Does this sound feasible in a reasonable time?

The alternative option would be to introduce new characters, the "observed" states (Ao, Co, Go, To) distinct from the "unobserved" states (Au, Cu, Gu, Tu), allow only unobserved characters in the inner phylogeny, add "dummy" branches, one below each tree tip, and on those dummy branches allow only substitutions from unobserved to observed states. This should work fine, but should also slow down everything (as the state space would have increased dimension, and I am already basically working with codons).

Alternatively, is there a way to allow at most 1 substitution on a branch?


Thanks, ad best regards,

Nicola De Maio

Hi Nicola,

Yes, there is a way to define the transition matrix directly. All you do is include a flag in the Model definition as in:

Code:

Model 		MG = (TransitionMatrix, EqFreqs ,EXPLICIT_FORM_MATRIX_EXPONENTIAL);
 

Sergei

Hi Sergei,

Thank you very much!
I think that I should be able to solve my problems with the command that you suggested me! I will let you know!

One small complication is that I need to enforce the molecular clock on the original tree, but I guess it should be fine once I fix the "dummy" branches to a constant, uniform, length.

Thank you very much again! It was an incredibly fast and helpful reply!

Nicola De Maio

Hi again!

I found the error, I just had to use "matrix_probs" instead of matrix_probs.

But now I have more problems.
With the MolecularClock() function, If I set:

Code:

Tree   myTree = "((s3{MLeaf})u3,((s2{MLeaf})u2,(s1{MLeaf})u1)));";
myTree.s1.t := 1.0;
myTree.s2.t := 1.0;
myTree.s3.t := 1.0;
LikelihoodFunction theLF = (myFilter, myTree);
MolecularClock (myTree,t);
 

I get the (reasonable) error message: "Molecular clock constraint has failed, since variable t is not an independent member of the node testTree.s3".

I am sure that I can get around this in particular instances by defining the molecular clock constraints by hand, but, as I am running a topology search algorithm, I would like to have something that works on any tree.

Cheers,
Nicola

Hi Nicola,

It's a HyPhy bug (basically, the 16x16 matrix you specify gets stored in a sparse format, but the LF evaluator assumes that the matrix is dense, which it will be if derived from the rate matrix). I'll push a fix in shortly.

Sergei

PS Fix pushed to github (Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login , Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login). Please check out, recompile HyPhy and the issue should disappear.