/* * HypermutantErrorModel.java * * Copyright (c) 2002-2016 Alexei Drummond, Andrew Rambaut and Marc Suchard * * This file is part of BEAST. * See the NOTICE file distributed with this work for additional * information regarding copyright ownership and licensing. * * BEAST is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with BEAST; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA */ package dr.evomodel.tipstatesmodel; import dr.evolution.alignment.HypermutantAlignment; import dr.evolution.datatype.Nucleotides; import dr.inference.model.Parameter; import dr.inference.model.Statistic; import dr.inference.model.Variable; import dr.util.Author; import dr.util.Citable; import dr.util.Citation; import dr.xml.*; import java.util.Arrays; import java.util.HashMap; import java.util.List; import java.util.Map; import java.util.logging.Logger; /** * @author Andrew Rambaut * @version $Id$ */ public class HypermutantErrorModel extends TipStatesModel implements Citable { public static final String HYPERMUTANT_ERROR_MODEL = "hypermutantErrorModel"; public static final String HYPERMUTATION_RATE = "hypermutationRate"; public static final String HYPERMUTATION_INDICATORS = "hypermutationIndicators"; public static final String UNLINKED_RATES = "unlinkedRates"; public HypermutantErrorModel(HypermutantAlignment hypermutantAlignment, Parameter hypermutationRateParameter, Parameter hypermuationIndicatorParameter, boolean unlinkedRates) { super(HYPERMUTANT_ERROR_MODEL, null, null); this.hypermutantAlignment = hypermutantAlignment; this.unlinkedRates = unlinkedRates; this.hypermutationRateParameter = hypermutationRateParameter; addVariable(this.hypermutationRateParameter); this.hypermutationIndicatorParameter = hypermuationIndicatorParameter; addVariable(this.hypermutationIndicatorParameter); addStatistic(new TaxonHypermutatedStatistic()); addStatistic(new TaxonHypermutationRateStatistic()); addStatistic(new HypermutatedProportionStatistic()); } protected void taxaChanged() { if (hypermutationIndicatorParameter.getDimension() <= 1) { this.hypermutationIndicatorParameter.setDimension(tree.getExternalNodeCount()); } if (unlinkedRates && hypermutationRateParameter.getDimension() <= 1) { this.hypermutationRateParameter.setDimension(tree.getExternalNodeCount()); } } @Override public Type getModelType() { return Type.PARTIALS; } @Override public void getTipStates(int nodeIndex, int[] tipStates) { throw new IllegalArgumentException("This model emits only tip partials"); } @Override public void getTipPartials(int nodeIndex, double[] partials) { int[] states = this.states[nodeIndex]; boolean isHypermutated = hypermutationIndicatorParameter.getParameterValue(nodeIndex) > 0.0; double rate = (unlinkedRates ? hypermutationRateParameter.getParameterValue(nodeIndex) : hypermutationRateParameter.getParameterValue(0)); int k = 0; for (int j = 0; j < patternCount; j++) { switch (states[j]) { case Nucleotides.A_STATE: // is an A partials[k] = 1.0; partials[k + 1] = 0.0; partials[k + 2] = 0.0; partials[k + 3] = 0.0; break; case Nucleotides.C_STATE: // is an C partials[k] = 0.0; partials[k + 1] = 1.0; partials[k + 2] = 0.0; partials[k + 3] = 0.0; break; case Nucleotides.G_STATE: // is an G partials[k] = 0.0; partials[k + 1] = 0.0; partials[k + 2] = 1.0; partials[k + 3] = 0.0; break; case Nucleotides.UT_STATE: // is an T partials[k] = 0.0; partials[k + 1] = 0.0; partials[k + 2] = 0.0; partials[k + 3] = 1.0; break; case Nucleotides.R_STATE: // is an A in a APOBEC context if (isHypermutated) { partials[k] = 1.0 - rate; partials[k + 1] = 0.0; partials[k + 2] = rate; partials[k + 3] = 0.0; } else { partials[k] = 1.0; partials[k + 1] = 0.0; partials[k + 2] = 0.0; partials[k + 3] = 0.0; } break; default: // is an ambiguity partials[k] = 1.0; partials[k + 1] = 1.0; partials[k + 2] = 1.0; partials[k + 3] = 1.0; } k += stateCount; } } @Override protected final void handleVariableChangedEvent(Variable variable, int index, Parameter.ChangeType type) { if (variable == hypermutationIndicatorParameter) { fireModelChanged(tree.getTaxon(index)); } else if (variable == hypermutationRateParameter) { if (!unlinkedRates) { fireModelChanged(); } else if (hypermutationIndicatorParameter.getValue(index) > 0.5) { // only fire an update if the indicator is on.... fireModelChanged(tree.getTaxon(index)); } } else { throw new RuntimeException("Unknown parameter has changed in HypermutantErrorModel.handleVariableChangedEvent"); } } public static XMLObjectParser PARSER = new AbstractXMLObjectParser() { public String getParserName() { return HYPERMUTANT_ERROR_MODEL; } public Object parseXMLObject(XMLObject xo) throws XMLParseException { boolean unlinkedRates = false; if (xo.hasAttribute(UNLINKED_RATES)) { unlinkedRates = xo.getBooleanAttribute(UNLINKED_RATES); } HypermutantAlignment hypermutantAlignment = (HypermutantAlignment)xo.getChild(HypermutantAlignment.class); Parameter hypermutationRateParameter = null; if (xo.hasChildNamed(HYPERMUTATION_RATE)) { hypermutationRateParameter = (Parameter)xo.getElementFirstChild(HYPERMUTATION_RATE); } Parameter hypermuationIndicatorParameter = null; if (xo.hasChildNamed(HYPERMUTATION_INDICATORS)) { hypermuationIndicatorParameter = (Parameter)xo.getElementFirstChild(HYPERMUTATION_INDICATORS); } HypermutantErrorModel errorModel = new HypermutantErrorModel(hypermutantAlignment, hypermutationRateParameter, hypermuationIndicatorParameter, unlinkedRates); Logger.getLogger("dr.evomodel").info("Using APOBEC error model"); return errorModel; } //************************************************************************ // AbstractXMLObjectParser implementation //************************************************************************ public String getParserDescription() { return "This element returns a model that allows for APOBEC-type RNA editing."; } public Class getReturnType() { return HypermutantErrorModel.class; } public XMLSyntaxRule[] getSyntaxRules() { return rules; } private XMLSyntaxRule[] rules = new XMLSyntaxRule[] { AttributeRule.newBooleanRule(UNLINKED_RATES, true), new ElementRule(HypermutantAlignment.class), new ElementRule(HYPERMUTATION_RATE, Parameter.class, "The hypermutation rate per target site per sequence"), new ElementRule(HYPERMUTATION_INDICATORS, Parameter.class, "A binary indicator of whether the sequence is hypermutated"), }; }; public class TaxonHypermutatedStatistic extends Statistic.Abstract { public TaxonHypermutatedStatistic() { super("isHypermutated"); } public int getDimension() { return hypermutationIndicatorParameter.getDimension(); } public String getDimensionName(int dim) { return taxonMap.get(dim); } public double getStatisticValue(int dim) { return hypermutationIndicatorParameter.getParameterValue(dim); } } public class TaxonHypermutationRateStatistic extends Statistic.Abstract { public TaxonHypermutationRateStatistic() { super("hypermutationRate"); } public int getDimension() { return hypermutationRateParameter.getDimension(); } public String getDimensionName(int dim) { return taxonMap.get(dim) + ".rate"; } public double getStatisticValue(int dim) { return hypermutationRateParameter.getParameterValue(dim) * hypermutationIndicatorParameter.getParameterValue(dim); } } public class HypermutatedProportionStatistic extends Statistic.Abstract { public HypermutatedProportionStatistic() { super("proportionHypermutated"); } public int getDimension() { return 1; } public String getDimensionName(int dim) { return "P(hypermutated)"; } public double getStatisticValue(int dim) { if (mutatedContextCounts == null) { mutatedContextCounts = new HashMap<Integer, Integer>(); unmutatedContextCounts = new HashMap<Integer, Integer>(); for (int index : taxonMap.keySet()) { String name = taxonMap.get(index); int i = hypermutantAlignment.getTaxonIndex(name); mutatedContextCounts.put(index, hypermutantAlignment.getMutatedContextCounts()[i]); unmutatedContextCounts.put(index, hypermutantAlignment.getUnmutatedContextCounts()[i]); } } double mutatedCount = 0; double totalCount = 0; for (int i = 0; i < hypermutationIndicatorParameter.getDimension(); i++) { if (hypermutationIndicatorParameter.getParameterValue(i) > 0.5) { mutatedCount += mutatedContextCounts.get(i); totalCount += mutatedCount + unmutatedContextCounts.get(i); } } double r = hypermutationRateParameter.getParameterValue(0); return (r * mutatedCount) / totalCount; } } Map<Integer, Integer> mutatedContextCounts = null; Map<Integer, Integer> unmutatedContextCounts = null; private final HypermutantAlignment hypermutantAlignment; private final Parameter hypermutationRateParameter; private final Parameter hypermutationIndicatorParameter; private final boolean unlinkedRates; @Override public Citation.Category getCategory() { return Citation.Category.SUBSTITUTION_MODELS; } @Override public String getDescription() { return "Sequence error model"; } @Override public List<Citation> getCitations() { return Arrays.asList(new Citation( new Author[]{ new Author("A", "Rambaut"), new Author("SYW", "Ho"), new Author("AJ", "Drummond"), new Author("B", "Shapiro"), }, "Accommodating the effect of ancient DNA damage on inferences of demographic histories", 2008, "Mol Biol Evol", 26, 245, 248, "10.1093/molbev/msn256" ), new Citation( new Author[]{ new Author("J", "Felsenstein"), }, "Inferring Phylogenies", 2004, "Sinauer Associates", "" )); } }