Repository navigation
Activity 2
In this activity, we will
- Simulate a DNA sequence alignments for a fixed tree using Seq-Gen under the Jukes-Cantor, K2P, and HKY models. showing how to generate a multiple sequence alignment
- Generate multilocus datasets by combining gene-tree simulation with sequence simulation.
- Use SeqKit and ape to inspect simulated alignments and explore some implications and extensions of the models presented here.
For this activity, we will use Seq-Gen to run the Jukes-Cantor process on the primate gene tree from Exercise 2. First, save the Newick tree as a file using the following shell command:
echo "(((H:.052,C:.052):.107,G:.167):.052,O:.233);" > primate-gene.tree
We now use primate-gene.tree as the input. To simulate the Jukes-Cantor
model, we will run seq-gen with the options -mHKY (actually, the -mHKY
option specifies a more general model HKY85, for which the Jukes-Cantor model
is a special case; by not specifying the additional parameters of that model,
the defaults give the Jukes-Cantor model). The following code generates four
aligned sequences of length 60.
seq-gen -mHKY -l60 primate-gene.tree > primate-gene.fasta
cat primate-gene.fasta
Example output file
>H
GGCCTCTGGCGTATGGAGGGGCATGAGGCTATTCTGGTACGGCGCGTGCTTGTGGGAAAA
>C
GTCCTCTGGGGTAAGTAGGGGCATGAGGCTATTCTGGTACGGCGCGTACTTGTAGGAAAA
>G
GTCTACTGGCGTAAGTAGGGGCATGAGACTATGCTGGTGCTAGGCACGCTTGTACGAAAG
>O
GCCTACTCGCGTAAGTAGCAGTATGAGGCTATGCTTGTCCCAGGCATGAATGTACCAAAA
The file primate-gene.fasta is an example of a FASTA file, which is a format
for representing a multiple sequence alignment. Each of the fours sequences is
preceded by a line (e.g., with >H or >C), which identifies the sequence below
it. The sequences are represented by strings of A,C,G, and Ts. In this case,
the strings are of length 60, but for whole genomes, these may be many
millions of letters long, and are often broken up into lines of between 60-80
characters each. For a full specification of the FASTA file format, see
here.
-
Generate a multiple sequence alignment with 500 base pairs for this gene tree.
Code
seq-gen -mHKY -l500 primate-gene.tree > primate-gene.fastaInspecting the output file with
cat primate-gene.fastagives>H CTCTACCCTACGGCTGGCGGAGTGAGCCATCACTGGGTGAGTCATTTCGATGCAGGCTGACCCGACCCCTAT TCATTCGCCCTATTCTTCCATCCACTGCCGACCATTCCTTTGTTCAGTCTGGATTAACGAAGAGCATCGGGC GGCCACAAGTTGAACCGCTCCTATGCCTCTATATATAAAGAGGTATTTTCACAGTCTTGATCTAGTACGCTC CACAGCCCCTCCCAAGAAGAACTACGACTTGGACGAACCAGCCTACCTGCTTCCCGCTGGGCGGATGACACC CGGGGGAGTACGCTTCGTCGGTCTTCTGGTACGTCTGTACCTGTCAGAATTAAGTATACGGCTCTGCTGTTC CACATGGCACAGATTGTGGGAGTTTGGCCCTATTGGTGTCTTGCCAGATCCCTCTAGATACTATACTGGCGT CTCACTATCACGAAAGGCTCCAACTACGGTTCCGAACTCGCTCCCTAGTGAGAATTAACGCGACGTTC >C CTCTATCCTACGGCTTGCGGAGGGAGCCCTCACTGGGTGAGTCATTTCGATGCAGGCTGACCCGACCCCTAT TCATTCGCCCTTTTCTTCCATCCAGTGCCGACCATTTCTTTGTTCTGTCTGGATTAACGAAGAGCATCGGGC GGCCACAAGTTGAGCCGCTCCTATGCCACTATATATAAATAGGTAATTTCACAGTCTTCAACTAGGACGCCC GACATGCCTTCCCAAGAAGAACTAGGACTTGGACTTACCAGCCTACCCGCTTCCCGCAGGGCGGATGACACC CGGGGGAGTACGCTTCGTCGGTCTTCTGGTACGTCTGCACCTGTCATTATTAAGTATACGGCTATGCTGTTA CACATGGCACAGATTGTGGGAGTTTATCCCTATTGGTGTCTTGCCAGACCCCTCAAGAGACTTTACTGGCGT CTCACTACCACGAAAGGCTCCAACTACGGTTCCAAACTCGCTCCCTAGTGAGAATTAACGCGTCGTTG >G ATCTACCGTTCGACTGGCGGAGAGAGCCATCACTGGGGGTGCCATCTCGATGTAGGCTGACCTGACCCCTAA TCATTCATCCGCTACTTCCATCGACTGTCGACCGTTCCTTTTTGCAGTCAGGATCAACGAAGATCATCGGGC GTCCACAAGGGGTACCGCTCGTATGCCAGTGAATTTAAAGAGGTATTGTCATAATCTTTATCCGGGACTCTC CACACCCCCTTCCACGCCGAAAGAAGACTTGTACGTACAAACATAACAGCCTCCCGCAGCGGAAATGTCACT GGGGTTAGTACGATTCGTCGACTCTCTAGTACAACTGGACCTATCATTATTAAGTTAACTGCGATGCTGTTA CACCTGGCACAGATCGTGGGAATTTATCCCTATTGGTGTCCTGACAGAGCACGCTTGACACGTTACTGTCCT CTCTCTATCACTACCGGCTACAACTCCTGTGCCTAACTCCAGCCCAACTTAGAATTAAGGCAGCGGTT >O ATCTACCCTACGACGGGGGCAGTGAGAGATCACTTAGTGACTCATCTCTAGGTAGGCTGACGTGACTCCTAT TCATCCGGACTTTACTAGCGTCCATGGACGACCATTCCTCTTTTAGGGGAGGACCCACGAAGAGCATCGAAT GCCGACAATTTTAATCTCTCCGAGGCCAGTGTACTTAACGAGTTATTATCACAATCACCATCTGGAGCGCTC CACACCCCCTCTCACGAAGAAATAAGACTTTTCACGACCATCCCACCTGAGTCCCACAGCGAGGATGGCACG CAGGTTAGAACTCTTCCTCGTTAGTCTGTTACGGATGCCTATGTGATAGGACAGTCGACTACTTTGCTGTTC CACCTAGAACCGCTTGGCGGAAATTGACCATATTGTTCACTGGACAGAGCATTCTAGACACTTACCTGGCAT GTCACTATCACTAACGGCTACAACTCCGGTGCCAAATCTCCGGGCTATTGAGAGTTAACGCAGCGTGA -
We can scale all the branches of the gene tree my adding option
-s#, where#is a number representing a scaling factor by which all branch lengths are multiplied. (For example,-s2would double the length of each branch before generating sequences.) This is useful for simulating variable mutation rates. It is also useful for unit conversion (seq-genassumes branch lengths measured in expected number of mutations per site, but gene trees may have branches measured in coalescent units). What patterns do you observe in your sequences when the scaling factor is very close to zero? When it is very large?Answer
Small scaling:
seq-gen -mHKY -s.01 -l60 primate-gene.treereturns sequences that are highly correlated (or even identical):>H TCGCCACCCTATCCTTTCCCCACCACAGTATGGTAAGGTCAAATCTGCATACGGACTGAC >C TCGCCACCCTATCCTTTCCCCACCACAGTATGGTAAGGTCAAATCTGCATACGGACTGAC >G TCGCCACCCTATCCTTTCCCCACCACAGTATGGTAAGGTCAAATCTGCATACGGACTGAC >O TCGCCACCCTATCCTTTCCCCACCACAGTATGGTAAGGTCAAATCTGCATACGGACTGACLarge scaling:
seq-gen -mHKY -t1 -s10 -l60 primate-gene.treereturns sequences that are nearly independent.>H CATTTTTACGAGACTTAGTTGCGCGACTGGGCAACCTTATCGACTTATCTGGGGCTGCAG >C CATAGCGCCGAGTCGCATGCGCCAGACAGGTTGTGGGAGACGAGTGTTCGTAGGAGCGAG >G CTCGATGTTGGGAACGCGCCATCAGCTACGTGGAAAACAGCATAATTGGGACTTGGCAGC >O CAAGTTCAGATTGATGTTGGGCGAGGTAGGCTTATTATTGCGGCCCTTGAGTTTAATCAA -
Adding the options
-opor-onwill change the format of the output to PHYLIP or NEXUS respectively. Try these options and see how the output changes.NEXUS Format
To generate NEXUS file, run
seq-gen -mHKY -l60 -on primate-gene.tree > primate-gene.nex cat primate-gene.nexgives
#NEXUS [ Generated by seq-gen Version 1.3.5 Simulations of -1 taxa, 60 nucleotides for 1 tree(s) with 1 dataset(s) per tree Branch lengths assumed to be number of substitutions per site Rate homogeneity of sites. Model = HKY: Hasegawa, Kishino & Yano (1985) Rate of transitions and transversions equal: transition/transversion ratio = 0.5 (K=1) with nucleotide frequencies equal. Output format: NEXUS Writing all datasets to standard output ] Begin DATA; Dimensions NTAX=4 NCHAR=60; Format MISSING=? GAP=- DATATYPE=DNA; Matrix H AGCGCAGGAACGTCACACACAAATCGCTAAGCCGCGCCCGCTGGCTTTTTTATGTAGGTG C AGCGCAGGAACGTCACACACAAATCGCTAACCCGCGACCACTGGCTTTTTAATGTAGGTG G GGGGCCGGCACGTCACGCACAAATCGCAGAGCCGCGCCACCTGGATTTTTTCCGTAAGTG O TGGGCCGGCATGTCATGCCCCAATCCCTCAGGTGCGCCCCCTGGATTTGTTCCGTAGGTG ; END;PHYLIP Format
To generate PHYLIP file format, run
seq-gen -mHKY -l60 -op primate-gene.tree > primate-gene.phy cat primate-gene.phygives
4 60 H TGGTATGCACTTTAACGGCAAGACATTCTACGAGGTTCGTCTTACAGTAGTGGGTGAATG C TGGTATGCACTTTAACGGAAAGACATTCTACGAGGTACGTCTGACAGTAGTGGGTGCATG G TGGTAGGCACGTTAACGGATAGATATTCTGCGAGGAACGTCTGACAGGAGTGGGGGAACC O TGGTATTCACCTTCTCAGAAAGACTGTCGGCGGGGGATGGCTGACAGCAGGTGCGGATTG
In this activity, we will generate a synthetic genome consisting of
ms 4 2 -T -I 4 1 1 1 1 -ej 1 4 3 -ej 2 3 2 -ej 4 2 1 | grep "^(" > two-gene-trees.trees
We can then inspect the files using
cat two-gene-trees.trees
Sample output
(1:4.980,(2:2.286,(3:1.909,4:1.909):0.377):2.693);
(1:5.261,(2:3.848,(3:1.008,4:1.008):2.840):1.413);
The following code will output two FASTA files, one for each of the two genes.
seq-gen -mHKY -l144 -of -y "primate-gene" < two-gene-trees.trees
The option -l144 tells seq-gen to generate genes of length 144 bp. The option -y determines the names of the output files.
Observe that the output tells us that two gene sequences were generated, one for each tree:
Simulations of 4 taxa, 144 nucleotides
for 2 tree(s) with 1 dataset(s) per tree
The two output files are named primate-gene_1.fasta and
primate-gene_2.fasta.
Sample output
The content of the two files will look something like this:
>1
CTTTACGTGACGGCGTCCCTACTCGCCCGGACCATGCTCTTACTGTTCAGAAGATTCGAATTGATCAATGTT
AGCAGCGCCTTTAATGACGTAAAGGCACGGCTTTCACAACTTACAAGCATTGGAAACCCTCTTCAGTCTCAG
>2
AAAGAAAGTGGCGTGGACAAGGTTAAAGTCAGTTCGCTAAAAGTGGCGCCGGGCTCTCAGATACGTATAGAC
TCTTCCCACCCGTAGGGATGCATTATCGCCCTGTATTCGATCTAAAACAGCTTTGTCTTCACTAATATATTG
>3
CCGTCTCGCCTCATTCCCGGTCTCATAAGTGACGACGCTCGCCGGGGGCAAATCCTATGGAGGTCATCAGTC
AGATCCGCCTGTGGTGGGCATGATGCGAGGAAACCCCGACGTTGTGCCGGAAACAAGAACCATCGCTCCACA
>4
CGCGCACGTCCCCTTTCATAGTTAATGAACCCATGGGAGAGAGCGGCCCCACGAGTATACTAAACGGTGCTT
GTTCGCATGATGACCCGATACGACGTAACGTCTCCAGAGCTTCCGTTTTGGCCTTACCGTAGGCCCGCTATC
and
>1
CAGTCCTACTAACTGTCGTGCTTGTTTACCAAGGTGAATATAGTTCCCACTATTCTGAATCTTATAGGTACG
TCTTGTAACGAGATGATCGGTTATGTCAGGTCCGCCCACTGTGATGGTTTACCTTGAACTCAATGCCGCAGG
>2
AGACACAAGCGTTTCTAGCTTCCTATGTCGACGATCGAATAACATAATCTTTGGCGCGGCTAGTATTTTAAC
GTTGGTGGACATACGGAATAATCATGCTCGGTAGACTTTCAATCTGAATGAACTGAATACGAGCCGGATGCC
>3
AGAGTTTCGGTTTTGAAAGAGTCAGCAAGAAACGCATTAGTTTGAGTTTGGCTATAGCTACACGAATGAGGG
GAATTCTATTTGCTTGGATCTTGCCACAACCATAAACTTGAGCAGGATTCGCATAAAGCGGCGGGTCGCTCT
>4
GGTCCGGAGATTGACATCTGTTTCAAGTGGCGAGGCGATCACCCGCACACACTGCTATAACTCAATAACGTC
GTGGGGCGACAGGGGGTGTCGAAAGAAGAACATCGGTTACACTTAATGACGTCGTGAATAGGAAGTGCTTGG
Next, we concatenate the two FASTA files primate-gene_1.fasta and
primate-gene_1.fasta together into a single FASTA file using the software
SeqKit, which has a function concat designed for this purpose.
seqkit concat primate-gene_1.fasta primate-gene_2.fasta --out-file primate-gene.fasta
Alternatively, we could run the following, which uses shell expansion:
seqkit concat primate-gene_*.fasta --out-file primate-combined.fasta
This writes the output to the file primate-combined.fasta. Note that the
order of the taxa has been rearranged, but this is harmless.
Next, we will generate a larger alignment consisting of
(1:1,(2:.5,(3:.25,4:.25):.25):.5);
It will be helpful to inspect this tree now for later comparison. We can plot it in R using
speciesTree = read.tree(text ="(1:1,(2:.5,(3:.25,4:.25):.25):.5);")
plot(speciesTree, cex=2)
edgelabels(speciesTree$edge.length, cex=2)
See species tree
The following shell code will generate a our larger multisequence alignment.
First, it will use the multispecies coalescent to generate primate-combined-msa.fasta. The
output will all be stored in the directory primate-genes/.
mkdir -p primate-genes
ms 4 1000 -T -I 4 1 1 1 1 -ej 0.25 4 3 -ej 0.5 3 2 -ej 1.0 2 1 | grep "^(" > thousand-gene-trees.trees
seq-gen -mHKY -l500 -of -y "primate-genes/gene" < thousand-gene-trees.trees
cd primate-genes
seqkit concat gene_*.fasta --out-file primate-combined-msa.fasta
we can view basic summary statistics about the combined dataset with
seqkit stats primate-combined-msa.fasta
Sample output
file format type num_seqs sum_len min_len avg_len max_len
primate-combined-msa.fasta FASTA DNA 4 2,000,000 500,000 500,000 500,000
We can also view the sequence IDs using the command
seqkit seq --name primate-combined-msa.fasta
or the command
cat primate-combined-msa.fasta | grep ">"
This is useful for very large FASTA files, where visual inspection may not be feasible.
Sample output
3 2 1 4
In R, we can use the package ape to read additional data from the alignment
primate-combined-msa.fasta.
First, we load the alignment using the command read.dna(), and inspecting
the output shows the alignment length, IDs, and nucleotide base composition of
the dataset:
# library(ape)
# setwd("PATH/primate-genes")
M <- read.dna("primate-combined-msa.fasta", format = "fasta")
M
Sample output
4 DNA sequences in binary format stored in a matrix.
All sequences of same length: 500000
Labels:
2
1
4
3
Base composition:
a c g t
0.25 0.25 0.25 0.25
(Total: 2 Mb)
More precise base composition in the dataset can be obtained using
base.freq(M), and running this, we see that the frequencies of
Sample output
a c g t
0.2497295 0.2503420 0.2500945 0.2498340
-
Use the R function
dist.dna(from the ape library) to estimate pairwise distances between the four speices 1,2,3, and 4. Compare the estimated branch lengths with those of the species tree. (Here, evolutionary distance is measured in expected number of mutations per site).Show answer
The command
dist.dna(M, model="JC69")returns a distance matrix like this:
2 1 4 1 2.590000 4 1.615641 2.653829 3 1.628510 2.681735 1.137024 </details>The estimated pairwise distances are longer than the distances on the species tree. For example, one the species tree, the distance between species 1 and 2 is 2, but our estimated distance is 2.45.
-
Generate an MSA with 500,000 sites by running the Jukes Cantor process directly on the species tree
(1:1,(2:.5,(3:.25,4:.25):.25):.5);(i.e., without generating gene trees). Estimate distances again, and compare your results with the answer to the previous problem.Show answer
Run the following shell code from
primate-genes/echo "(1:1,(2:.5,(3:.25,4:.25):.25):.5);" > speciesTree.tree seq-gen -mHKY -l500000 -of -y "no-ILS-msa" < speciesTree.treeThis will return a single fasta file containing four aligned sequences of 500,000bp each, generated using the Jukes-Cantor model. (To check, run
seqkit stats no-ILS-msa_1.fasta).Switching to R (in the same directory
primate-genes/), runM1 <- read.dna("no-ILS-msa_1.fasta", format = "fasta") dist.dna(M1, model = "JC69")This returns
1 2 3 2 2.0002964 3 2.0134847 1.0006371 4 2.0027769 1.0014880 0.4992852The branch length estimates now closely match the branch lengths of the original species tree.
-
What explains the discrepancy between the distance estimates in the previous two problems?
Show answer
The overestimation of distances occurs when gene trees are incorporated into the model. This is because divergence times on the gene tree are always at least as old as the divergence times of the species tree.
For example, the divergence time of species A and B is when the two populations split. But if you sample one gene from each of the descendant species, the two gene copies do not coalesce exactly at the species split. Instead, they enter the ancestral population and then have to wait additional time before coalescing.
-
In the above section, we assumed that all genes evolve at the same rate (i.e., that their mutation rates are equal). We can relax this assumption by using the scaling factor option
-s. Generate a multisequence alignment in which each gene is scaled by an independent gamma-distributed random variable with shape and rate parameters both equal to 2.rgamma(1, shape = 2, rate = 2) -
Another option for modeling rate heterogeneity is the rates-across-sites (RAS) model, in which the mutation rate of each site is multiplied by an independent gamma-distributed scaling factor (with mean 1) [Yang 1993]. This can be implemented in the
seq-genfunction using the optional flag-a#where#is a number representing the shape parameter$\alpha$ of the gamma distribution. This results in mutation rates that are gamma-distributed with variance$1/\alpha$ , so smaller$\alpha$ corresponds to higher heterogeneity).On the other hand, some sites are conserved (e.g., due to selection) making their effective mutation rate 0. These are called invariant sites. This can be implemented with the optional flag
-i#, where#is the desired proportion of invariant sites.Simulate a multisequence alignment using the rates-across-sites model, in which the proportion of invariant sites is .2.
-
So far we have assumed that all genes have the same sequence length,
$k$ . Let's relax this assumption. Generate an alignment with$m=1000$ genes whose lengths (in bp) are approximately normally distributed with mean 500 and standard deviation 50 (e.g., for each gene,k = round(rnorm(1, mean=500, sd=1),0)).Hints
You can concatenate the genes with SeqKit, or alternatively, with the python script AMAS.py, using a command like this:
AMAS.py concat -i gene_*.fasta -f fasta -d dna -p partitions.txt -t output-file.fasta(This has the advantage over SeqKit of producing a partition file, which tells you where each gene starts and ends, which SeqKit does not support).
The assumptions made by the Jukes-Cantor model are quite restrictive. In this section, we explore models which incorporate more realistic assumptions about DNA mutation.
Let
The Jukes-Cantor model assumes that all nucleotides mutate to each other at
equal rates. Letting
Equivalently, we can represent this by the rate matrix:
In this matrix,
The structure of this matrix therefore encodes the the symmetry assumption
that each state is equally likely to mutate to any other state. This symmetry
is not realistic, due to the molecular structure of purines (
In this activity, we will relax the Jukes-Cantor symmetry assumption by considering a more realistic model called the K80 Model (Kimura 1980). This is also called the K2P (Kimura 2-parameter) model because it has rate parameters for both transitions and transversion.
If we assume that transitions occur at rate
The transition/transversion ratio is
Note that when
Using Seq-Gen, we can generate data according to the K2P model by passing the
function seq-gen with the options -mHKY -t#, where # specifies the
transition/transversion ratio.
-
What is the state diagram for the K2P Model?
Show answer
-
Recall that when we generated alignments for the Jukes-Cantor model, the base frequencies were very close to
Pr(A) = Pr(C) = P(G) = P(T) = \frac{1}{4}.Is that still the case for the K2P model? Generate a large alignment using the K2P substitution model with any input tree and any value value for
-tother than$0.5$ , and compute the base frequencies.Show answer
Note: here we will use the R command
system()to shell commands directly from within R.In R, we run
system('echo "(1:1,(2:.5,(3:.25,4:.25):.25):.5);" > speciesTree.tree') system('seq-gen -mHKY -l500000 -t3 -of "speciesTree.tree" -y "K2P-t3-msa"')The choice of species tree was arbitrary. Here, we've chosen
$t=3$ .To get base frequencies, we can runbase.freq(read.dna("K2P-t3-msa_1.fasta", format = "fasta"))which outputs
a c g t 0.250440 0.250011 0.249601 0.249948This is the same result we obtained under the earlier Jukes-Cantor model: the base frequencies are very close to
$1/4$ (and indeed, they converge to$1/4$ as the dataset grows).
In the previous exercise, we saw that the equilibrium base distribution under
both the Jukes-Cantor and K2P models is uniform: the nucleotide bases
The assumption that all base frequencies are close to
-
Generate to generate sequences using the HKY85 dataset, by using the seq-gen function with the
-foption to specify equilibrium frequency distributions.Example code
In R, run
system('echo "(1:1,(2:.5,(3:.25,4:.25):.25):.5);" > speciesTree.tree') system('seq-gen -mHKY -f.15,.25,.40,.20 -l500000 -t3 -of "speciesTree.tree" -y "K2P-t3-msa"') base.freq(read.dna("K2P-t3-msa_1.fasta", format = "fasta")) -
The most general "time-reversible" model is the General Time Reversible (GTR) model, which can be implemented with the
seq-genfunction using the flag-mGTR. The GTR substitution supports up to 6 rate parameters with the flag-r x1,x2,x3,x4,x5,x6, where-
x1:
$\mathtt{A} \leftrightarrow \mathtt{C}$ (rate of mutation from$\mathtt{A}$ to$\mathtt{C}$ and vis-versa) -
x2:
$\mathtt{A} \leftrightarrow \mathtt{G}$ -
x3:
$\mathtt{A} \leftrightarrow \mathtt{T}$ -
x4:
$\mathtt{C} \leftrightarrow \mathtt{G}$ -
x5:
$\mathtt{C} \leftrightarrow \mathtt{T}$ -
x6:
$\mathtt{G} \leftrightarrow \mathtt{T}$
along with arbitrary base frequencies (e.g., set with the
-foption).When used along with the
-aflag (for variable rate sites), it is called the GTR+Γ model, and when also used with the-ioption for invariant sites, it is called the GTR+Γ+I model. This model is sufficiently general for many purposes, and indeed all of the other models we have considered are submodels of it.Simulate a multisequence alignment consisting of
$m=10000$ genes, each with$k=500$ base pairs under the GTR+Γ+I model, using arbitrary parameter values for the rates and base frequencies. -
In this section, we introduce the a standard model of hybridization known as the Network Multispecies Coalescent Model (NMSC).
When hybridization occurs between species, as in the case of the recent hybridization between polar bears and grizzly bears (which gave rise to a new hybrid species known as "pizzly bears"), the result is a non-treelike evolutionary history, like shown in the following diagram:
The evolutionary history in the above figure is called a phylogenetic network (or species network), since it is not a tree. The topology of this network can be represented in "extended Newick format" as
((Brown_Bear,(Grizzly_Bear,(Pizzly_Bear)#H1)),(#H1,Polar_Bear));
Here, the tag #H1 indicates the hybridization.
The dashed lines in this network are called reticulation edges. The parameter
-
What are the two trees displayed in this network, and what are their probabilities? Assume
$\gamma=0.1$ Show answer
The left tree has probability
$0.9$ . The right tree has probability$0.1$ .
-
Now suppose that the pizzly bear network has the following branch lengths (in coalescent units)
What are the branch lengths of the displayed trees, and their newick strings?
Show answer
Left tree:
((Brown:0.60,(Grizzly:0.35,Pizzly:0.35):0.25):0.45,Polar:1.05);Right tree:((Brown:0.60,Grizzly:0.60):0.45,(Pizzly:0.60,Polar:0.60):0.45);
Assuming the branch lengths are measured in coalescent units, we can generate 9000 gene trees from the left tree (the 'grizzly-pizzly' tree) using the following code:
ms 4 9000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 0.35 3 2 \
-ej 0.6 2 1 \
-ej 1.05 1 4 | grep "^(" > grizzly-pizzly.trees
In the above code, the leaf labels 1 = Brown bear, 2 = Grizzly bear, 3 = Pizzly bear, and 4 = Polar bear. Hence, the option -ej 0.35 3 2 says to merge pizzly bears into grizzly bears at time -n 1 2 -n 2 2 -n 3 2 -n 4 2 is a technicality, since we are measuring branch lengths in coalescent units, while ms measures time in twice coalescent units.
Similarly, we can generate 1000 gene trees from the right tree (the 'pizzly-polar' tree) with
ms 4 1000 -T \
-I 4 1 1 1 1 \
-n 1 2 -n 2 2 -n 3 2 -n 4 2 \
-ej 0.6 2 1 \
-ej 0.6 4 3 \
-ej 1.05 3 1 | grep "^(" > pizzly-polar.trees
We then combine and shuffle the gene trees using
cat grizzly-pizzly.trees pizzly-polar.trees | shuf > bear-gene.trees
Next we wish to generate sequences for these genes. At this point, it would be natural to try the following, but the following command probably won't work due to shell limitations:
DONT RUN THIS
mkdir -p bear-gene-sequences
seq-gen -mHKY -l1000 -of -y "bear-gene-sequences/sequence" < bear-gene.trees
Instead, we will create sequences in batches:
First create 100 batches of 1000 gene trees each
mkdir -p bear-gene-trees
split -a 3 --numeric-suffixes=1 -l 1000 bear-gene.trees bear-gene-trees/bear-batch-
Next, to simulate sequences for the first batch of 1000 genes, run
mkdir -p bear-gene-sequences/batch-001
seq-gen -mHKY -l1000 -of -y "bear-gene-sequences/batch-001/gene" < bear-gene-trees/bear-batch-001
and to generate sequences for all batches, we can run a shell loop:
for N in {001..100}; do
mkdir -p "bear-gene-sequences/batch-${N}"
seq-gen -mHKY -l1000 -of -y "bear-gene-sequences/batch-${N}/gene" < "bear-gene-trees/bear-batch-${N}"
done
The Julia package PhyloCoalSimulations.jl provides a possiblly simpler way to simulate the coalescent on a network.
using PhyloCoalSimulations, PhyloNetworks
net = readnewick("((Brown_Bear:0.60,(Grizzly_Bear:0.35,(Pizzly_Bear:0.25)#H1:0.10::0.90):0.25):0.45,(#H1:0.35::0.10,Polar_Bear:0.60):0.45);")
trees = simulatecoalescent(net,1000,1) # simulate 1000 gene trees, 1 sample per species
writemultinewick(trees, "bear-gene-trees.phy")
Explore what aspects of this activity can be done using PhyloCoalSimulations.jl.
Matthew William Hahn. Molecular population genetics. Oxford University Press (2018).
David A. Baum and Stacy D. Smith, Tree Thinking, W.H. Freeman and Company (2013).
Sebastián Duchêne, Simon Y.W. Ho and Edward C. Holmes, Declining transition/transversion ratios through time reveal limitations to the accuracy of nucleotide substitution models, BMC Evolutionary Biology (2015) DOI: 10.1186/s12862-015-0312-6
Motoo Kimura, A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences, Journal of Molecular Evolution (1980) DOI: 10.1007/BF01731581
Masami Hasegawa, Hirohisa Kishino, and Taka-aki Yano, Dating of the human-ape splitting by a molecular clock of mitochondrial DNA, Journal of Molecular Evolution (1985) DOI: 10.1007/BF02101694
Ziheng Yang, Maximum-likelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites, Molecular Biology and Evolution (1993) https://doi.org/10.1093/oxfordjournals.molbev.a040082