Examples with SimpleAlignment dr.evolution.alignment.SimpleAlignment used on opensource projects

Example 11 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class SelectionMapping method main.

public static void main(String[] args) throws java.io.IOException {
    if (args.length != 1)
        throw new RuntimeException("Expect a file name containing parameters!");
    FileInputStream fi = new FileInputStream(args[0]);
    Properties props = new Properties();
    props.load(fi);
    fi.close();
    int genomeLength = Integer.parseInt(props.getProperty("genomeLength", "1000"));
    int populationSize = Integer.parseInt(props.getProperty("populationSize", "250"));
    int replicates = Integer.parseInt(props.getProperty("replicates", "1000"));
    double mu = Double.parseDouble(props.getProperty("mu", "1e-4"));
    double alpha = Double.parseDouble(props.getProperty("alpha", "0.06"));
    double pInv = Double.parseDouble(props.getProperty("pinv", "0.5"));
    double beta = Double.parseDouble(props.getProperty("beta", "1"));
    int stateSize = Integer.parseInt(props.getProperty("stateSize", "4"));
    int n = Integer.parseInt(props.getProperty("sampleSize", "40"));
    boolean randomFittest = Boolean.valueOf(props.getProperty("randomizeFittest", "false"));
    boolean outputAlignments = Boolean.valueOf(props.getProperty("outputAlignments", "false"));
    String alignmentFileName = props.getProperty("alignment.filename", "alignment.fasta");
    String fileName = props.getProperty("output.filename", "out.txt");
    int burninFactor = Integer.parseInt(props.getProperty("burninFactor", "10"));
    PrintWriter writer = new PrintWriter(new FileWriter(fileName));
    writer.println("// WMD v1.0");
    writer.println("// genomeLength: " + genomeLength);
    writer.println("// populationSize: " + populationSize);
    writer.println("// alpha: " + alpha);
    writer.println("// beta: " + beta);
    writer.println("// replicates: " + replicates);
    writer.println("// output.filename: " + fileName);
    writer.println("// mu: " + mu);
    writer.println("// outputAlignments: " + outputAlignments);
    writer.println("// randomFittest: " + randomFittest);
    writer.println("// alignment.filename: " + alignmentFileName);
    int[] totalMutationalDensity = new int[200];
    int[] totalLineages = new int[200];
    int[] nMutants = new int[genomeLength * 3];
    double[] mrcaFitness = new double[1];
    final int generations = populationSize;
    ArrayList[] unfoldedSites = new ArrayList[populationSize + 1];
    for (int i = 0; i < populationSize + 1; i++) {
        unfoldedSites[i] = new ArrayList();
    }
    long start = System.currentTimeMillis();
    int dnaSingles = 0;
    int aaSingles = 0;
    int dnaSegs = 0;
    int aaSegs = 0;
    for (int reps = 0; reps < replicates; reps++) {
        //			FitnessFunction f = null;
        //			if (alpha == 0.0) {
        //				f = new NeutralModel();
        //			}
        FitnessFunction f = new CodonFitnessFunction(genomeLength, alpha, beta, pInv);
        Population p = new Population(populationSize, genomeLength * 3, new SimpleMutator(mu, stateSize), f, randomFittest);
        Genome master = new SimpleGenome(genomeLength * 3, f, randomFittest);
        p = Population.forwardSimulation(p, burninFactor * populationSize);
        int age = -1;
        while (age == -1) {
            p = Population.forwardSimulation(p, generations);
            age = p.getAgeOfMRCA(mrcaFitness);
        }
        Population children = Population.forwardSimulation(p, 1);
        writer.println(reps + "\t" + age + "\t" + p.getMeanParentFitness() + "\t" + mrcaFitness[0] + "\t" + p.getProportionAsFit(mrcaFitness[0]));
        if (reps % 10 == 0) {
            System.out.print(reps + "\t" + age + "\t" + p.getMeanParentFitness() + "\t" + mrcaFitness[0] + "\t" + p.getProportionAsFit(mrcaFitness[0]));
            long millis = System.currentTimeMillis() - start;
            double seconds = millis / 1000.0;
            if (reps != 0) {
                double seconds2go = (replicates - reps) * seconds / reps;
                System.out.println(" -- " + (Math.round(seconds2go / 36) / 100.0) + " hours");
            } else {
                System.out.println();
            }
        }
        //p.unfoldedSiteFrequencies(unfoldedSites);
        // ************************************************************************************
        // get the mutational density per lineage as a function of time
        // ************************************************************************************
        List<Integer>[] mutations = new ArrayList[200];
        for (int i = 0; i < mutations.length; i++) {
            mutations[i] = new ArrayList<Integer>(200);
        }
        int sampleSize = Math.min(n, populationSize);
        p.getMutationDensity(sampleSize, mutations);
        for (int i = 0; i < mutations.length; i++) {
            for (int j = 0; j < mutations[i].size(); j++) {
                totalMutationalDensity[i] += mutations[i].get(j);
            }
            totalLineages[i] += mutations[i].size();
        }
        // ************************************************************************************
        // output alignments if requested
        // ************************************************************************************
        SimpleAlignment alignment = new SimpleAlignment();
        SimpleAlignment aaAlignment = new SimpleAlignment();
        alignment.setDataType(Nucleotides.INSTANCE);
        aaAlignment.setDataType(AminoAcids.INSTANCE);
        if (outputAlignments) {
            PrintWriter pw = new PrintWriter(new FileWriter(alignmentFileName));
            //pw.println(sampleSize + "\t" + p.getGenomeLength());
            for (int i = 0; i < sampleSize; i++) {
                pw.print(">sequence_");
                if (i < 10) {
                    pw.print("0");
                }
                pw.println("" + i);
                String dnaSequence = p.getGenome(i).getDNASequenceString();
                String aaSequence = p.getGenome(i).getAminoAcidSequenceString();
                pw.println(dnaSequence);
                alignment.addSequence(new Sequence(dnaSequence));
                aaAlignment.addSequence(new Sequence(aaSequence));
            }
            pw.println();
            pw.close();
        }
        // ************************************************************************************
        for (int i = 0; i < sampleSize; i++) {
            nMutants[master.hammingDistance(p.getGenome(i))] += 1;
        }
        dnaSingles += countSingletons(alignment);
        aaSingles += countSingletons(aaAlignment);
        dnaSegs += countSegregating(alignment);
        aaSegs += countSegregating(aaAlignment);
    }
    for (int i = 0; i < unfoldedSites.length; i++) {
        double meanFitness = 0.0;
        double proportionNeutral = 0.0;
        double proportionPositive = 0.0;
        double proportionNegative = 0.0;
        for (int j = 0; j < unfoldedSites[i].size(); j++) {
            double relativeFitness = (Double) unfoldedSites[i].get(j);
            meanFitness += relativeFitness;
            if (relativeFitness == 1.0) {
                proportionNeutral += 1.0;
            } else if (relativeFitness > 1.0) {
                proportionPositive += 1.0;
            } else {
                proportionNegative += 1.0;
            }
        }
        meanFitness /= (double) unfoldedSites[i].size();
        proportionNeutral /= (double) unfoldedSites[i].size();
        proportionPositive /= (double) unfoldedSites[i].size();
        proportionNegative /= (double) unfoldedSites[i].size();
        writer.println(i + "\t" + unfoldedSites[i].size() + "\t" + meanFitness + "\t" + proportionNeutral + "\t" + proportionPositive + "\t" + proportionNegative);
    }
    writer.println("--------------------");
    writer.println("SINGLETON COUNTS");
    writer.println("--------------------");
    writer.println("dna singletons = " + dnaSingles);
    writer.println("aa singletons = " + aaSingles);
    int dnaNons = dnaSegs - dnaSingles;
    int aaNons = aaSegs - aaSingles;
    System.out.println("dna singletons = " + dnaSingles);
    System.out.println("aa singletons = " + aaSingles);
    System.out.println("dna segregating = " + dnaSegs);
    System.out.println("aa segregating = " + aaSegs);
    System.out.println("dna non-singles = " + dnaNons);
    System.out.println("aa non-singles = " + aaNons);
    System.out.println("ratio = " + ((double) dnaSingles / (double) aaSingles));
    System.out.println("ratio(non) = " + ((double) dnaNons / (double) aaNons));
    writer.println("--------------------");
    writer.println("MUTATIONAL DENSITIES");
    writer.println("--------------------");
    for (int i = 0; i < 200; i++) {
        writer.println(totalMutationalDensity[i] + "\t" + totalLineages[i]);
    }
    // ************************************************************************************
    // output hamming distances
    // ************************************************************************************
    writer.println("--------------------");
    writer.println("Hamming distance distribution");
    writer.println("--------------------");
    for (int i = 0; i < nMutants.length; i++) {
        writer.println(i + "\t" + nMutants[i]);
    }
    writer.close();
}

Example 12 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class BeagleSequenceSimulator method compileAlignment.

// END: SimulatePartitionCallable class
private SimpleAlignment compileAlignment() {
    SimpleAlignment simpleAlignment = new SimpleAlignment();
    simpleAlignment.setReportCountStatistics(false);
    simpleAlignment.setDataType(dataType);
    LinkedHashMap<Taxon, int[]> alignmentMap = new LinkedHashMap<Taxon, int[]>();
    // compile the alignment
    for (Partition partition : partitions) {
        Map<Taxon, int[]> sequenceMap = partition.getTaxonSequencesMap();
        Iterator<Entry<Taxon, int[]>> iterator = sequenceMap.entrySet().iterator();
        while (iterator.hasNext()) {
            Entry<Taxon, int[]> pairs = (Entry<Taxon, int[]>) iterator.next();
            Taxon taxon = pairs.getKey();
            int[] partitionSequence = pairs.getValue();
            if (alignmentMap.containsKey(taxon)) {
                int j = 0;
                for (int i = partition.from; i <= partition.to; i += partition.every) {
                    alignmentMap.get(taxon)[i] = partitionSequence[j];
                    j++;
                }
            // END: i loop
            } else {
                int[] sequence = new int[siteCount];
                // dirty solution for gaps when taxa between the tree
                // topologies don't match
                Arrays.fill(sequence, gapFlag);
                int j = 0;
                for (int i = partition.from; i <= partition.to; i += partition.every) {
                    sequence[i] = partitionSequence[j];
                    j++;
                }
                // END: i loop
                alignmentMap.put(taxon, sequence);
            }
        // END: key check
        }
    // END: iterate seqMap
    }
    // END: partitions loop
    Iterator<Entry<Taxon, int[]>> iterator = alignmentMap.entrySet().iterator();
    while (iterator.hasNext()) {
        Entry<Taxon, int[]> pairs = (Entry<Taxon, int[]>) iterator.next();
        Taxon taxon = (Taxon) pairs.getKey();
        int[] intSequence = (int[]) pairs.getValue();
        Sequence sequence = //
        Utils.intArray2Sequence(//
        taxon, //
        intSequence, //
        gapFlag, dataType);
        //			sequence.setDataType(dataType);
        simpleAlignment.addSequence(sequence);
        iterator.remove();
    }
    return simpleAlignment;
}

Example 13 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class DataModelImporter method importNexusFile.

// nexus
private void importNexusFile(File file, DateGuesser guesser, Map dataModel) throws IOException, ImportException {
    TaxonList taxa = null;
    SimpleAlignment alignment = null;
    List<Tree> trees = new ArrayList<Tree>();
    List<NexusApplicationImporter.CharSet> charSets = new ArrayList<NexusApplicationImporter.CharSet>();
    try {
        FileReader reader = new FileReader(file);
        NexusApplicationImporter importer = new NexusApplicationImporter(reader);
        boolean done = false;
        while (!done) {
            try {
                NexusBlock block = importer.findNextBlock();
                if (block == NexusImporter.TAXA_BLOCK) {
                    if (taxa != null) {
                        throw new MissingBlockException("TAXA block already defined");
                    }
                    taxa = importer.parseTaxaBlock();
                    dataModel.put("taxa", createTaxonList(taxa));
                } else if (block == NexusImporter.CALIBRATION_BLOCK) {
                    if (taxa == null) {
                        throw new MissingBlockException("TAXA or DATA block must be defined before a CALIBRATION block");
                    }
                    importer.parseCalibrationBlock(taxa);
                } else if (block == NexusImporter.CHARACTERS_BLOCK) {
                    if (taxa == null) {
                        throw new MissingBlockException("TAXA block must be defined before a CHARACTERS block");
                    }
                    if (alignment != null) {
                        throw new MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    alignment = (SimpleAlignment) importer.parseCharactersBlock(taxa);
                } else if (block == NexusImporter.DATA_BLOCK) {
                    if (alignment != null) {
                        throw new MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    // A data block doesn't need a taxon block before it
                    // but if one exists then it will use it.
                    alignment = (SimpleAlignment) importer.parseDataBlock(taxa);
                    if (taxa == null) {
                        taxa = alignment;
                    }
                } else if (block == NexusImporter.TREES_BLOCK) {
                    // I guess there is no reason not to allow multiple trees blocks
                    //                        if (trees.size() > 0) {
                    //                            throw new MissingBlockException("TREES block already defined");
                    //                        }
                    Tree[] treeArray = importer.parseTreesBlock(taxa);
                    trees.addAll(Arrays.asList(treeArray));
                    if (taxa == null && trees.size() > 0) {
                        taxa = trees.get(0);
                    }
                } else if (block == NexusApplicationImporter.ASSUMPTIONS_BLOCK) {
                    importer.parseAssumptionsBlock(charSets);
                } else {
                // Ignore the block..
                }
            } catch (EOFException ex) {
                done = true;
            }
        }
        reader.close();
        // Allow the user to load taxa only (perhaps from a tree file) so that they can sample from a prior...
        if (alignment == null && taxa == null) {
            throw new MissingBlockException("TAXON, DATA or CHARACTERS block is missing");
        }
    } catch (IOException e) {
        throw new IOException(e.getMessage());
    } catch (ImportException e) {
        throw new ImportException(e.getMessage());
    //        } catch (Exception e) {
    //            throw new Exception(e.getMessage());
    }
    setData(dataModel, guesser, file.getName(), taxa, null, alignment, charSets, null, trees);
}

Example 14 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class BeagleSequenceSimulatorParser method parseXMLObject.

// END: getSyntaxRules
@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    String msg = "";
    boolean parallel = false;
    boolean outputAncestralSequences = false;
    if (xo.hasAttribute(PARALLEL)) {
        parallel = xo.getBooleanAttribute(PARALLEL);
    }
    if (xo.hasAttribute(OUTPUT_ANCESTRAL_SEQUENCES)) {
        outputAncestralSequences = xo.getBooleanAttribute(OUTPUT_ANCESTRAL_SEQUENCES);
    }
    SimpleAlignment.OutputType output = SimpleAlignment.OutputType.FASTA;
    if (xo.hasAttribute(OUTPUT)) {
        output = SimpleAlignment.OutputType.parseFromString(xo.getStringAttribute(OUTPUT));
    }
    int siteCount = 0;
    int to = 0;
    for (int i = 0; i < xo.getChildCount(); i++) {
        Partition partition = (Partition) xo.getChild(i);
        to = partition.to + 1;
        if (to > siteCount) {
            siteCount = to;
        }
    }
    // END: partitions loop
    ArrayList<Partition> partitionsList = new ArrayList<Partition>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Partition partition = (Partition) xo.getChild(i);
        if (partition.from > siteCount) {
            throw new XMLParseException("Illegal 'from' attribute in " + PartitionParser.PARTITION + " element");
        }
        if (partition.to > siteCount) {
            throw new XMLParseException("Illegal 'to' attribute in " + PartitionParser.PARTITION + " element");
        }
        if (partition.to == -1) {
            partition.to = siteCount - 1;
        }
        if (partition.getRootSequence() != null) {
            //            	TODO: what about 'every'?
            int partitionSiteCount = (partition.to - partition.from) + 1;
            if (partition.getRootSequence().getLength() != 3 * partitionSiteCount && partition.getFreqModel().getDataType() instanceof Codons) {
                throw new RuntimeException("Root codon sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + 3 * partitionSiteCount + " characters");
            } else if (partition.getRootSequence().getLength() != partitionSiteCount && partition.getFreqModel().getDataType() instanceof Nucleotides) {
                throw new RuntimeException("Root nuleotide sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + partitionSiteCount + " characters");
            }
        // END: dataType check
        //                System.exit(-1);
        }
        // END: ancestralSequence check
        partitionsList.add(partition);
    }
    // END: partitions loop
    msg += "\n\t" + siteCount + ((siteCount > 1) ? " replications " : " replication");
    if (msg.length() > 0) {
        Logger.getLogger("dr.app.beagle.tools").info("Using Beagle Sequence Simulator: " + msg);
    }
    BeagleSequenceSimulator s = new BeagleSequenceSimulator(partitionsList);
    SimpleAlignment alignment = s.simulate(parallel, outputAncestralSequences);
    alignment.setOutputType(output);
    return alignment;
}

Example 15 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class BEAUTiImporter method importNexusFile.

// nexus
private void importNexusFile(File file) throws IOException, ImportException {
    TaxonList taxa = null;
    SimpleAlignment alignment = null;
    List<Tree> trees = new ArrayList<Tree>();
    PartitionSubstitutionModel model = null;
    List<NexusApplicationImporter.CharSet> charSets = new ArrayList<NexusApplicationImporter.CharSet>();
    try {
        FileReader reader = new FileReader(file);
        NexusApplicationImporter importer = new NexusApplicationImporter(reader);
        boolean done = false;
        while (!done) {
            try {
                NexusBlock block = importer.findNextBlock();
                if (block == NexusImporter.TAXA_BLOCK) {
                    if (taxa != null) {
                        throw new MissingBlockException("TAXA block already defined");
                    }
                    taxa = importer.parseTaxaBlock();
                } else if (block == NexusImporter.CALIBRATION_BLOCK) {
                    if (taxa == null) {
                        throw new MissingBlockException("TAXA or DATA block must be defined before a CALIBRATION block");
                    }
                    importer.parseCalibrationBlock(taxa);
                } else if (block == NexusImporter.CHARACTERS_BLOCK) {
                    if (taxa == null) {
                        throw new MissingBlockException("TAXA block must be defined before a CHARACTERS block");
                    }
                    if (alignment != null) {
                        throw new MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    alignment = (SimpleAlignment) importer.parseCharactersBlock(taxa);
                } else if (block == NexusImporter.DATA_BLOCK) {
                    if (alignment != null) {
                        throw new MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    // A data block doesn't need a taxon block before it
                    // but if one exists then it will use it.
                    alignment = (SimpleAlignment) importer.parseDataBlock(taxa);
                    if (taxa == null) {
                        taxa = alignment;
                    }
                } else if (block == NexusImporter.TREES_BLOCK) {
                    // I guess there is no reason not to allow multiple trees blocks
                    //                        if (trees.size() > 0) {
                    //                            throw new MissingBlockException("TREES block already defined");
                    //                        }
                    Tree[] treeArray = importer.parseTreesBlock(taxa);
                    trees.addAll(Arrays.asList(treeArray));
                    if (taxa == null && trees.size() > 0) {
                        taxa = trees.get(0);
                    }
                } else if (block == NexusApplicationImporter.PAUP_BLOCK) {
                    model = importer.parsePAUPBlock(options, charSets);
                } else if (block == NexusApplicationImporter.MRBAYES_BLOCK) {
                    model = importer.parseMrBayesBlock(options, charSets);
                } else if (block == NexusApplicationImporter.ASSUMPTIONS_BLOCK || block == NexusApplicationImporter.SETS_BLOCK) {
                    importer.parseAssumptionsBlock(charSets);
                } else {
                // Ignore the block..
                }
            } catch (EOFException ex) {
                done = true;
            }
        }
        reader.close();
        // Allow the user to load taxa only (perhaps from a tree file) so that they can sample from a prior...
        if (alignment == null && taxa == null) {
            throw new MissingBlockException("TAXON, DATA or CHARACTERS block is missing");
        }
    } catch (IOException e) {
        throw new IOException(e.getMessage());
    } catch (ImportException e) {
        throw new ImportException(e.getMessage());
    //        } catch (Exception e) {
    //            throw new Exception(e.getMessage());
    }
    setData(file.getName(), taxa, alignment, charSets, model, null, trees);
}