You will probably not need to directly look at or edit LAMARC's XML input files, because these are generated for you by our file conversion programs. However, in rare cases, you may need or want to write or edit these files by hand. Here is the necessary information for doing so.
We are not currently supporting all the nice abilities that would really justify using XML, but we hope to do so in the future.
Please be careful not to save the XML file in Microsoft Word or another word processor format. This will fill it with invisible editing characters, and it will not work. Always use the "save as plain text" option of your word processor. If a previously successful input file stops working, check whether you have accidentally saved it as formatted text. This is the #1 cause of input-file problems in both PHYLIP and LAMARC today.
Also, be aware that some word processors (including at least one version of MS Word) will automatically parse XML text, reducing your input data file to a disorganized heap of data. If this happens, you can try giving your file an extension of .txt. If that doesn't work, we have no suggestions other than using a different word processor.
XML is a very simple minded language, so you have to be methodical as it is easily confused. Best approach is to start with a functioning script, either the one generated for you by LAMARC or one of examples included in this web site. (If you find one that doesn't work, please tell us. We've tested them, but they're easily broken.) Starting from an appropriate script, incrementally modify it to your needs, testing as you go.
Looking at our examples, you'll note that each level of tag is indented 4 spaces further than the next higher level tag. This not required by the XML compiler (it ignores all those spaces) but is invaluable when trying to debug a malfunctioning script. We highly recommend you follow this convention. It will save you many hours of pain (don't ask how we know that :-). One thing to particularly avoid is XML that is indented differently from its meaning. This leads to exceptionally hard-to-find errors.
If you are unfortunate enough to start with an XML script that is not properly indented, it is worth your time to go through and indent it properly. Not only will you understand it better, the indenting it will make modifications much easier.
Each bit of information in the input file is surrounded by a beginning and ending tag. For example, the whole input file is surrounded by <lamarc> and </lamarc>, and the data section is surrounded by <data> and </data>. Tags must come in <tag> </tag> pairs and must be strictly nested, so these are legal:
<blue> <red> </red> </blue>
<red> </red> <blue> </blue>
but this is illegal:
<blue> <red> </blue> </red>
The symbol "<!--" starts a comment, and the symbol "-->" ends it. Anything inside a comment is ignored. You can use comment symbols to temporarily remove parts of your input file that you don't want, as long as what remains is legal.
Some tags require additional information, such as the name of a population. This is done with "attributes":
<population name="Washington DC">
The quotes are required.
The data section contains your actual molecular data, and additional information used to interpret it. It is required to be present, and is enclosed in <data> tags.
LAMARC divides molecular data into "regions". A region is all the available genetic information that is closely linked on the same chromosome and has a known map. Use a single region for data which is one contiguous stretch, so that it would be meaningful to calculate recombination rates along it. Use multiple regions for data composed of several disconnected bits or bits whose connections are not known. Regions are enclosed in the <region> tag, and at least one must be present. The region's name can be added in an optional name attribute.
The optional tag <effective-popsize> may be used to specify a different relative effective population size for each <region>. For example, data from nuclear chromosomes of a diploid organism reflect an effective population size four times larger than data from the mitochondrion of the same organism. Data from sex chromosomes also have unique effective population sizes--the relative effective population size ratios for a non-sex chromosome to an X chromosome to a Y chromosome is 4:3:1. Be aware that the parameter estimates produced by LAMARC will be scaled proportional to an effective population size of 1. That is, if you tell LAMARC that you have two regions, one with an effective population size of 4 and one with an effective population size of 3, your final overall estimate of Θ will be lowered to correspond to an effective population size of 1. If you are combining mitochondrial and autosomal data in diploids into a joint analysis, set your relative effective population sizes to 0.25 and 1, respectively, if you want the joint Θ estimate to be reported using the autosomal scale; otherwise, set these to 1 and 4, respectively, to obtain Θ using the mitochondrial scale.
Within a region there may be several "contiguous segments." Segments are stretches of genetic information which are linked to one another but need separate handling. For instance, if we have a stretch of DNA (modelled by a nucleotide substitution model) and an adjacent microsatellite (modelled by a stepwise model) they need to be in the same region as they are linked, but cannot be in the same segment as they require different mutational models. Information about the relative position of segments is placed in a <spacing> tag.
Each segment is indicated by a <block> tag, which can also give information about the position of the segment itself, and (particularly for non-contiguous markers such as SNPs and microsatellites) the positions of the markers within the segment.
Within <block>, the <length> tag indicates the total length of the segment. This is important for SNPs in particular because correct interpretation of the markers requires knowledge of how many non-markers were surveyed. The <map-position> tag gives the position of this segment on an overall map of the region. The map position is the point at which sequencing or scanning began, even if that is not the zero point of your segment's internal numbering. For example, if you sequenced a gene from upstream position -45 to downstream position 500, your map position should tell where position -45 is on the overall map of the region. The <locations> tag encloses a list of marker positions within the region; for example, the tag <locations> -10 7 18 22 </locations> indicates that the four markers are at positions -10, 7, 18 and 22 with respect to the segment's numbering system. The <offset> tag gives the origin of the segment's numbering system with respect to the boundaries of the region. For example, if you began sequencing at position -45 with regard to your chosen numbering system (perhaps you have begun your numbering at a gene's start codon but obtained upstream sequence), your <offset> would be -45. If no offset is given, the offset is assumed to be zero.
Within each region you can list various populations. If you list <population> tags under more than one region, they will be matched by means of their name attributes, so the names are not optional. Use of good names here will also make your output much easier to interpret!
If a SNP panel was used to generate your input data, the panel size will be specified under the <population> tags. Panel members are represented as additional tips in the coalescence tree that have unknown data. All that is needed is the number of tips, though optionally one can name the panel. Here is a fragment that shows the panel correction information for one of two populations.
<lamarc>
...
<data>
<region name="Alcohol dehydrogenase">
<population name="Seattle">
<panel name="adhpanel">
<panel-size> 6 </panel-size>
</panel>
<individual> ... </individual>
...
</population>
<population>
...
</population>
</region>
</data>
Within each population you can list various individuals. An individual represents all the data for that region that comes from a single biological individual, which may be one, two, or more sets of data depending on how you obtained your data. For example, an individual might consist of one mtDNA sequence, or two nuclear DNA sequences. Individuals can have a name attribute, but it is optional; the <individual> tag itself is required.
Within each individual, you can have one or more <sample> tags indicating the actual sequences. For example, an individual with one mtDNA sequence would have a single <sample> tag for it, which would contain this sequence. An individual with two nuclear DNA sequences would contain two <sample> tags. While the two sequences are treated as separate tips of the tree, their identity as a single individual is important if haplotypes are to be considered.
Optionally, an individual can have a <phase> tag, indicating uncertainty about the phase of certain sites. For example, if you have only genotypic data, you will want to indicate uncertainty about the phase of all sites. The <phase> tag has an obligatory attribute, "type," which can be either "known" or "unknown." If the type is "known," the list of sites which follows is the list of all sites whose phase is known, and therefore need not be reconsidered during the run. If the type is "unknown," the list which follows is the list of all sites whose phase is unknown, and thus should be reconsidered. A simple way to code a sequence of genotypic data with no phase information is <phase type="known"> </phase>. It is not necessary to specify homozygous sites as to their phase-known or phase-unknown status, as they will not be reconsidered anyway.
Valid values for the phase tag are site numbers between the value of the offset for that segment (which defaults to 1) and the length of the segment plus the offset. If the segment is longer than the number of markers you have (as is the case for SNP data), valid values here are the same values used for the 'locations' tag in the 'block' section (above). (Note that versions of LAMARC prior to 2.1 did not use this numbering scheme, and instead required you to indicate the particular marker number, starting from 0, that was phased or unphased. Old LAMARC infiles will give errors or fail to run properly using the 2.1 system.)
The sequences themselves are enclosed in <datablock> tags, one per segment per sample. Currently we support DNA, RNA, SNP, and microsatellite data. Each datablock must have an attribute indicating the type of data it contains. Use type="DNA" for full DNA or RNA sequences, type="SNP" for SNP sequences, and type="Microsat" for microsatellites. Please do not mislabel SNPs as full DNA, because the estimates of population size will become vastly overblown.
Sequence data must be aligned and of the same length for all samples within a region. "Unknown nucleotide" codes (X, N or -) can be used to fill in missing or unknown sequence. There is no point in including individuals for whom the entire sequence is unknown, as they add nothing to the analysis (and will slow it down). The full IUPAC nucleotide ambiguity code is available, and DNA and RNA are both accepted (and treated identically). Upper- and lowercase nucleotide symbols are treated equivalently. Deletions should be coded as unknown, and will be treated as unknown; no attempt is made to model the insertion/deletion process.
Here is a minimal DNA data block describing a single region, a single segment, a single population, and two individuals with a single haplotype each. Note that while the two blocks of data are differently formatted, they contain the same number of bases; this is required since all blocks corresponding to a single segment must contain the same number of markers. If your sequences for a given segment are of different lengths, they must be padded out with unknown-nucleotide codes.
<data>
<region name="Alcohol dehydrogenase">
<population name="Seattle">
<individual name="Mary">
<sample>
<datablock type="DNA">
CTTGTAACCTAATGGCTTCCGAGATGGACTAGTGAGCCGCTTTCTC
TACACCAACGCAGCACATGACGGTCTTACATGCGGAGCCCGCTCAA
</datablock>
</sample>
</individual>
<individual name="Jon">
<sample>
<datablock type="DNA">
CTTGTAACCTAATGGCTTCCGA
GATGGACTAGTGAGCCGCTTTCTC
TACACCAACGCAGCACATGACG
GTCTTACATGCGGAGCCCGCTCAA
</datablock>
</sample>
</individual>
</population>
</region>
</data>
Microsatellite data are coded as the number of repeats, with "?" standing for unknown data. Successive microsatellites within the same region are separated by blank spaces. Here is a microsatellite data block which also illustrates the use of multiple samples per individual. In this example, "Mary" is a heterozygote for the second microsatellite and a homozygote for the other five.
<data>
<region name="Alcohol dehydrogenase">
<population name="Seattle">
<individual name="Mary">
<sample>
<datablock type="Microsat">
7 8 14 7 9 21
</datablock>
</sample>
<sample>
<datablock type="Microsat">
7 9 14 7 9 21
</datablock>
</sample>
</individual>
<individual name="Jon">
<sample>
<datablock type="Microsat">
7 9 14 7 10 23
</datablock>
</sample>
<sample>
<datablock type="Microsat">
8 9 13 7 ? 23
</datablock>
</sample>
</individual>
</population>
</region>
</data>
A <region> can contain models <model> to be used to interpret its data. One model should be provided per segment within the region, in the same order that the segments are listed. You can also provide a global <model> (inside the <lamarc> tag) which will be used for any region that does not provide its own models. Only models in these two locations will be used by the program; if both are present, the regional models will be used in preference to the global one. The model chosen must be appropriate for the data type chosen. If no model is provided, then a default model of the appropriate datatype will be used by the program--nucleotides use the Felsenstein '84 model, and microsatellites use the Brownian motion model.
Several tags are common to all models.
The <relative-murate> tag is used to specify segment-specific mutation rates. This is essential if very different data types are combined in one analysis. For example, including DNA and microsat data without taking account of the fact that microsats mutate thousands of times more rapidly than the single-base substitution rate will lead to a nonsense overall estimate. If all of your segments have the same expected mutation rate, this tag is not needed. When it is needed, you should chose some data type as the standard and set its relative rate to 1, and give all other data types in relation to that standard. The estimated parameters will then also be in terms of your standard.
Be aware that the parameter estimates produced by LAMARC will be scaled proportional to a segment of relative mutation rate 1, even if no such segment is included in the data. That is, if you tell LAMARC that you have two segments, one with a relative mutation rate of 5 and the other with a relative mutation rate of 50, your final estimate of Θ will describe a fictional segment with a relative mutation rate of 1, and you will need to multiply by 5 or 50 to find the Θ of your actual segments.
If you believe that your segments vary in mutation rate according to a gamma distribution, you may wish to use the gamma-estimation facilities of LAMARC instead of the <relative-murate> tag. Relative mutation rate is designed for the case where different segments or regions fall into clearly distinct groups, known in advance, such as DNA versus microsats or introns versus exons.
The <categories> tag encloses information about variable rate categories per site (or microsat). Within it, <num-categories> gives the number of rate categories, <rates> gives the relative mutation rate of each category, and <probabilities> gives the probability that a site is in each category. The <probabilities> must add up to 1.0, and there must be a rate and a probability for each category. The <autocorrelation> tag encloses information about the average length of a run of sites (or microsats) with the same rate. For example, if you believe your data consists of runs averaging 100 sites with the same rate, you would set it to 100.0. If you wish to assume that there is no autocorrelation of rates, set this value to 1.0. (We expect 1.0 to be an appropriate value for microsatellites in most cases.)
The <normalize> tag
controls internal normalization of data likelihood values. Data likelihoods can be extremely small, and may be subject to underflow (that is, rounded down to zero, resulting in a loss of information). Normalization attempts to buffer them against underflow at the cost of making the program run more slowly. Versions 1.2 and above of LAMARC turn on normalization automatically if they encounter underflow, so this tag should almost always be set to "false," and we may remove it in a future version. If you are sure your data set will need normalization (extremely polymorphic data, especially microsatellites, or extremely large data sets), then there is no harm in turning it on to free the program from the task of diagnosing the need for normalization.The remaining tags are specific to the particular data model. The current selection of models is described below.
Lamarc offers a choice of two mutational models for DNA, RNA, and SNP data. The simpler (and slightly quicker) model is the model of Felsenstein 1984 (F84), which allows differing nucleotide frequencies and differing rates of transition versus transversion. The more general model is the General Time-Reversible model (GTR), which allows every pair of nucleotides to have a characteristic mutation rate, as well as allowing differing nucleotide frequencies.
Simpler models such as Jukes-Cantor or Kimura Two-Parameter can be obtained by correct choice of parameters to the two models given. If a model can be expressed as a simplification of F84, it will run faster than the same model expressed as a simplification of GTR.
The difference between F84 and GTR is that F84 expresses its mutation rate information as a single parameter, the ttratio (ratio of transitions to transversions) whereas GTR uses six parameters (relative mutation rates between each pair of bases).
Here is a sample F84 data model, which could appear either globally or within a region.
<model name="F84">
<base-freqs> 0.25 0.25 0.25 0.25 </base-freqs>
<per-base-error-rate> 0.001 </per-base-error-rate>
<ttratio> 2.0 </ttratio>
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 5.0 </autocorrelation>
</categories>
<normalize> false </normalize>
</model>
To change this to a GTR model, remove the ttratio, change the model name, and add a line for the mutational rates. We recommend using an external tool such as PAUP* with Modeltest to estimate appropriate GTR rates.
<model name="GTR">
<base-freqs> 0.25 0.25 0.25 0.25 </base-freqs>
<per-base-error-rate> 0.001 </per-base-error-rate>
<gtr-rates> 5.4050 147.7765 4.2745 3.5801 96.3678 1.0 </gtr-rates>
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 5.0 </autocorrelation>
</categories>
<normalize> false </normalize>
</model>
The <base-freqs> tag sets the frequencies of the nucleotides A, C, G, and T (or U) in that order. They must be strictly greater than zero, and must add to 1. Alternatively, instead of four frequencies, you can enter the keyword "calculated," which will cause the program to calculate nucleotide frequencies based on your input data. This will not work if one or more nucleotides are missing from your input data; you must explicitly set four non-zero frequencies in such cases. We also recommend that you always set the base frequencies explicitly when using the GTR model.
The <per-base-error-rate> tag gives the rate at which each individual nucleotide should be assumed to have been miss-called. A value of 0 indicates that all were sequenced correctly. A value of 0.001 indicates one in one thousand is incorrect. If not present, the value is assumed to be 0. This functionality is in beta test as of December, 2009.
The <ttratio> tag (F84 model only) gives the ratio of transitions to transversions. The Jukes-Cantor model (no transition bias) would correspond to a <ttratio> of 0.5, but due to a limitation in the algorithm this (and lower) values are illegal. If you wish to use a Jukes-Cantor model, set <ttratio> to a very slightly larger number such as 0.500001.
The <gtr-rates> tag (GTR model only) gives the relative mutation rate of each pair of nucleotides, in the order AC, AG, AT, CG, CT, GT. Only relative values are important. These are rates as output by PAUP*, before consideration of nucleotide frequencies.
For microsatellite data we offer four models: the stepwise model of Beerli and Felsenstein (1999), the Brownian model of Beerli and Felsenstein (in preparation), a simple K-Allele model similar in concept to the Jukes-Cantor DNA model but for arbitrary K>1, and a mixed K-Allele/Stepwise model where both stepwise and K-Allele-type mutations are allowed at a relative ratio.
The stepwise model assumes that microsatellites evolve via stepwise changes and are constrained not to go below one repeat. This model currently has no unique user-settable parameters; it deduces its required number of bins from the data, and always considers a window of 10 steps on either side of the most extreme data elements unless this is found to overlap zero.
Here is a sample data model which could appear either globally or within a region:
<model name="Stepwise">
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 1.0 </autocorrelation>
</categories>
<normalize> false </normalize>
</model>
The Brownian model assumes that the changes in microsatellite length can be approximated by a continuous distribution (we use a Normal). This model currently has no unique user-settable parameters. It is much faster than the stepwise model, and appears to work well, except for genealogies with very short branches (such as those associated with very small population sizes) on which it shows a significant upward bias. When using this model, be on the lookout for data log-likelihoods of zero in the runtime reports (these are labelled "Data lnL"). If many of these appear, they are an indication that your population sizes are too small for safe use of the Brownian approximation.
<model name="Brownian">
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 1.0 </autocorrelation>
</categories>
<normalize> false </normalize>
</model>
The K-Allele model assumes that the observed alleles represent all possible alleles and that mutation is equally likely among any pair of alleles. It is probably not appropriate for most microsatellite data, and is provided mainly for its ability to handle data types otherwise not analyzable with Lamarc, such as elecrophoretic or indel data. The K-Allele model can also be used to assess how severe an effect violation of the Stepwise model's assumptions might have on the results, since it is essentially the opposite of the Stepwise model. If both have similar results, the results are probably quite insensitive to the details of the microsatellite mutational process.
A K-Allele model block looks just like the Brownian one with a different name:
<model name="KAllele">
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 1.0 </autocorrelation>
</categories>
<normalize> false </normalize>
</model>
The Mixed K-Allele/Stepwise considers both Stepwise mutational possibilities and K-Allele mutational possibilities. The relative weight of the two types of mutation is given by the parameter 'percent_stepwise'. The allele range considered for the K-Allele changes is the same as that for the Stepwise model (that is, all alleles within a certain range of any observed allele) and will therefore include some alleles never observed in the data. The initial percent_stepwise is set by the user, and if the 'optimize' option is set, it is reset at the end of every chain to its optimum value based on the final genealogy of that chain, as calculated using the bisection approach.
Because this model incorporates both a K-Allele and Stepwise approach, it should only be used for data for which both of those models are legal, namely microsatellite data.
A Mixed K-Allele/Stepwise model block includes two new tags, 'alpha' (meaning the percent_stepwise parameter), and 'optimize', which should be on the same level as 'categories' and 'normalize'.
<model name="MixedKS">
<categories>
<num-categories> 2 </num-categories>
<rates> 1.0 10.0 </rates>
<probabilities> 0.8 0.2 </probabilities>
<autocorrelation> 1.0 </autocorrelation>
</categories>
<normalize> false </normalize>
<alpha> 0.3 </alpha>
<optimize> false </optimize>
</model>
The <forces> tag encloses information about each evolutionary force to be considered in the analysis. The <forces> section should not specify forces that make no sense--for example, migration is not allowed if there is only one population, and recombination is not allowed if there is only one site.
The force tags are <coalescence>, <migration>, <recombination>, <growth>, <gamma-over-regions> and the pair <divergence-migration> and <divergence>. In this version of Lamarc. <coalescence> is required, the others are optional (though if multiple populations are specified, either <migration> or the <divergence-migration>, <divergence> pair is required).
The gamma-over-regions "force" may only be applied to data spread over multiple, unlinked genomic regions. This "force" assumes the relative mutation rates over unlinked genomic regions are gamma-distributed, and causes Lamarc to simultaneously infer the shape of the gamma distribution which best fits the data. More information about this can be found here. Please note that Lamarc is unable to co-estimate the shape of the gamma distribution and population growth rates. Lamarc will accept input files containing either the tag <growth> or the tag <gamma-over-regions>, but not both.
Divergence can only be defined if there are two or more populations in the data. If the user wishes to estimate Divergence, they must define the relationships between the populations in the input data by defining parents. This is most easily done in the data file converter, but can be done directly in the XML as shown below.
For each force tag, the following information can be included:
<start-values> contains a space-delimited list of the starting values of the parameters for that force. For example, in a 3-population case, you would provide 3 starting values for population Thetas, and 9 starting values for immigration rates. In the migration case, diagonal entries (meaningless values for migration from a population to itself) can be indicated with dashes instead of zeros, if desired.
<method> indicates the algorithms for computing starting values for each parameter. This can be "User," meaning that the user-specified values should be used. Other options are "FST" to set migration parameters using the FST algorithm, and "Watterson" to set Thetas using the method of Watterson, We do not currently provide algorithms to estimate starting recombination or growth rates, or a starting value for the single parameter of the gamma-over-regions "force." (None of the available algorithms seemed to perform well enough to be helpful.) If you are going to specify any methods, you need to specify one for each parameter.
It is sometimes helpful to set certain parameters to "User" even if you mainly intend to use FST and Watterson. FST, in particular, will fail in certain cases. We use an arbitrary default value when FST fails, but you may be able to provide a better value than this default.
<max-events> gives a maximum for the number of events in the tree that are generated by this force. For example, for the migration force it gives the maximum total number of migrations. If a tree violates any of these maxima, that tree will be discarded. Discarded trees are noted in the runtime reports and, if verbose output is requested, in the output file. This option is not useful (though it is harmless) for coalescence, growth, and gamma-over-regions.
You may wish to set the maxima relatively low if the program is running out of space or slowing down tremendously. However, if the maxima are encountered often (look at the runtime report to check), the estimate of parameters for that force will tend to be biased downward.
<profiles> gives the type of profile likelihood to be computed for each parameter. The options are "percentile," which will compute profile likelihoods at selected percentiles of the distribution; "fixed," which will compute profile likelihoods at fixed multiples of the maximum-likelihood parameter value; and "none," which will compute no profiles. The profiles of a given run may be a mix of "percentile," "fixed," and "none," but you cannot use both "percentile" and "fixed" percentiles for the same force.
In a likelihood analysis, percentile profiles are very time-consuming. If you are not interested in a particular parameter, consider turning off its profiling, and if the overall run is still too slow, consider fixed rather than percentile profiles. If you are only interested in the 95% support interval, you may also change the output file verbosity to 'concise', which causes only those intervals to be calculated and output (decreasing the time spent profiling by approximately a factor of 5). In a Bayesian analysis, profiling is simply a matter of reading off values from the produced posterior probability curve, so the more-informative percentile profiles should be the profile type of choice for most users.
The profile line, if present, must have one entry per parameter for that force (for example, a three-population case must have nine entries for migration profiling). It does not matter what profile type you specify for the "diagonal" migration rates.
An example <forces> block for a case with two populations and migration:
<forces>
<coalescence>
<start-values> 0.01 0.03 </start-values>
<method> Watterson Watterson </method>
<profiles> fixed fixed </profiles>
</coalescence>
<migration>
<start-values> - 1.2 1.8 - </start-values>
<method> - FST User User - </method>
<max-events> 1000 </max-events>
<profiles> - none fixed - </profiles>
</migration>
<recombination>
<start-values> 0.04 </start-values>
<method> User </method>
<max-events> 1000 </max-events>
<profiles> none </profiles>
</recombination>
</forces>
Divergence introduces population divergence into migration, so things get much more complex. The populations are combined pairwise, so if there are three measured populations, there will be two ancestors. Once the first ancestor is defined, migration cannot happen between the two populations subsumed into that ancestor and the third population, as those two populations no longer exist. This is shown graphically on the Divergence page. As a result you need to define two forces <divergence> and <divergence-migration>.
Besides all the standard force tags <divergence> also specifies how the measured populations and their ancestors are related. This is done creating the <population-tree> which explicitly defines, using <epoch-boundary> tags, a pair of <new-populations> and their <ancestor>. That <ancestor> name can then be used as a <new-populations> member for an earlier <ancestor>.
<divergence-migration> defines the migrations between the populations and ancestors (it has exactly the same format as <migration> but is, of course, much larger because of the additional rows and columns required to specify migration to and from the ancestor populations). Below is the XML for the 3 population example on the Divergence page.
Note that in the example both <divergence> and <divergence-migration> have a <prior> defined. This is because Divergence only functions in Bayesian analysis currently. Likelihood analysis is not available.
<divergence>
<prior type="linear">
<paramindex> default </paramindex>
<lower> 0.0 </lower>
<upper> 0.01 </upper>
</prior>
<method> USER USER </method>
<start-values> 0.002000 0.004000 </start-values>
<population-tree>
<epoch-boundary>
<new-populations> North South </new-populations>
<ancestor> Parent_1 </ancestor>
</epoch-boundary>
<epoch-boundary>
<new-populations> East Parent_1 </new-populations>
<ancestor> Parent_2 </ancestor>
</epoch-boundary>
</population-tree>
</divergence>
<divergence-migration>
<start-values> 0 50.000000 50.000000 0 0 50.000000 0 50.000000 0 0
50.000000 50.000000 0 50.000000 0 0 0 50.000000 0 0 0 0 0 0 0 </start-values>
<method> USER USER USER USER USER USER USER USER USER USER USER
USER USER USER USER USER USER USER USER USER USER USER USER USER USER </method>
<max-events> 10000 </max-events>
<profiles> None None None None None None None None None None None
None None None None None None None None None None None None None None </profiles>
<constraints> Invalid Unconstrained Unconstrained Invalid Invalid
Unconstrained Invalid Unconstrained Invalid Invalid Unconstrained Unconstrained
Invalid Unconstrained Invalid Invalid Invalid Unconstrained Invalid Invalid
Invalid Invalid Invalid Invalid Invalid </constraints>
<prior type="linear">
<paramindex> default </paramindex>
<lower> 0.0 </lower>
<upper> 100.0 </upper>
</prior>
</divergence-migration>
The <chains> tag contains information controlling the search strategy, such as number and length of chains, sampling interval, heating, and rearrangement strategy. For details on selecting a search strategy, see the file "Search Strategies."
Chains come in two kinds, "initial" and "final;" the program will run the requested number of initial chains, and then the requested number of final chains. It is often useful to make the initial chains shorter, giving a quick-and-dirty estimate which the final chains can refine, but the program does not dictate any particular relationship between initial and final.
The <initial> and <final> tags lay out parameters for the two chain types. Within them, <number> is the number of chains of that type, <samples> is the number of genealogies that will be sampled for parameter estimation, <interval> is the number of genealogies generated for each one that is sampled, and <discard> is the length of the burn-in period before any genealogies are sampled. For example, if <samples> is 35, <interval> is 10, and <discard> is 200, the program will first produce and discard 200 genealogies, and then produce 350 more, sampling every 10th one for a total of 35 sampled genealogies.
If you wish, for some reason, to run only one chain, you can set the samples of the other chain type to zero. For likelihood-based parameter estimation we believe that one chain is always too few. One chain is a reasonable choice for Bayesian analysis, however. (If only one chain is run, burn-in should probably be quite long.)
Note that this convention for indicating how many genealogies to sample is the same as the one used in MIGRATE, but different from the one used in COALESCE, FLUCTUATE, and RECOMBINE. The latter three programs take the total number of genealogies to be produced, not the number to be sampled. Be careful of this point when comparing runs of the different programs.
The <replicates> tag gives the number of replications for each chain (one or more). If more than one replicate is requested, a joint parameter estimate over all replicates will be produced.
The <heating> tag controls the heating strategy. It must contain two tags. <temperatures> gives a list of temperatures for the various searches (the number of entries in this tag determines the number of searches that will be run). The lowest temperature should always be 1.0. (We know of no use for a chain colder than that.) The <swap-interval> tag gives the number of chain steps that will pass between each attempt to swap trees among the different temperatures. 1 is a reasonable default; higher values may run a little faster, but are less effective.
A third, optional tag <adaptive> can be used to turn on or off "adaptive" heating, in which the program will adjust the temperatures in use so as to try to optimize acceptance rates. Adaptive heating is off by default. If you enable it (<adaptive> true </adaptive>), keep a close eye on your searches to make sure that they are performing well. In theory adaptive heating should be superior to fixed-temperature heating, but we have little experience with it so far.
The <strategy> tag indicates the rearrangement strategy. It contains a relative frequency (a number between 0 and 1) for use of this particular strategy; the frequencies of all strategies should add up to 1.
In version 2.0, four strategies are possible.
The <resimulating> strategy rearranges the genealogy. Every run must have a <resimulating> strategy, since only this strategy allows free movement through the genealogy search space. Its frequency should always be fairly high; we recommend at least 80% in a likelihood run, and at least 40% in a Bayesian run.
The <haplotyping> strategy reconsiders haplotype assignments. This is only useful if you have some phase-unknown sites in your data, and will be silently disabled, even if you specify it, on data with no phase-unknown sites. A reasonable frequency for it, if it is needed, might be around 20%.
The <bayesian> strategy allows a Bayesian analysis to search the space of population parameters. It is only useful in a Bayesian run (including a Bayesian arranger will turn the run into a Bayesian run). We have found setting its frequency equal to the resimulating arranger's frequency to be satisfactory. Please note that a Bayesian run with a 50/50 resimulating/Bayesian arranger strategy will have to take twice as many steps to consider as many trees as a pure likelihood run, and probably should. Bayesian rearrangements are relatively quick, so this will not slow the program inordinately.
The <trait-arranger> strategy is useful only when trying to map a trait to a location on the chromosome using the "jumping" algorithm. It allows the search to reconsider the location of the trait, and should probably have a high frequency such as 30-50%. It will be silently turned off if not needed.
Here is a sample search-strategy block.
<chains>
<replicates>1</replicates>
<heating>
<temperatures>1.0</temperatures>
<swap-intervals>1 1</swap-intervals>
</heating>
<strategy>
<resimulating> 1.0 </resimulating>
</strategy>
<initial>
<number>5</number>
<samples>50</samples>
<discard>100</discard>
<interval>20</interval>
</initial>
<final>
<number>1</number>
<samples>100</samples>
<discard>100</discard>
<interval>20</interval>
</final>
</chains>
The <format> tag encloses various options for formatting and detail level of results, as well as the random number seed.
The <convert-output-to-eliminate-zero> tag indicates whether you think of your data in terms of it having a 'site 0' or not. The traditional biologist assumption is that 'site 1' is right next to 'site -1', with nothing inbetween. As you might imagine, such a scheme would cause no end of trouble computationally, so a LAMARC input file scoots all such negative numbers up by one such that 'site -1' becomes 0, 'site -2' become -1, and so on. (The converter lam_conv does not make this assumption, but produces LAMARC input files that do.) If you are using the mapping option, you will see numbers for sites in both the menu and in the output file; this option can be toggled ('true' or 'false') by hand here so that the data output will match your expectations. See "Does LAMARC use 'site 0'?" in the FAQ.
The <verbosity> tag indicates how lengthy the output report should be; options are concise, normal, and verbose. When in doubt, try verbose--you can always discard unneeded parts later.
Similarly, the <progress-reports> tag indicates what kind of feedback should be given while the program is running. Options are none, concise, normal and verbose. For exploratory runs, we strongly recommend verbose.
The <seed> tag initializes the random number generator. The number provided should be an integer of the form 4N+1, such as 101 or 105. If the program is run twice with the same parameters, data, and seed, it will produce the same results.
If no <seed> tag is given and no seed is set from the menu, the program will try to use the system clock to seed the random number generator. If your system is peculiar enough to have no system clock, you'll want to prevent this by always specifying a seed.<results-file> gives the name of the file which will receive the output report when the program finishes (destroying any previous contents).
<use-in-summary-file> can be set to either "true" (if data reading from a summary file is desired) or "false" (if not). <in-summary-file> gives the name of the file from which to read in data from a previous run of LAMARC.
<use-out-summary-file> can be set to either "true" (if data writing to a summary file is desired) or "false" (if not). <out-summary-file> gives the name of the file which will contain the output summary file.
<use-curvefiles> can be set to "true" in a Bayesian run for which curvefiles should be saved, or "false" otherwise. <curvefile-prefix> gives a prefix which will form the first part of the name of all such curvefiles (the remainder of the name is formed from the name of the parameter whose curve is being recorded).
<use-reclocfile> can be set to "true" if you want to dump out the locations of all recombination events in sampled trees in the last final chain. This is "false" by default since the files can be quite large. <reclocfile-prefix> gives a prefix which will form the first part of the name of all recombination location files.
<use-tracefile> can be set to "true" if Tracer-readable summaries of the run are desired, or "false" otherwise. <tracefile-prefix> gives a prefix which will form the first part of the name of all Tracer files.
<use-newicktreefile> can be set to "true" if trees from this run should be written out in Newick format, or "false" otherwise. <newicktreefile-prefix> gives a prefix which will form the first part of the name of all Newick-format tree files.
<out-xml-file> gives the name of the file which will receive a copy of the input file, as modified by options selected in the menu.
There is no entry for the name of the input data file, since that's the file this XML is in.
Here is a sample <format> block:
<format> <verbosity>verbose</verbosity> <progress-reports>normal</progress-reports> <seed>1005</seed> <use-in-summary-file>false</use-in-summary-file> <in-summary-file>insumfile.xml</in-summary-file> <use-out-summary-file>false</use-out-summary-file> <out-summary-file>outsumfile.xml</out-summary-file> <out-xml-file>menusettings_infile.xml</out-xml-file> </format>
Trait mapping is enabled by presence of a <traits> block nested within a <region>. Within it there should be a <trait> block for each trait being mapped. Although you are unlikely to need it, if you have multiple traits that map to the same genomic region, you may map them simultaneously by giving each a 'trait' tag. The trait is assumed to be located somewhere in this region; LAMARC will not consider the possibility that it is elsewhere.
The trait is given a name with a <name> block, and the type of analysis (floating or jumping) is given with an <analysis> block whose legal values are "float" or "jump".
If information is available about possible locations for the trait, they can be indicated in a <possible-locations> block containing one or more <range> blocks. Each <range> should contain a <start> and an <end> block giving the starting and ending sites of the range, numbering the first site of the region as 1. At least one site must be a legal location for the trait (and really, it is silly to attempt mapping unless more than one site is legal).
The trait block may also contain an appropriate data model for trait allele mutations, which must (at present) be a K-Allele model. It is inappropriate to have multiple categories for a trait allele (since the categories involve differences in mutation rate among markers, and there will always be exactly one marker for trait alleles), and normalization can usually be ignored, but you may want to input a relative mu rate. For example, if you believe your alleles are the result of a dysfunctional gene, and you are working with SNP data, you might estimate that there are approximately 500 sites in an average gene that lead to that gene's disruption, and observed SNPs approximately one every 100 sites, so you would tell LAMARC that the relative mutation rate of your trait allele is 5 times more than that for your SNP data. (If you had DNA data instead, you would put in '500' as the relative mutation rate.)
As an example, here is XML input for the trait 'funny-nose', which we know might be mappable to site 1 or somewhere within sites 35-68, and which has a relative mutation rate of 5:
<data>
<region name="Region1">
<traits>
<trait>
<name> funny-nose </name>
<analysis> float </analysis>
<possible-locations>
<range>
<start> 1 </start>
<end> 1 </end>
</range>
<range>
<start> 35 </start>
<end> 68 </end>
</range>
</possible-locations>
<model name="KAllele">
<normalize>false</normalize>
<categories>
<num-categories>1</num-categories>
<rates> 1</rates>
<probabilities> 1</probabilities>
<autocorrelation>1</autocorrelation>
</categories>
<relative-murate>5</relative-murate>
</model>
</trait>
</traits>
Additionally, within the <individual> tag of each individual from whom data was sampled, you will need to include a <genotype-resolutions> block specifying the possible genotypes which could correspond to that individual's phenotype. The first block within the <genotype-resolutions> block should be a <trait-name> matching the name given earlier for the trait. This is followed by a collection of <haplotypes> blocks which give each haplotype possible for that individual, and the penetrance. Be careful with these penetrances! The penetrance of a haplotype is the chance that, if the individual has this haplotype, they will show the phenotype they did in fact show; it is not the chance that the individual has this haplotype. Only haplotypes which the individual could have should be listed. The penetrance, which should be a number greater than 0 and less than or equal to 1, goes in a <penetrance> tag and the two alleles, in order, which make up the haplotype go in a <alleles> tag.
Here is an example from an individual whose nose is bent. Based on our previous analysis of this trait, we know that 100% of individuals who are homozygous for the 'broken' allele have a bent nose, that 50% of individuals who are heterozygous and have one 'broken' allele and one 'normal' allele have a bent nose, and that no individual who is homozygous for the 'normal' allele has a bent nose:
<individual name="Mary">
<genotype-resolutions>
<trait-name> funny-nose </trait-name>
<haplotypes>
<penetrance> 1 </penetrance>
<alleles> broken broken </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.5 </penetrance>
<alleles> normal broken </alleles>
</haplotypes>
<haplotypes>
<penetrance> 0.5 </penetrance>
<alleles> broken normal </alleles>
</haplotypes>
</genotype-resolutions>
Note that we must include both 'normal broken' as well as 'broken normal' for the heterozygote.
Every individual will require such a block giving all possible haplotype resolutions of the trait data.