BAli-Phy Tutorial (for version 4.0-beta11)

% cd ~/alignment_files/examples/5S-rRNA

% bali-phy 5d-clustalw.fasta --smodel 'gtr +> Rates.gamma(4) +> inv' --name 5d-freeA &

Yay! Now you have started your first copy of the analysis. BAli-Phy will create a directory called 5d-freeA-1 to hold output files from this analysis, as directed by the --name 5d-freeA option.

Let's run another copy of the same analysis:

% bali-phy 5d-clustalw.fasta --smodel 'gtr +> Rates.gamma(4) +> inv' --name 5d-freeA &

This command which will create another directory called 5d-freeA-2 for its own output files. This additional run will take advantage of a second processor. It will also allow us to determine if the two runs have performed enough iterations to agree.

% jobs

% top          // press 'q' to exit top

% wc -l 5d-*/C1.log

% bp-analyze 5d-freeA-1/ 5d-freeA-2/
% firefox Results/index.html &

% wc -l 5d-*/C1.log
% bp-analyze 5d-freeA-1/ 5d-freeA-2/
% mv Results 5d-freeA.Results
% firefox 5d-freeA.Results/index.html &

% top         // use top to find the PID (process id) for each of the bali-phy processes
% kill pid1
% kill pid2
% top         // check that the jobs have stopped.

% figtree 5d-freeA.Results/c50.PP.tree &

% cd ~/alignment_files/examples
% bali-phy 5d-clustalw.fasta -S 'gtr +> Rates.gamma(4) +> inv' -I none -n 5d-fixedA &
% bali-phy 5d-clustalw.fasta -S 'gtr +> Rates.gamma(4) +> inv' -I none -n 5d-fixedA &

Summarize the analysis after about 2000 iterations:

% wc -l 5d-*/C1.log
% bp-analyze 5d-fixedA-1/ 5d-fixedA-2

Examine the majority consensus tree for the fixed alignment analysis:

% figtree Results/c50.PP.tree &

% killall bali-phy

% cd ~/alignment_files/examples/ferns/
% alignment-cat -c 5-8       -e cleaned.fasta > exon1.fasta
% alignment-cat -c 9-453     -e cleaned.fasta > intron1.fasta
% alignment-cat -c 454-596   -e cleaned.fasta > exon2.fasta
% alignment-cat -c 597-864   -e cleaned.fasta > intron2.fasta
% alignment-cat -c 865-948   -e cleaned.fasta > exon3.fasta
% alignment-cat -c 949-1328  -e cleaned.fasta > intron3.fasta
% alignment-cat -c 1329-1393 -e cleaned.fasta > exon4.fasta

% alignment-cat exon1.fasta intron1.fasta exon2.fasta intron2.fasta exon3.fasta intron3.fasta exon4.fasta > combined.fasta

We will add -S 'gtr+>Rates.free(n=4)' to use the free-rates model of substitution with 4 categories:

% bali-phy exon1.fasta intron1.fasta exon2.fasta intron2.fasta exon3.fasta intron3.fasta exon4.fasta -S 'gtr+>Rates.free(n=4)' --test

This model is more general than Rates.gamma, since it estimates 4 free rates and the frequencies of the 4 rates. These 8 parameters have 6 degrees of freedom, which is a lot more than the 1 degree of freedom for the Rates.gamma model. Furthermore, each partition by default gets a separate copy of the model, which leads to 42 degrees of freedom just for the Rates.free part of the model!

We will now add -I 1,3,5,7:none to disable alignment estimation, but only for the exons:

% bali-phy exon1.fasta intron1.fasta exon2.fasta intron2.fasta exon3.fasta intron3.fasta exon4.fasta -S 'gtr+>Rates.free(n=4)' -I 1,3,5,7:none --test

% bali-phy exon1.fasta intron1.fasta exon2.fasta intron2.fasta exon3.fasta intron3.fasta exon4.fasta -I 1,3,5,7:none -S '1,3,5,7:gtr+>Rates.free(n=4)' -S '2,4,6:gtr+>Rates.free(n=4)' --test

% bali-phy exon1.fasta intron1.fasta exon2.fasta intron2.fasta exon3.fasta intron3.fasta exon4.fasta -I 1,3,5,7:none -S '1,3,5,7:gtr+>Rates.free(n=4)' -S '2,4,6:gtr +> Rates.free(n=4)' -I 2,4,6: --test

Instead of creating separate files for each intron and exon, we can also refer to character sets of the combined file:

% bali-phy cleaned.fasta:5-8 cleaned.fasta:454-596 cleaned.fasta:865-948 cleaned.fasta:1329-1393 cleaned.fasta:9-453 cleaned.fasta:597-864 cleaned.fasta:949-1328 -I 1,2,3,4:none -S '1,2,3,4:gtr +> Rates.free(n=4)' -S '5,6,7:gtr +> Rates.free(n=4)' -I 5,6,7: --test

However, the command line is getting so long that it is difficult to manage.

Instead of writing all options on the command line, We can put some of them into a text file that is easier to manage and edit. Let's make a text file called options1.txt:

align = cleaned.fasta:5-8       # exon1
align = cleaned.fasta:9-453     # intron1
align = cleaned.fasta:454-596   # exon2
align = cleaned.fasta:597-864   # intron2
align = cleaned.fasta:865-948   # exon3
align = cleaned.fasta:949-1328  # intron3
align = cleaned.fasta:1329-1393 # exon4

# exons
smodel = 1,3,5,7:gtr +> Rates.free(n=4)
imodel = 1,3,5,7:none

# introns
smodel = 2,4,6:gtr +> Rates.free(n=4)
imodel = 2,4,6:

% bali-phy -c options1.txt --test

Let's first extract the coding data:

% alignment-cat exon{1,2,3,4}.fasta -c2-295 > coding.fasta     // Concatenate exons and chop off partial codons at beginning and end
% sed -i 's/-/N/g' coding.fasta                                // Change codons like A-- to ANN at end of sequence

We can then construct an option file options2.txt:

align = intron1.fasta
align = intron2.fasta
align = intron3.fasta
align = coding.fasta

smodel = 1,2,3:gtr +> Rates.free(n=4)
imodel = 1,2,3:rs07

smodel = 4:gy94
#smodel = 4:gy94_ext(gtr_sym)     # A GTR-based version of GY94 a.k.a. M0.
#smodel = 4:gtr +> fMutSel0
#smodel = 4:function(w:fMutSel0(omega=w)) +> m3
imodel = 4:none

% bali-phy -c options2.txt --test

% alignment-translate < coding.fasta > protein.fasta

We can then construct an option file options3.txt:

align = intron1.fasta
align = intron2.fasta
align = intron3.fasta
align = protein.fasta

smodel = 1,2,3:gtr +> Rates.free(n=4)
imodel = 1,2,3:rs07

smodel = 4:lg08 +> Rates.log_normal(n=4)
imodel = 4:none

% bali-phy -c options3.txt --test

When you don't specify values for parameters like imodel, bali-phy uses sensible defaults. For example, these two commands are equivalent:

% cd ~/alignment_files/examples/
% bali-phy 5S-rRNA/25-muscle.fasta --test
% bali-phy 5S-rRNA/25-muscle.fasta --test --alphabet=RNA --smodel=tn93 --imodel=rs07

You can change the substitution model from the Tamura-Nei model to the General Time-Reversible model:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S gtr

Here the -S gtr is a short form of --smodel=gtr, where smodel means the substitution model.

You can also allow different sites to evolve at 5 different rates using the gamma[4] + INV model of rate heterogeneity:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Rates.gamma(4) +> inv'

You can also use a log_normal instead of a gamma. log_normal+INV is sometimes better behaved than gamma+INV, because the smallest rate bin under log_normal is not quite as close to 0.

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Rates.log_normal(4) +> inv'

You can allow 5 different rates that are all independently estimated:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Rates.free(n=5)'

Heterotachy occurs when the evolutionary rate changes across the tree. One version of heterotachy is the Tuffley-Steel (1998) model, where each state alternatives between an ON and OFF state.

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Covarion.ts98'

The Huelsenbeck (2002) model allows combining rate variation with Tuffley-Steel heterotachy:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Rates.gamma +> Covarion.hb02'

You can also use the Huelsenbeck model with the free-rates model:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'gtr +> Rates.free +> Covarion.hb02'

We can also conduct codon-based analyses using the Nielsen and Yang (1998) model of diversifying positive selection (dN/dS):

% bali-phy Globins/bglobin.fasta --test -S 'gy94(pi=f1x4)'

The gy94_ext model extends the gy94 model by taking a nucleotide exchange model as a parameter. The gy94 model is equivalent to gy94_ext(hky85_sym). However, if you want a more flexible codon model, you could use gtr_sym:

% bali-phy Globins/bglobin.fasta --test -S 'gy94_ext(gtr_sym,pi=f1x4)'

You can make the codon frequencies to be generated from a single set of nucleotide frequencies:

% bali-phy Globins/bglobin.fasta --test -S 'mg94_ext(gtr)'

The M7 model allows different sites to have different dN/dS values, where the probability of dN/dS values follows a beta distribution:

% bali-phy Globins/bglobin.fasta --test -S m7

The M7 model has parameters as well. Here are the defaults:

% bali-phy Globins/bglobin.fasta --test -S 'function(w:gy94(omega=w)) +> m7'

The M3 model allows different sites to have different dN/dS values, but directly estimates what these values are:

% bali-phy Globins/bglobin.fasta --test -S 'm3(n=3)'

The M8a_Test model allows testing for positive selection in some fraction of the sites:

% bali-phy Globins/bglobin.fasta --test -S 'function(w:gy94(omega=w,pi=f3x4)) +> m8a_test'

We can use the TN93+Gamma(4)+INV model without specifying parameters:

% bali-phy Globins/bglobin.fasta --test -S 'tn93 +> Rates.gamma +> inv'

However, we can also fix parameter values:

% bali-phy Globins/bglobin.fasta --test -S 'tn93 +> Rates.gamma(n=4,alpha=1) +> inv(p_inv=0.2)'

Here we have set the shape parameter for the Gamma distribution to 1, and the fraction of invariant sites to 20%. Since these parameters are fixed, they will not be estimated and their values will not be shown in the log file.

You can see the parameters for a model by using the help command, as in:

% bali-phy help Rates.gamma

This will show the default value or default prior for each parameter, if there is one.

By default the fraction of invariant sites follows a uniform(0,1) distribution:

% bali-phy help inv

However, we can specify an alternative prior:

% bali-phy Globins/bglobin.fasta --test -S 'tn93 +> Rates.gamma(n=4) +> inv(p_inv~uniform(0,0.2))'

We can also specify parameters as positional arguments instead of using variable names:

% bali-phy Globins/bglobin.fasta --test -S 'tn93 +> Rates.gamma(4) +> inv(~uniform(0,0.2))'

Here "~" indicates a sample from the uniform distribution instead of the distribution itself.

The insertion-deletion model also has parameters.

% bali-phy help rs07

Here the default value for rs07:mean_length is exponential(10,1). This indicates a random value that is obtained by sampling an Exponential random variable with mean 10 and then adding 1 to it.

Now lets specify different substitution models for different partitions.

% cd ~/alignment_files/examples/ITS
% bali-phy {ITS1,5.8S,ITS2}.fasta -S 1:gtr -S 2:hky85 -S 3:tn93 --test

We can also disable alignment estimation for some, but not all, partitions:

% bali-phy {ITS1,5.8S,ITS2}.fasta -I 1:rs07 -I 2:none -I 3:rs07 --test

Specifying -I none removes the insertion-deletion model and parameters for partition 2 and also disables alignment estimation for that partition.

Note that there is no longer an I3 indel model. Partition #3 now has the I2 indel model.

We can also specify that some partitions with the same model also share the same parameters for that model:

% bali-phy {ITS1,5.8S,ITS2}.fasta -S 1,3:gtr -I 1,3:rs07 -S 2:tn93 -I 2:none --test

This means that the information is pooled between the partitions to better estimate the shared parameters.

Take a look at the model parameters, and the parentheticals after the model descriptions. You should see that there is no longer an S3 substitution model or an I3 indel model. Instead, partitions #1 and #3 share the S1 substitution model and the I1 indel model.

We can also specify that some partitions share the same scaling factor for branch lengths:

% bali-phy {ITS1,5.8S,ITS2}.fasta -S 1,3:gtr -I 1,3:rs07 -S 2:tn93 -I 2:none --scale=1,3: --test

This means that the branch lengths for partitions 1 and 3 are the same, instead of being independently estimated.

Take a look at the model parameters. There is no longer a Scale[3] parameter. Instead, partitions 1 and 3 share Scale[1].

You might want to split concatenated gene regions:

% cd ~/alignment_files/examples/ITS-many/
% alignment-cat -c1-223 ITS-region.fasta > 1.fasta
% alignment-cat -c224-379 ITS-region.fasta > 2.fasta
% alignment-cat -c378-551 ITS-region.fasta > 3.fasta

% alignment-cat 1.fasta 2.fasta 3.fasta > 123.fasta

You might want to reduce the number of taxa while attempting to preserve phylogenetic diversity:

% alignment-thin --cutoff=5 ITS-region.fasta -v --preserve=Csaxicola420 > ITS-region-thinned.fasta

This removes sequences with 5 or fewer differences to the closest other sequence.

You can also specify that data sets should be shrunk to a specific number of sequences:

% alignment-thin --down-to=30 ITS-region.fasta -v > ITS-region-thinned.fasta

To see all options, run:

% alignment-thin --help

% alignment-thin --longer-than=250 ITS-region.fasta > ITS-region-long.fasta

% alignment-thin --remove-crazy=10 ITS-region.fasta > ITS-region-sane.fasta

% alignment-thin --min-letters=5 ITS-region.fasta > ITS-region-censored.fasta