Tutorial 4: Evaluating the accuracy of a multiple sequence alignment

Objectives

  • Complete the evaluating the accuracy of a multiple sequence alignment workflow.

  • Interpret the reported values and generated artifacts in their scientific context.

  • Identify the canonical command references for each analysis step.

Prerequisites and working directory

Install the current PhyKIT release and create a dedicated working directory. Download the data linked in this tutorial into that directory before running the commands. All paths below are relative to this directory.

mkdir phykit-tutorial-04
cd phykit-tutorial-04

Workflow

Evaluating the accuracy of multiple sequence alignments is an appropriate way to benchmark multiple sequence alignment strategies. Two popular methods to assess multiple sequence alignment accuracy are sum-of-pairs score and column score, which were introduced by Thompson et al., Nucleic Acids Research (1999), doi: 10.1093/nar/27.13.2682. Sum-of-pairs is calculated by summing the correctly aligned residue pairs over all pairs of sequences. Column score is calculated by summing the correctly aligned columns over all columns in an alignment. Both metrics range from 0 to 1 and higher values indicate more accurate alignments. Correctly aligned pairs or columns require knowing some ground truth of what the correct alignment is. Thus, a reference alignment that is perfectly (or near-perfectly) aligned is required. A reference alignment can be generated using simulations or be obtained from publicly available databases such as BAliBASE 4 (http://www.lbgi.fr/balibase/). For this tutorial we will use a reference alignment from BAliBASE.

Download test data: msa_accuracy.tar.gz


Step 0: Generate query alignments

In the msa_accuracy directory, there are two fasta files: BBA0001_query.faa, an unaligned set of sequences and BBA0001_reference.faa, the reference alignment. We will align BBA0001_query.faa using three different strategies implemented in Mafft, v.7.475 (https://mafft.cbrc.jp/alignment/software/), and evaluate the accuracy of each strategy.

To do so, please ensure Mafft is installed and then execute the following commands:

# first alignment
mafft --localpair BBA0001_query.faa > BBA0001_query.localpair.faa

# second alignment
mafft --genafpair BBA0001_query.faa > BBA0001_query.genafpair.faa

# third alignment
mafft --globalpair BBA0001_query.faa > BBA0001_query.globalpair.faa

To gain some initial insight as to whether the alignments differ, we can look at the length of each alignment using the aln_len function

for i in $(ls *pair.faa) ; do phykit aln_len $i ; done

1560
1464
1497

However, alignment length is not a measurement of accuracy. Thus, we will score each alignment using the sum_of_pairs_score and column_score functions.


Step 1: Score each alignment

For both functions, the first argument is the query alignment and the -r/--reference argument specifies the reference alignment. We will programmatically score each alignment using the same for loop that was used to calculate alignment length.

echo -e "BAliBASE_id_and_aln_strategy\tsop\tcs"
for i in $(ls *pair.faa)
do
   sop=$(phykit sum_of_pairs_score $i -r BBA0001_reference.faa)
   cs=$(phykit column_score $i -r BBA0001_reference.faa)
   echo -e "$i\t$sop\t$cs"
done

BAliBASE_id_and_aln_strategy sop     cs
BBA0001_query.genafpair.faa  0.8964  0.2943
BBA0001_query.globalpair.faa 0.9025  0.292
BBA0001_query.localpair.faa  0.8992  0.2959

Examination of the output reveals that the globalpair strategy has more correctly aligned pairs because it has a higher sum-of-pairs score whereas the localpair strategy has more correctly aligned columns. Of note, the column scores are generally low, which reflects a potential limitation of column score wherein column score is sensitive to alignment errors.

In summary, calculating sum-of-pairs score and column score can help assess the accuracy of multiple sequence alignment strategies.


Expected artifacts

Each step identifies its expected terminal output or generated files. Confirm that those artifacts exist before continuing to the next step; filenames are relative to the tutorial working directory unless an absolute path is shown.

Troubleshooting

  • Run phykit <command> --help to compare an invocation with the live interface.

  • Confirm that downloaded files are in the current working directory and retain the filenames shown in the tutorial.

  • For parsing errors, compare taxon names exactly across alignments, trees, and trait tables, including capitalization and underscores.

  • See Troubleshooting for installation, format, and error-reporting guidance.