The Data¶
To use phycoeval, you will need biallelic genetic characters aligned across the species for which you wish to infer a phylogeny. Phycoeval assumes these characters are unlinked, but, depending on your data, you might be better off violating this assumption of the model (more about this below). Biallelic means that each character has two possible states. But, most genomic data are composed of nucleotides, which have four states. This can be accommodated by coding the nucleotides as biallelic.
Coding nucleotide characters¶
Phycoeval will read your nucleotide characters, and recode them as biallelic. If phycoeval finds a polyallelic character with more than two states across your samples, it will report an error and stop running. However, you can tell phycoeval to recode these polyallelic characters as either being the first state (0) or a different state (1). Your other option is to exclude these sites.
Note
If you are analyzing nucleotide data as biallelic characters, we strongly recommend that you do not try and estimate the frequencies of the two states. Instead, fix the frequencies to be equal. This is because there are many ways to recode the 4 nucleotide states to two states. Thus, if you try to estimate frequencies of the two states, your results can be sensitive to how you decided to code your nucleotides as binary.
Unlinked¶
Unlinked means that each character is assumed to have evolved along a gene tree that is independent from the others (conditional on the population history). In other words, it assumes that the characters are far enough apart from one another in the genome that they segregate independently during meiosis.
Many genomic data sets consist of many loci, each of which is a stretch of linked nucleotides. Examples include “RADseq” or sequence-capture data. What should we do with such data?
What to do with linked characters?¶
What if you have loci comprising sequences of linked nucleotides? Based on simulations [14][18][16][20], we recommend that you analyze all of your characters (including the constant ones) and violate the assumption of unlinked characters. Ecoevolity and phycoeval perform better when you use all of the sites (including linked and constant ones) compared to reducing the data to only one variable character per locus. So, your best bet is to include all of the sites from loci.
Formatting your data¶
Currently, phycoeval accepts two input formats for genetic characters: Nexus and YAML.
Nexus format¶
Standard haploid data¶
You can represent your data in a “standard” 0/1 format. Here’s an example of three species from which we’ve sampled 4 genomes each (2 diploid individuals):
#NEXUS
BEGIN TAXA;
DIMENSIONS NTAX=12;
TAXLABELS
lizard-953-a-speciesA
lizard-953-b-speciesA
lizard-954-a-speciesA
lizard-954-b-speciesA
lizard-152-a-speciesB
lizard-152-b-speciesB
lizard-154-a-speciesB
lizard-154-b-speciesB
lizard-331-a-speciesC
lizard-331-b-speciesC
lizard-338-a-speciesC
lizard-338-b-speciesC
;
END;
BEGIN CHARACTERS;
DIMENSIONS NCHAR=273658;
FORMAT DATATYPE=STANDARD SYMBOLS="01" MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA 001010...
lizard-953-b-speciesA 001010...
lizard-954-a-speciesA 101010...
lizard-954-b-speciesA 001011...
lizard-152-a-speciesB 101110...
lizard-152-b-speciesB 101110...
lizard-154-a-speciesB 001011...
lizard-154-b-speciesB 101010...
lizard-331-a-speciesC 001010...
lizard-331-b-speciesC 001010...
lizard-338-a-speciesC 011010...
lizard-338-b-speciesC 001011...
;
END;
Note, we don’t need separate TAXA and CHARACTER blocks like above. Instead, we can specify a DATA block:
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=12 NCHAR=273658;
FORMAT DATATYPE=STANDARD SYMBOLS="01" MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA 001010...
lizard-953-b-speciesA 001010...
lizard-954-a-speciesA 101010...
lizard-954-b-speciesA 001011...
lizard-152-a-speciesB 101110...
lizard-152-b-speciesB 101110...
lizard-154-a-speciesB 001011...
lizard-154-b-speciesB 101010...
lizard-331-a-speciesC 001010...
lizard-331-b-speciesC 001010...
lizard-338-a-speciesC 011010...
lizard-338-b-speciesC 001011...
;
END;
Both examples above would be equivalent for phycoeval, but the Nexus Class Library used by phycoeval will report a message about an implicit TAXA block if you use the latter format. Either way, in your phycoeval config file, you need to tell phycoeval that the states, or genotypes, are haploid by declaring:
genotypes_are_diploid: false
Standard diploid data¶
Above, each cell in our matrix represented which state was present for the character in a particular haploid genome. We can also represent the same data where each cell represents the genotype of a diploid individual:
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=6 NCHAR=273658;
FORMAT DATATYPE=STANDARD SYMBOLS="012" MISSING=? GAP=-;
MATRIX
lizard-953-speciesA 002020...
lizard-954-speciesA 102021...
lizard-152-speciesB 202220...
lizard-154-speciesB 102021...
lizard-331-speciesC 002020...
lizard-338-speciesC 012021...
;
END;
Now, “0” represents that the individual has two copies with the 0 state, “2” represents two copies of the 1 state, and “1” represents a heterozygote. Again, in your phycoeval config file, you need to tell phycoeval that the states, or genotypes, are diploid by declaring:
genotypes_are_diploid: true
Nucleotide data¶
If you have nucleotide data, the easiest thing is provide the nucleotide characters to phycoeval as is, and let it recode them as biallelic. Here’s an example where we are providing nucleotides as haploid (each cell is a haploid genotype):
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=12 NCHAR=273658;
FORMAT DATATYPE=DNA MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA ACGTAG...
lizard-953-b-speciesA ACGTAG...
lizard-954-a-speciesA GCGTAG...
lizard-954-b-speciesA ACGTAA...
lizard-152-a-speciesB GCGCAG...
lizard-152-b-speciesB GCGCAG...
lizard-154-a-speciesB ACGTAA...
lizard-154-b-speciesB GCGTAG...
lizard-331-a-speciesC ACGTAG...
lizard-331-b-speciesC ACGTAG...
lizard-338-a-speciesC ATGTAG...
lizard-338-b-speciesC ACGTAA...
;
END;
This is sometimes referred to as “phased” data. Again, if we are providing a matrix where each cell represents a haploid genotype, we need to tell phycoeval this is so via the config file:
genotypes_are_diploid: false
We can also represent the same data as “unphased”, where each cell represents a diploid genotype:
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=6 NCHAR=273658;
FORMAT DATATYPE=DNA MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA ACGTAG...
lizard-954-a-speciesA RCGTAR...
lizard-152-a-speciesB GCGCAG...
lizard-154-a-speciesB RCGTAR...
lizard-331-a-speciesC ACGTAG...
lizard-338-a-speciesC AYGTAR...
;
END;
We need to indicate this in the config file:: accordingly:
genotypes_are_diploid: true
Population labels¶
In our nexus character matrix, we need to indicate which species (or population) each row corresponds to. We can do this with either using a prefix or suffix in the row (or taxon) labels. For example, in this nexus data file:
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=6 NCHAR=273658;
FORMAT DATATYPE=DNA MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA ACGTAG...
lizard-954-a-speciesA RCGTAR...
lizard-152-a-speciesB GCGCAG...
lizard-154-a-speciesB RCGTAR...
lizard-331-a-speciesC ACGTAG...
lizard-338-a-speciesC AYGTAR...
;
END;
we are using the suffixes to indicate that the first two samples came
from a species we are calling speciesA,
the next to samples came from a species called speciesB,
and
the last two samples came from a species we are calling
speciesC.
In our phycoeval config file
we have to indicate this with:
population_name_delimiter: "-"
population_name_is_prefix: false
This tells phycoeval to look for the last bit of each row label that is separated by a “-” to figure out the population label.
Note
If you like to use underscores as a population label delimiter, just watch out for a gotcha related to how the nexus format treats underscores
YAML format¶
The YAML-formatted data format is a lot more efficient (much smaller file sizes). Instead of a full alignment, it only contains a list of allele count patterns, followed by a list of the weight of each pattern (i.e., how many times the allele pattern occurs in the alignment). For a very small example let’s convert the following nexus alignment and convert it to the YAML format that ecoevolity accepts as input:
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=12 NCHAR=6;
FORMAT DATATYPE=DNA MISSING=? GAP=-;
MATRIX
lizard-953-a-speciesA ACGTAG
lizard-953-b-speciesA ACGTAG
lizard-954-a-speciesA GCGTAG
lizard-954-b-speciesA ACGTAA
lizard-152-a-speciesB GCGCAG
lizard-152-b-speciesB GCGCAG
lizard-154-a-speciesB ACGTAA
lizard-154-b-speciesB GCGTAG
lizard-331-a-speciesC ACGTAG
lizard-331-b-speciesC ACGTAG
lizard-338-a-speciesC ATGTAG
lizard-338-b-speciesC ACGTAA
;
END;
To convert these data to YAML format, we will assume the first nucleotide (from the top) in each column is state “0”, and the second nucleotide (if any) is state “1”. Doing so gives us the following allele-count patterns in YAML format:
---
markers_are_dominant: false
population_labels:
- speciesA
- speciesB
- speciesC
allele_count_patterns:
- [[1,4], [3,4], [0,4]]
- [[0,4], [0,4], [1,4]]
- [[0,4], [0,4], [0,4]]
- [[0,4], [2,4], [0,4]]
- [[1,4], [1,4], [1,4]]
pattern_weights:
- 1
- 1
- 2
- 1
- 1