Phylogenetics: Compositional Heterogeneity Lab

From eebedia

Jump to: navigation, search
Adiantum.png EEB 5349: Phylogenetics



Simulated data and the true tree

The program p4, written by Peter Foster, specializes in simulating and analyzing data in which nucleotide composition varies across the tree. See Foster (2004) for an introduction to the potential problems created by nucleotide frequency heterogeneity. I used p4 to simulate data on the tree shown on the right. The black-colored lineages were characterized by an AT-biased nucleotide composition very different from the red-colored lineages, which were strongly GC-biased. My goal was to generate a data set that would be very susceptible to nucleotide compositional attraction under ordinary substitution models, and in that I succeeded (as you will see). Taxa C and H share many similarities due to the large number of G and C bases they have independently acquired from their AT-rich ancestors, and it will be very tempting for an ordinary nucleotide model such as GTR to place C and H together.

Note: do not use qlogin initially. If you are in a qlogin shell, you will not be allowed to submit a job using qsub (which can only be done from the head node, which is the machine that you are on when you first login to the cluster). So, if you qlogin now, you will need to get out of your qlogin shell using Ctrl-d before running nhPhyloBayes (which uses sneakily uses qsub to submit your job for you).


mkdir nhlab
cd nhlab
curl > sim500.nex


Running nhPhyloBayes






Ordinary models yield the compositional attraction tree



begin paup;
  log file=pauplog.txt start replace;
  set crit=like;
  exe sim500.nex;
  lscores 1;
  savetrees file=gtrig-alltrees.tre brlens;
  log stop;


  • Why is a neighbor joining tree generated? answer
  • What model is used? answer
  • What is the purpose of the second lset command? answer
  • What type of search will be conducted? answer
  • How many tree topologies will PAUP examine? answer

Now run the paup.nex file in PAUP*. Ordinarily you would create a script and feed that script to the qsub command to start the analysis, but it is possible to avoid creating a script by using the concept of a unix pipe - use echo to just write out the command you want to run and pipe this to qsub using the pipe symbol "|" (If you happen to have issued qlogin previously, you will first need to get out of the qlogin shell using Ctrl-d).




  • How is this tree different from the true tree? answer

Non-homogeneous models yield the true tree


/common/nh/monitor -f sim500


error while reading MCParameters : do not recognise the program version
read -1.-1 instead of 2.3


data file	sim500a.nex
Nsite		500
Ntaxa		8

save every	10
sample size	157

total time            : 	1	 hour, 	12	 min, 	11	 sec.
time per sample point : 	27.585	 seconds

StatInfCount		0






FYI: Things I had to do to get nhPhyloBayes to read my NEXUS data file

nhPhyloBayes is extremely picky about some aspects of the nexus file format. Here are three things I had to do in order to get nhPhyloBayes to read the sim500.nex data file:

  • Make all nucleotides uppercase letters (TextWrangler can help with this)
  • Should be "Matrix", not "matrix"
  • Should be "Begin data", not "begin data"

You will probably find more issues when you try to modify your own nexus files. One symptom that you haven't done enough to your file to make nhPhyloBayes happy is that nhPhyloBayes will say something like this (in the case of "matrix" vs "Matrix"):

expected to read "Matrix"
terminate called without an active exception 

or this (in the case of a lowercase nucleotide):

taxon 0  site  0	c
PhyloParameters::RegisterWithData : error in zip data making
error in data matrix : -2

FYI: How I generated sim500.nex

To (re)create the file sim500.nex used in this exercise, create the following python script (p4 is written as an extension to Python):

from p4 import *


read("( (C:0.4,(A:0.1,B:0.1):0.2):0.2, (D:0.1,E:0.1):0.4, ((F:0.1,G:0.1):0.2,H:0.4):0.2 );")
t = var.trees[0]
a = func.newEmptyAlignment(dataType='dna', taxNames=['A','B','C','D','E','F','G','H'], length=500) = Data([a])
x = t.newComp(free=1, spec='specified', symbol='A', val=[0.4, 0.1, 0.1, 0.4])
y = t.newComp(free=1, spec='specified', symbol='B', val=[0.01, 0.49, 0.49, 0.01])
t.setModelThing(x, node=t.root, clade=1)
t.setModelThing(y, node=2, clade=1)
t.setModelThing(y, node=13, clade=1)
t.newRMatrix(free=1, spec='2p', val=4.0)
t.newGdasrv(free=1, val=0.5)
t.setPInvar(free=1, val=0.0)
t.simulate(calculatePatterns=True, resetSequences=True, resetNexusSetsConstantMask=True)"sim500.nex", writeDataBlock=1, interleave=0, flat=1, append=0)

The p4 manual describes most of the commands in this file, so I will only go into detail on a few key commands.

  • The read command provides a means for creating the true tree used as the basis for simulating data
  • The line defining a creates an empty data matrix that we will later fill with simulated data
  • The lines defining x and y create new compositions (sets of nucleotide frequencies) which we will assign to different sections of the tree before simulating data
  • The setModelThing lines assign composition x to most of the tree and composition y to just 2 edges (you can use t.draw() to find out what numbers are assigned to each node)
  • The simulate command actually simulates the data
  • The writeNexus command stores the simulated data in the file sim500.nex
  • The last line shows a summary of the nucleotide composition of each taxon

Simulating just one data set takes so little time that you can just run this file directly - no need to use qsub:


Here, is the name of the Python script you just created, and note that I used the full path to python. The reason that the full path is needed is that there is more than one python version installed on the cluster (this one is 2.7.3, while the default one is 2.6.6), and only version 2.7.3 has p4 installed. In fact, if you want to always use this version (a good idea), add this line to the bottom of your ~/.bashrc file:

export PATH="/opt/python/bin:$PATH"

The ~/.bashrc file is read each time you login, and the line we added causes /opt/python/bin to be inserted at the beginning of your PATH, which is an environmental variable that tells the operating system where to look for programs that you ask it to run. Now when you ask for python, it will find the python in /opt/python/bin before it finds the default version, which is in /usr/bin. You can add the line above either by editing ~/.bashrc in nano, or you can append the line using the cat command (be sure that the hypthen and ">>" each are separated by spaces, and there is no space between the two ">" symbols):

cat - >> ~/.bashrc

After entering the first line above, the cat command will sit and wait for you to type the second line. After entering the second line, press (don't type out) Ctrl-d - this is a signal to the cat command that you have finished entering the text that you want to append to the file ~/.bashrc. The symbol >> means "append" whereas just ">" would mean "overwrite", so be careful to use ">>" if you go this route.

Literature Cited





Personal tools