Mercurial > repos > bcclaywell > argo_navis

annotate bin/format_beastfile.py @ 2:7eaf6f9abd28 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b-dirty
author bcclaywell
date Mon, 12 Oct 2015 17:57:38 -0400
parents d67268158946
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
1 #!/usr/bin/env python
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
2 """
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
3 Formats a given BEAST XML file (possibly all ready to run) and respecifies the information needed to run the
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
4 classic Discrete trait.
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
5
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
6 Some things that would be nice:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
7 * specify output files/formats (could let you run from root instead of the dir)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
8 """
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
9
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
10 from Bio import SeqIO
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
11 import xml.etree.ElementTree as ET
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
12 import argparse
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
13 import copy
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
14 import csv
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
15
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
16
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
17 def clear_children(node):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
18 "Element.remove doesn't seem to work the way it's supposed to, so we're doing this"
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
19 node_attrib = copy.copy(node.attrib)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
20 node.clear()
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
21 node.attrib = node_attrib
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
22
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
23
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
24 def set_alignment(xmldoc, alignment):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
25 """This function replaces the alignment data in xmldoc with that from sequences in alignment."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
26 aln_node = xmldoc.find('data')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
27 # First clear out the old alignment sequences
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
28 clear_children(aln_node)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
29 print "seqs"
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
30 for seq in aln_node:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
31 print seq
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
32 # Next, construct and throw in the new sequence nodes
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
33 for seq_record in alignment:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
34 seqid = seq_record.name
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
35 ET.SubElement(aln_node, 'sequence',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
36 attrib=dict(id="seq_" + seqid,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
37 taxon=seqid,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
38 totalcount="4",
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
39 value=str(seq_record.seq)))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
40
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
41
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
42 def get_data_id(xmldoc):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
43 """The data set will have a given name, assigned by BEAUti, typically based on the named of the data file
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
44 loaded into it. This name gets referred to in a number of places (presumably so there can be a number of
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
45 partitions/datasets in an analysis), and is needed by other bits of code that do their thing."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
46 return xmldoc.find(".//data[@name='alignment'][@id]").attrib['id']
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
47
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
48
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
49 def default_deme_getter(metarow):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
50 """A default function for getting the deme data from a given metadata row. Specifically defaults to 'deme'
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
51 first, then to 'community' next. Returns none if it doesn't find either."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
52 return metarow.get('deme') or metarow.get('community')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
53
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
54
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
55 def set_deme(xmldoc, metadata, deme_getter=default_deme_getter):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
56 """Sets the deme information of the xmldoc based on metadata, and using the deme_getter (by default the
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
57 `default_deme_getter` function above."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
58 trait_node = xmldoc.iter('traitSet').next()
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
59 trait_string = ",\n".join([row['sequence'] + "=" + deme_getter(row) for row in metadata])
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
60 trait_node.text = trait_string
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
61
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
62
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
63 def build_date_node(date_spec, data_id):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
64 """Builds a node of date traits, given the date_spec string which is the actual string representation of
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
65 the sequence -> date mapping. Has to create a `taxa` subnode, and a `data` subnode of that, which points
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
66 to the data set in question via `idref`."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
67 date_node = ET.Element('trait',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
68 id='dateTrait.t:' + data_id,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
69 spec='beast.evolution.tree.TraitSet',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
70 traitname='date')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
71 date_node.text = date_spec
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
72 taxa_node = ET.SubElement(date_node, 'taxa',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
73 id='TaxonSet.' + data_id,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
74 spec='TaxonSet')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
75 _ = ET.SubElement(taxa_node, 'data',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
76 idref=data_id,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
77 name="alignment")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
78 return date_node
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
79
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
80
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
81 def set_date(xmldoc, metadata, date_attr='date'):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
82 """Builds a dateTrait node via `build_date_node` above, and inserts into the `.//state/tree` node.
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
83 However, this `tree` node already contains a `taxonset` node which has a `data` node, and this
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
84 `taxonset` node has the same id as the `taxa` node in the the date `trait` node. As such, the node that
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
85 _was_ present must be removed, so that we don't get a duplicate id error. Instead, we replace the old
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
86 taxonset node with one which has an `idref` pointing to the `taxa` node inside the `trait` node. This is
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
87 rather convoluted, and I'm not possible that some file with multiple datasets wouldn't break on this, but
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
88 this described strategy seems to work for now."""
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
89 # First get our tree node; we'll be adding the date data to this
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
90 tree_node = xmldoc.find('.//state/tree')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
91 # Construct our trait string, just as we do for `set_deme`
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
92 trait_string = ",\n".join([row['sequence'] + "=" + row[date_attr] for row in metadata])
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
93 # Build the date trait node, and carry out all the weird mucking to get the new `taxonset` node in, as
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
94 # described in the docstring
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
95 data_id = get_data_id(xmldoc)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
96 date_node = build_date_node(trait_string, data_id)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
97 old_taxonset = tree_node.find("./taxonset")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
98 tree_node.insert(0, date_node)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
99 tree_node.remove(old_taxonset)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
100 new_taxonset = ET.SubElement(tree_node, "taxonset", idref="TaxonSet."+data_id)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
101
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
102
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
103 def get_current_interval(xmldoc):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
104 run_node = xmldoc.find('run')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
105 loggers = run_node.findall('logger')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
106 intervals = list(set([int(l.get('logEvery')) for l in loggers if l.get('id') != 'screenlog']))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
107 if len(intervals) > 1:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
108 raise "Cannot get an interval for this xml doc; there are multiple such values"
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
109 return intervals[0]
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
110
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
111
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
112 def set_mcmc(xmldoc, samples, sampling_interval):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
113 "Sets the MCMC chain settings (how often to log, how long to run, etc"
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
114 run_node = xmldoc.find('run')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
115 # XXX Should really make it so that you only have to specify _one_, and it will find current value of
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
116 # other so that chain length doesn't break.
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
117 chain_length = samples * sampling_interval + 1
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
118 run_node.set('chainLength', str(chain_length))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
119 loggers = run_node.findall('logger')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
120 for logger in loggers:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
121 logevery = sampling_interval * 10 if logger.get('id') == 'screenlog' else sampling_interval
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
122 logger.set('logEvery', str(logevery))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
123
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
124
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
125 def normalize_filenames(xmldoc, logger_filename="posterior.log", treefile_filename="posterior.trait.trees"):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
126 run_node = xmldoc.find('run')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
127 logfile_node = run_node.find('logger[@id="tracelog"]')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
128 treefile_node = run_node.find('logger[@id="treeWithTraitLogger.deme"]')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
129 logfile_node.set('fileName', logger_filename)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
130 treefile_node.set('fileName', treefile_filename)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
131
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
132
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
133 def set_deme_count(xmldoc, metadata, deme_getter=default_deme_getter):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
134 "Updates the model specs based onthe number of demes in the data set."
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
135 demes = list(set(map(deme_getter, metadata)))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
136 demes.sort()
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
137 deme_count = len(demes)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
138 mig_dim = (deme_count - 1) * deme_count / 2
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
139 for xpath in ['.//parameter[@id="relativeGeoRates.s:deme"]', './/stateNode[@id="rateIndicator.s:deme"]']:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
140 xmldoc.find(xpath).set('dimension', str(mig_dim))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
141 code_map = map(lambda ix: ix[1] + "=" + str(ix[0]), enumerate(demes))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
142 code_map = ",".join(code_map) + ",? = " + " ".join(map(str, range(deme_count))) + " "
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
143 user_data_type_node = xmldoc.find('.//userDataType')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
144 user_data_type_node.set('codeMap', code_map)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
145 user_data_type_node.set('states', str(deme_count))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
146 trait_frequencies_param = xmldoc.find('.//frequencies/parameter[@id="traitfrequencies.s:deme"]')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
147 trait_frequencies_param.set('dimension', str(deme_count))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
148 trait_frequencies_param.text = str(1.0/deme_count)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
149
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
150
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
151
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
152 def get_args():
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
153 def int_or_floatify(string):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
154 return int(float(string))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
155 parser = argparse.ArgumentParser()
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
156 parser.add_argument('template', type=argparse.FileType('r'),
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
157 help="""A template BEAST XML (presumably created by Beauti) ready insertion of an alignment and
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
158 discrete trait.""")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
159 parser.add_argument('-a', '--alignment',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
160 help="Replace alignment in beast file with this alignment; Fasta format.")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
161 parser.add_argument('-m', '--metadata', type=argparse.FileType('r'),
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
162 help="Should contain 'community' column referencing the deme.")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
163 parser.add_argument('-s', '--samples', type=int_or_floatify,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
164 help="Number of samples in output log file(s).")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
165 parser.add_argument('-d', '--deme-col',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
166 help="""Specifies the deme column for metadata; defaults to deme or community (whichever is present)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
167 if not specified.""")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
168 parser.add_argument('-D', '--date-col',
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
169 help="If specified, will add a date specification to the output BEAST XML file.")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
170 parser.add_argument('-i', '--sampling-interval', type=int_or_floatify,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
171 help="""Number of chain states to simulate between successive states samples for logfiles. The
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
172 total chain length is therefor samples * sampling_interval.""")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
173 parser.add_argument('beastfile', type=argparse.FileType('w'),
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
174 help="Output BEAST XML file.")
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
175 return parser.parse_args()
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
176
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
177
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
178 def main(args):
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
179 # Read in old data
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
180 xmldoc = ET.parse(args.template)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
181
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
182 # Modify the data sets
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
183 if args.alignment:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
184 alignment = SeqIO.parse(args.alignment, 'fasta')
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
185 set_alignment(xmldoc, alignment)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
186 if args.metadata:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
187 metadata = list(csv.DictReader(args.metadata))
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
188 # Set the deme getter
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
189 deme_getter = lambda row: row[args.deme_col] if args.deme_col else default_deme_getter(row)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
190 set_deme(xmldoc, metadata, deme_getter)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
191 # _could_ do something smart here where we look at which sequences in the XML file traitset that match
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
192 # alignment passed in if _only_ alignment is passed in. Probably not worth it though...
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
193 set_deme_count(xmldoc, metadata, deme_getter)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
194 if args.date_col:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
195 set_date(xmldoc, metadata, args.date_col)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
196
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
197 if args.samples or args.sampling_interval:
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
198 interval = args.sampling_interval if args.sampling_interval else get_current_interval(xmldoc)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
199 set_mcmc(xmldoc, args.samples, interval)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
200
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
201 # Make sure that we always have the same file names out. These are specified as defaults of the function,
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
202 # but could be customized here or through the cl args if needed.
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
203 normalize_filenames(xmldoc)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
204
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
205 # Write the output
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
206 xmldoc.write(args.beastfile)
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
207
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
208
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
209 if __name__ == '__main__':
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
210 main(get_args())
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
211
d67268158946 planemo upload commit a3f181f5f126803c654b3a66dd4e83a48f7e203b
bcclaywell
parents:
diff changeset
212