-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
fix bug in read count simulator and skeleton code for beast likelihood …
- Loading branch information
1 parent
0bdf1c8
commit 81cda54
Showing
3 changed files
with
167 additions
and
6 deletions.
There are no files selected for viewing
12 changes: 12 additions & 0 deletions
12
phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/RandomAlignment.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,12 @@ | ||
| package phylonco.beast.evolution.readcountmodel; | ||
|
|
||
| import beast.base.inference.State; | ||
|
|
||
|
|
||
| // readcount datatype (string -> to array of observations) | ||
| // random alignment matrix of observations (generic observations) | ||
| // sequences as a state, or alignment as a state | ||
|
|
||
| public class RandomAlignment extends State { | ||
|
|
||
| } |
149 changes: 149 additions & 0 deletions
149
...t/src/main/java/phylonco/beast/evolution/readcountmodel/TreeLikelihoodWithReadCounts.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,149 @@ | ||
| package phylonco.beast.evolution.readcountmodel; | ||
|
|
||
| import beast.base.core.Input; | ||
| import beast.base.evolution.likelihood.TreeLikelihood; | ||
| import beast.base.evolution.tree.Node; | ||
| import beast.base.evolution.tree.Tree; | ||
| import beast.base.inference.parameter.RealParameter; | ||
| import phylonco.beast.evolution.datatype.ReadCount; | ||
|
|
||
| import java.util.List; | ||
|
|
||
| public class TreeLikelihoodWithReadCounts extends TreeLikelihood { | ||
|
|
||
| // input data is an alignment with the GT16 datatype | ||
|
|
||
| // epsilon, allelic dropout, ... parameters | ||
| final public Input<RealParameter> epsilonInput = new Input<>("epsilon", "sequencing error"); | ||
|
|
||
| final public Input<RealParameter> deltaInput = new Input<>("delta", "allelic dropout probability"); | ||
|
|
||
| // other parameters of read count model | ||
|
|
||
| final public Input<ReadCount> readCountdataInput = new Input<>("readcount", "read count data for the beast.tree", Input.Validate.REQUIRED); | ||
|
|
||
| private RealParameter epsilon; | ||
|
|
||
| @Override | ||
| public void initAndValidate() { | ||
| super.initAndValidate(); | ||
| // your initialization (input checking for validity) | ||
|
|
||
| // check Alignment type is the correct datatype (GT16) | ||
|
|
||
| // set up likelihood core for memory allocation | ||
|
|
||
| this.epsilon = epsilonInput.get(); | ||
| } | ||
|
|
||
| @Override | ||
| public double calculateLogP() { | ||
| // get the values from input parameters | ||
| traverse(treeInput.get().getRoot()); | ||
|
|
||
| // calculate log likelihood | ||
|
|
||
| return 0.0; | ||
| } | ||
|
|
||
| /* Assumes there IS a branch rate model as opposed to traverse() */ | ||
| protected int traverse(final Node node) { | ||
|
|
||
| int update = (node.isDirty() | hasDirt); | ||
|
|
||
| final int nodeIndex = node.getNr(); | ||
|
|
||
| final double branchRate = branchRateModel.getRateForBranch(node); | ||
| final double branchTime = node.getLength() * branchRate; | ||
|
|
||
| // First update the transition probability matrix(ices) for this branch | ||
| //if (!node.isRoot() && (update != Tree.IS_CLEAN || branchTime != m_StoredBranchLengths[nodeIndex])) { | ||
| if (!node.isRoot() && (update != Tree.IS_CLEAN || branchTime != m_branchLengths[nodeIndex])) { | ||
| m_branchLengths[nodeIndex] = branchTime; | ||
| final Node parent = node.getParent(); | ||
| likelihoodCore.setNodeMatrixForUpdate(nodeIndex); | ||
| for (int i = 0; i < m_siteModel.getCategoryCount(); i++) { | ||
| final double jointBranchRate = m_siteModel.getRateForCategory(i, node) * branchRate; | ||
| substitutionModel.getTransitionProbabilities(node, parent.getHeight(), node.getHeight(), jointBranchRate, probabilities); | ||
| //System.out.println(node.getNr() + " " + Arrays.toString(m_fProbabilities)); | ||
| likelihoodCore.setNodeMatrix(nodeIndex, i, probabilities); | ||
| } | ||
| update |= Tree.IS_DIRTY; | ||
| } | ||
|
|
||
| // If the node is internal, update the partial likelihoods. | ||
| if (!node.isLeaf()) { | ||
|
|
||
| // Traverse down the two child nodes | ||
| final Node child1 = node.getLeft(); //Two children | ||
| final int update1 = traverse(child1); | ||
|
|
||
| final Node child2 = node.getRight(); | ||
| final int update2 = traverse(child2); | ||
|
|
||
| // If either child node was updated then update this node too | ||
| if (update1 != Tree.IS_CLEAN || update2 != Tree.IS_CLEAN) { | ||
|
|
||
| final int childNum1 = child1.getNr(); | ||
| final int childNum2 = child2.getNr(); | ||
|
|
||
| likelihoodCore.setNodePartialsForUpdate(nodeIndex); | ||
| update |= (update1 | update2); | ||
| if (update >= Tree.IS_FILTHY) { | ||
| likelihoodCore.setNodeStatesForUpdate(nodeIndex); | ||
| } | ||
|
|
||
| if (m_siteModel.integrateAcrossCategories()) { | ||
| likelihoodCore.calculatePartials(childNum1, childNum2, nodeIndex); | ||
| } else { | ||
| throw new RuntimeException("Error TreeLikelihood 201: Site categories not supported"); | ||
| //m_pLikelihoodCore->calculatePartials(childNum1, childNum2, nodeNum, siteCategories); | ||
| } | ||
|
|
||
| if (node.isRoot()) { | ||
| // No parent this is the root of the beast.tree - | ||
| // calculate the pattern likelihoods | ||
|
|
||
| final double[] proportions = m_siteModel.getCategoryProportions(node); | ||
| likelihoodCore.integratePartials(node.getNr(), proportions, m_fRootPartials); | ||
|
|
||
| // get constant pattern | ||
| List<Integer> constantPattern = getConstantPattern(); | ||
|
|
||
| if (constantPattern != null) { // && !SiteModel.g_bUseOriginal) { | ||
| proportionInvariant = m_siteModel.getProportionInvariant(); | ||
| // some portion of sites is invariant, so adjust root partials for this | ||
| for (final int i : constantPattern) { | ||
| m_fRootPartials[i] += proportionInvariant; | ||
| } | ||
| } | ||
|
|
||
| double[] rootFrequencies = substitutionModel.getFrequencies(); | ||
| if (rootFrequenciesInput.get() != null) { | ||
| rootFrequencies = rootFrequenciesInput.get().getFreqs(); | ||
| } | ||
| likelihoodCore.calculateLogLikelihoods(m_fRootPartials, rootFrequencies, patternLogLikelihoods); | ||
| } | ||
|
|
||
| } | ||
| } else { | ||
| // leaf node | ||
| // call method to calculate leaf likelihood with read count model | ||
| } | ||
| return update; | ||
| } // traverseWithBRM | ||
|
|
||
| @Override | ||
| protected boolean requiresRecalculation() { | ||
|
|
||
| boolean requiresRecal = super.requiresRecalculation(); | ||
|
|
||
| // check whether the other parameters are "dirty" or have been changed | ||
| return true; | ||
| } | ||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters