BAli-Phy Tutorial (for version 4.3)

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

% bali-phy 5d-clustalw.fasta --smodel 'GTR +> ASRV.Gamma(4) +> Inv' --name 5d-freeA &

You have started the 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 +> ASRV.Gamma(4) +> Inv' --name 5d-freeA &

This command 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-freeA-*/C1.log

% bpy-summarize 5d-freeA-1/ 5d-freeA-2/

Open Results/index.html in your web browser.

% wc -l 5d-freeA-*/C1.log
% bpy-summarize 5d-freeA-1/ 5d-freeA-2/
% mv Results 5d-freeA.Results

Open 5d-freeA.Results/index.html in your web browser.

% 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/5S-rRNA
% bali-phy 5d-clustalw.fasta -S 'GTR +> ASRV.Gamma(4) +> Inv' -I none -n 5d-fixedA &
% bali-phy 5d-clustalw.fasta -S 'GTR +> ASRV.Gamma(4) +> Inv' -I none -n 5d-fixedA &

Summarize the results while the analyses run. For this exercise, collect about 2000 iterations from each run before comparing the fixed and free alignments:

% wc -l 5d-fixedA-*/C1.log
% bpy-summarize 5d-fixedA-1/ 5d-fixedA-2

Examine the majority consensus tree for the fixed alignment analysis:

% figtree Results/c50.PP.tree &

% 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 +> ASRV.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 +> ASRV.Free(n=4)' --test

This model is more general than ASRV.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 ASRV.Gamma model. Furthermore, each partition by default gets a separate copy of the model, which leads to 42 degrees of freedom just for the ASRV.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 +> ASRV.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 +> ASRV.Free(n=4)' -S '2,4,6:GTR +> ASRV.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 +> ASRV.Free(n=4)' -S '2,4,6:GTR +> ASRV.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 +> ASRV.Free(n=4)' -S '5,6,7:GTR +> ASRV.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:

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

# exons
:smodel 1,3,5,7:GTR +> ASRV.Free(n=4)
:imodel 1,3,5,7:none

# introns
:smodel 2,4,6:GTR +> ASRV.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 '/^>/!s/-/N/g' coding.fasta > coding-with-N.fasta        # Replace gaps in sequences, preserving FASTA names
% mv coding-with-N.fasta coding.fasta

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 +> ASRV.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:|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 +> ASRV.Free(n=4)
:imodel 1,2,3:RS07

:smodel 4:LG +> ASRV.LogNormal(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 ASRV.Gamma(4) + Inv model of rate heterogeneity:

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'GTR +> ASRV.Gamma(4) +> Inv'

You can also use ASRV.LogNormal instead of ASRV.Gamma. ASRV.LogNormal+Inv is sometimes better behaved than ASRV.Gamma+Inv, because the smallest rate bin under ASRV.LogNormal is not quite as close to 0.

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'GTR +> ASRV.LogNormal(4) +> Inv'

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

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'GTR +> ASRV.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 alternates between an ON and OFF state.

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

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

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'GTR +> ASRV.Gamma +> Covarion.Huelsenbeck02'

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

% bali-phy 5S-rRNA/25-muscle.fasta --test -S 'GTR +> ASRV.Free +> Covarion.Huelsenbeck02'

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 '|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 '|w:GY94(omega=w,pi=F3x4)| +> M8a_test'

See the User Guide's positive-selection section for model choices, running analyses, and interpreting posterior support.

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

% bali-phy Globins/bglobin.fasta --test -S 'TN93 +> ASRV.Gamma +> Inv'

However, we can also fix parameter values:

% bali-phy Globins/bglobin.fasta --test -S 'TN93 +> ASRV.Gamma(n=4,alpha=1) +> Inv(pInv=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 ASRV.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 +> ASRV.Gamma(n=4) +> Inv(pInv~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 +> ASRV.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:meanLength is ShiftedExponential(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 let's 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 2: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 scale3 parameter. Instead, partitions 1 and 3 share scale1.

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 -c380-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 --verbose --protect=Csaxicola420 > ITS-region-thinned.fasta

This removes sequences with fewer than 5 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 --verbose > 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-gappy=10 ITS-region.fasta > ITS-region-sane.fasta

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