use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class EmpiricalCodonModelParser method createNewFreqModel.
// creates new FrequencyModel from XML frequencies
private FrequencyModel createNewFreqModel(DataType codons, EmpiricalCodonRateMatrix type) throws XMLParseException {
double[] freqs = type.getFrequencies();
double sum = 0;
for (int j = 0; j < freqs.length; j++) {
sum += freqs[j];
}
if (Math.abs(sum - 1.0) > 1e-8) {
throw new XMLParseException("Frequencies do not sum to 1 (they sum to " + sum + ")");
}
FrequencyModel fm = new FrequencyModel(codons, freqs);
Logger.getLogger("dr.evomodel").info("Using frequencies from data file");
return fm;
}
use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class SeqGen method main.
public static void main(String[] argv) {
String treeFileName = argv[0];
String outputFileStem = argv[1];
int length = 500;
double[] frequencies = new double[] { 0.25, 0.25, 0.25, 0.25 };
double kappa = 10.0;
double alpha = 0.5;
double substitutionRate = argv.length < 3 ? 1.0E-3 : Double.parseDouble(argv[2]);
int categoryCount = argv.length < 4 ? 8 : Integer.parseInt(argv[3]);
//1.56E-6;
double damageRate = argv.length < 5 ? 0 : Double.parseDouble(argv[4]);
System.out.println("substitutionRate = " + substitutionRate + "; categoryCount = " + categoryCount + "; damageRate = " + damageRate);
FrequencyModel freqModel = new FrequencyModel(dr.evolution.datatype.Nucleotides.INSTANCE, frequencies);
HKY hkyModel = new HKY(kappa, freqModel);
SiteModel siteModel = null;
if (categoryCount > 1) {
siteModel = new GammaSiteModel(hkyModel, alpha, categoryCount);
} else {
// no rate heterogeneity
siteModel = new GammaSiteModel(hkyModel);
}
List<Tree> trees = new ArrayList<Tree>();
FileReader reader = null;
try {
reader = new FileReader(treeFileName);
// TreeImporter importer = new NexusImporter(reader);
TreeImporter importer = new NewickImporter(reader);
while (importer.hasTree()) {
Tree tree = importer.importNextTree();
trees.add(tree);
System.out.println("tree height = " + tree.getNodeHeight(tree.getRoot()) + "; leave nodes = " + tree.getExternalNodeCount());
}
} catch (FileNotFoundException e) {
e.printStackTrace();
return;
} catch (Importer.ImportException e) {
e.printStackTrace();
return;
} catch (IOException e) {
e.printStackTrace();
return;
}
SeqGen seqGen = new SeqGen(length, substitutionRate, freqModel, hkyModel, siteModel, damageRate);
int i = 1;
for (Tree tree : trees) {
Alignment alignment = seqGen.simulate(tree);
FileWriter writer = null;
try {
// writer = new FileWriter(outputFileStem + (i < 10 ? "00" : (i < 100 ? "0" : "")) + i + ".nex");
// NexusExporter exporter = new NexusExporter(writer);
//
// exporter.exportAlignment(alignment);
//
// writer.close();
String outputFileName = outputFileStem + "-" + substitutionRate + ".fasta";
writer = new FileWriter(outputFileName);
BufferedWriter bf = new BufferedWriter(writer);
FastaExporter exporter = new FastaExporter(bf);
exporter.exportSequences(alignment.getSequenceList());
bf.close();
System.out.println("Write " + i + "th sequence file : " + outputFileName);
i++;
} catch (IOException e) {
e.printStackTrace();
return;
}
}
}
use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class SequenceSimulator method simulate.
//END: sequence2intArray
/**
* perform the actual sequence generation
* @return alignment containing randomly generated sequences for the nodes in the
* leaves of the tree
*/
public Alignment simulate() {
NodeRef root = m_tree.getRoot();
double[] categoryProbs = m_siteModel.getCategoryProportions();
int[] category = new int[m_sequenceLength];
for (int i = 0; i < m_sequenceLength; i++) {
category[i] = MathUtils.randomChoicePDF(categoryProbs);
}
int[] seq = new int[m_sequenceLength];
if (has_ancestralSequence) {
seq = sequence2intArray(ancestralSequence);
} else {
FrequencyModel frequencyModel = m_siteModel.getFrequencyModel();
for (int i = 0; i < m_sequenceLength; i++) {
seq[i] = MathUtils.randomChoicePDF(frequencyModel.getFrequencies());
}
}
if (DEBUG) {
synchronized (this) {
System.out.println();
System.out.println("root Sequence:");
Utils.printArray(seq);
}
}
SimpleAlignment alignment = new SimpleAlignment();
alignment.setReportCountStatistics(false);
alignment.setDataType(m_siteModel.getFrequencyModel().getDataType());
traverse(root, seq, category, alignment);
return alignment;
}
use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class GLMSubstitutionModelParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
DataType dataType = DataTypeUtils.getDataType(xo);
if (dataType == null)
dataType = (DataType) xo.getChild(DataType.class);
int rateCount = (dataType.getStateCount() - 1) * dataType.getStateCount();
LogLinearModel glm = (LogLinearModel) xo.getChild(GeneralizedLinearModel.class);
int length = glm.getXBeta().length;
if (length != rateCount) {
throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (rateCount) + " dimensions. However GLM dimension is " + length);
}
XMLObject cxo = xo.getChild(ComplexSubstitutionModelParser.ROOT_FREQUENCIES);
FrequencyModel rootFreq = (FrequencyModel) cxo.getChild(FrequencyModel.class);
if (dataType != rootFreq.getDataType()) {
throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its rootFrequencyModel.");
}
return new GLMSubstitutionModel(xo.getId(), dataType, rootFreq, glm);
}
use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class GeneralF81ModelParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
FrequencyModel freqModel = null;
if (xo.hasChildNamed(FREQUENCIES)) {
XMLObject cxo = xo.getChild(FREQUENCIES);
freqModel = (FrequencyModel) cxo.getChild(FrequencyModel.class);
}
Logger.getLogger("dr.evomodel").info(" General F81 Model from frequencyModel '" + freqModel.getId() + "' (stateCount=" + freqModel.getFrequencyCount() + ")");
return new GeneralF81Model(freqModel);
}
Aggregations