diff --git a/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/RandomAlignment.java b/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/RandomAlignment.java new file mode 100644 index 0000000..b288bcf --- /dev/null +++ b/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/RandomAlignment.java @@ -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 { + +} diff --git a/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/TreeLikelihoodWithReadCounts.java b/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/TreeLikelihoodWithReadCounts.java new file mode 100644 index 0000000..7f5a5ab --- /dev/null +++ b/phylonco-beast/src/main/java/phylonco/beast/evolution/readcountmodel/TreeLikelihoodWithReadCounts.java @@ -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 epsilonInput = new Input<>("epsilon", "sequencing error"); + + final public Input deltaInput = new Input<>("delta", "allelic dropout probability"); + + // other parameters of read count model + + final public Input 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 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; + } + + + + + +} diff --git a/phylonco-lphy/src/main/java/phylonco/lphy/evolution/readcountmodel/ReadCountModel.java b/phylonco-lphy/src/main/java/phylonco/lphy/evolution/readcountmodel/ReadCountModel.java index d432c18..6cdc43a 100644 --- a/phylonco-lphy/src/main/java/phylonco/lphy/evolution/readcountmodel/ReadCountModel.java +++ b/phylonco-lphy/src/main/java/phylonco/lphy/evolution/readcountmodel/ReadCountModel.java @@ -183,7 +183,7 @@ public RandomVariable sample() { readC[i][j] = multinomial.sample().value(); } }else { - pA = 0.5 - eps / 3; pC = 0.5 - eps / 3; pG = eps / 3; pT = eps / 3; + pA = 0.5 - eps / 6; pC = 0.5 - eps / 6; pG = eps / 6; pT = eps / 6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat); @@ -207,7 +207,7 @@ public RandomVariable sample() { readC[i][j] = multinomial.sample().value(); } }else { - pA = 0.5-eps/3; pC = eps/3; pG = 0.5-eps/3; pT = eps/3; + pA = 0.5-eps/6; pC = eps/6; pG = 0.5-eps/6; pT = eps/6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat); @@ -231,7 +231,7 @@ public RandomVariable sample() { readC[i][j] = multinomial.sample().value(); } }else { - pA = 0.5 - eps / 3; pC = eps / 3; pG = eps / 3; pT = 0.5 - eps / 3; + pA = 0.5 - eps / 6; pC = eps / 6; pG = eps / 6; pT = 0.5 - eps / 6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat); @@ -263,7 +263,7 @@ else if (genotype.equals("CG") || genotype.equals("GC")){ readC[i][j] = multinomial.sample().value(); } } else { - pA = eps / 3; pC = 0.5 - eps / 3; pG = 0.5 - eps / 3; pT = eps / 3; + pA = eps / 6; pC = 0.5 - eps / 6; pG = 0.5 - eps / 6; pT = eps / 6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat); @@ -287,7 +287,7 @@ else if (genotype.equals("CG") || genotype.equals("GC")){ readC[i][j] = multinomial.sample().value(); } } else { - pA = eps / 3; pC = 0.5 - eps / 3; pG = eps / 3; pT = 0.5 - eps / 3; + pA = eps / 6; pC = 0.5 - eps / 6; pG = eps / 6; pT = 0.5 - eps / 6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat); @@ -320,7 +320,7 @@ else if (genotype.equals("GT") || genotype.equals("TG")){ readC[i][j] = multinomial.sample().value(); } }else { - pA = eps / 3; pC = eps / 3; pG = 0.5 - eps / 3; pT = 0.5 - eps / 3; + pA = eps / 6; pC = eps / 6; pG = 0.5 - eps / 6; pT = 0.5 - eps / 6; proMatrix = new Double[]{pA, pC, pG, pT}; Value proMat = new Value<>("proMatrix", proMatrix); multinomial.setParam("p", proMat);