Multiple Sequence Alignments

Exercises and Demonstrations

Aidan Budd

3D Structure-Based Interpretation of Alignments

Teaching Objectives

After completing this section, you will hopefully appreciate the properties that we expect residues in the same column to share of a sequence alignment that is being interpreted in a "structural" context.

Notes

Comparing (and aligning) pairs of similar structures demonstrates what it means for a pair of residues to be "structurally equivalent" - the relationship that we want residues in the same column of an alignment to share if we are using the alignment in a "structural" context e.g. predicting secondary structure, or building a protein profile HMM to carry out a sensitive sequence similarity search.

At the same time, it demonstrates that there may be regions of two structures for which there is not any such equivalence.

Demonstration

This example is a pair of bacterial toxins 1ji6 and 1i5p with very similar structures

We have aligned the N-terminal regions of these structures using FATCAT

this link is to the aligned structures in a PyMOL session file

two examples of structurally equivalent residues are coloured in dark blue
side-chains of adjacent residues that do not share 1:1 structural equivalence with these residues are shown in the default colours of their protein chains
regions of the chains which do not have any 1:1 structural equivalence in the other structure are shown in lighter colours

notice that these are all in surface loops

this link is the sequence alignment implied by this structural alignment, with residues coloured as in the above PyMOL session file
this link provides the same alignment, but in Jalview format (allows us to easily change the format of regions of the alignment, and to examine the effect of introducing gaps in different positions in the alignment)

It should be clear, when looking at the structural alignment, that the structures are very similar, with most residues sharing 1:1 structural equivalence with a residue in the other structure.

By contrast, we provide below an alignment of two very different structures. Indeed, considered in terms of the kinds of secondary structure elements they contain, they are completely different (one is mainly alpha, the other mainly beta).

1PPR, chain M, alpha-solenoid
1GVK, chain B, beta barrel
PyMOL session file of pairwise structural alignment using CE pairwise structural alignment server
Pairwise sequence alignment indicated from above structural alignment

Note that we are still able to align these two structures, despite our opinion that they are global extremely dissimilar. The same is true for multiple sequences alignments - most software will report an alignment, whether or not the global similarity most of the software assumes the sequences share is indeed present within the sequences.

Given the global dissimilarity of the alpha solenoid and beta barrel structures, you would almost certainly want to avoid using such an alignment to make any inferences about similarity of function/structure between residues aligned in the same column.

In general, this illustrates the fact that it is important to be very confident that the sequences you include in an MSA indeed share the relationship you are interested in - this will typically be "structural" and/or "evolutionary" equivalence. Note, however, that we are avoiding a discussion of how to judge whether or not two sequences do indeed have such a relationship with each other - this involves using tools such as BLAST, that we we are not focusing on in this set of exercises.

If you are interested, you might like to try comparing some other pairs of structures and examining their structural and sequence alignments as calculated by FATCAT and CE.

Instructions on how to calculate pairwise structural alignments using CE and FATCAT and display the results in PyMOL

Note that all the structures used in the exercise below contain only one polypeptide chain, chain A. You need to supply this information when submitting the structures to CE and FATCAT.

Serine Proteases

1ton PDBsum Website PDB at RCSB Website
5ptp PDBsum Website PDB at RCSB Website
If you have trouble using the servers, here are pre-calculated results

FATCAT structure alignment PyMOL session file and sequence alignment
CE structure alignment PyMOL session file and sequence alignment

TIM Barrel hydrolases

1bf2 PDBsum Website PDB at RCSB Website
1smd PDBsum Website PDB at RCSB Website
If you have trouble using the servers, here are pre-calculated results

FATCAT structure alignment PyMOL session file and sequence alignment
CE structure alignment PyMOL session file and sequence alignment

Note that we are suggesting that you use FATCAT and CE as they both provide a pairwise sequence alignment along with easy-to-view structural alignments - not because they are known to provide on average the best quality alignments (although these seem to be often good).

Non-Equivalence Between Structural-Based and "Evolutionary"-/"Homology"-Based Interpretation of Alignments

Teaching Objectives

After completing this section, you will hopefully appreciate that while we often assume that "evolutionarily equivalent"/"homologous" residues are also structurally equivalent, this is not necessarily the case. This should highlight the importance of being aware of which of these contexts you are assuming when working with any given sequence alignment

Notes

During the accompanying presentation, you will probably have carried out an exercise discussing an alternative way of interpreting sequence alignments.

Remember that the previous section describes interpretation of alignments in a 3D-structural context - residues in the same column of an alignment should ideally share structurally-equivalent locations within the 3D structures of the molecules they describe. In the alternative "evolutionary" context, residues in the same column of an alignment should be able to trace their evolutionary history back to the same residues in the most recent common ancestor of the two sequences - it is assumed that any differences observed in a column (e.g. a column in which two sequences have a different amino acid residue e.g. glycine in one sequence and alanine in another) are due only to point mutation/substitution events.

In many cases it is appropriate to assume that closely-related sequences have very similar structures - in the previous section we looked at alignments of 3D protein structures between proteins that were evolutionarily related.

However, structural equivalence does not necessarily imply evolutionary relatedness - nor does evolutionary relatedness necessarily imply structural equivalence, as we will see in the following demonstration.

Notes - Short Protein Linear Motifs

Short protein linear motifs, predominantly found in eukaryotic proteomes (hence sometimes called ELMs for Eukaryotic Linear Motifs - follow this link to a resource developed by our group for exploring the presence of ELMs in your sequence of interest) are typically between 3 - 9 amino acid residues in length, and have many important roles in a wide range of different cellular processes. They typically function by mediating interactions with specific globular domains e.g. protein kinase motifs interact with protein kinase domains, cyclin-binding boxes interact specifically with cyclins etc.

While there is typically one or more position in a motif that must be present for the motif to be functional, other positions can often be occupied by many different possible amino acid residues. For example, a motif that is able to bind the SH3 domain has the sequence pattern ...[PV]..P where "." indicates any one of the 20 possible amino acids, "[PV]" indicates either the amino acids proline or valine, and "P" indicates that the position must be occupied by a proline.

The relatively weak constraints placed on the sequence of such a motif make them relatively likely to appear de novo in a protein sequence, for example, as the result of point mutations - in stark contrast to globular domains. Thus, where the same linear motif is present in two unrelated proteins (i.e. we have no evidence for the relatedness of the proteins to each other, but both possess a short sequence mediating an interaction with the same kind of globular domain), we assume that the two motifs are not evolutionarily related to each other. Thus, it would not make sense to analyse an alignment of such motifs in an evolutionary context - despite the fact that the 3D structures of the motifs will usually be very similar.

Demonstration - Structural Equivalence Without "Evolutionary" Equivalence

Both the HIV1 nef protein (UniProt:NEF_HV1BR) and the human neutrophil cytosol factor 1 (UniProt:NCF1_HUMAN) contain SH3-domain binding motifs that match the pattern "...[PV]..P" (note that peptides matching the patterns "P..P.[KR]" and "[RKY]..P..P" are also known to bind to the SH3 domain). These SH3-binding peptides are highlighted in red in the sequences shown below.

>NEF_HV1BR Nef
MGGKWSKSSVVGWPTVRERMRRAEPAADGVGAASRDLEKHGAITSSNTAATNAACAWLEA
QEEEEVGFPVTPQVPLRPMTYKAAVDLSHFLKEKGGLEGLIHSQRRQDILDLWIYHTQGY
FPDWQNYTPGPGVRYPLTFGWCYKLVPVEPDKVEEANKGENTSLLHPVSLHGMDDPEREV
LEWRFDSRLAFHHVARELHPEYFKNC

>NCF1_HUMAN Neutrophil cytosol factor 1
MGDTFIRHIALLGFEKRFVPSQHYVYMFLVKWQDLSEKVVYRRFTEIYEFHKTLKEMFPI
EAGAINPENRIIPHLPAPKWFDGQRAAENRQGTLTEYCSTLMSLPTKISRCPHLLDFFKV
RPDDLKLPTDNQTKKPETYLMPKDGKSTATDITGPIILQTYRAIADYEKTSGSEMALSTG
DVVEVVEKSESGWWFCQMKAKRGWIPASFLEPLDSPDETEDPEPNYAGEPYVAIKAYTAV
EGDEVSLLEGEAVEVIHKLLDGWWVIRKDDVTGYFPSMYLQKSGQDVSQAQRQIKRGAPP
RRSSIRNAHSIHQRSRKRLSQDAYRRNSVRFLQQRRRQARPGPQSPGSPLEEERQTQRSK
PQPAVPPRPSADLILNRCSESTKRKLASAV

The following PyMOL session files show the 3D protein structures of these SH3-motif:SH3-domain intereactions. SH3 domains are coloured in cyan, the SH3-binding motifs are coloured in pink/salmon, and are shown with their amino acid side chains drawn as sticks. The N-termini of the peptides are coloured in blue, the C-termini in red.

NCF1_HUMAN interacting with an SH3 domain also contained within the NCF1_HUMAN protein

Original PDB file (PDB:1w70)

local copy
link to PDB record at the RCSB
link to PDBsum record at the EBI

PyMOL session file
Multiple sequence alignemnt of nef proteins (via BLAST at the NCBI)

NEF_HV1BR interacting with the SH3 domain in human FYN1 (UniProt:FYN1_HUMAN)

Original PDB file (PDB:1efn):

local copy
link to PDB record at the RCSB
link to PDBsum record at the EBI

PyMOL session file
Multiple sequence alignment of NCF1-related proteins (very roughly adapted from the Ensembl GeneTree that includes human NCF1):

This PyMOL session file shows the two motifs aligned in 3D - there are clear 1:1 equivalencies between positions in the two structures, despite no evidence for the two sequences being evolutionarily related (i.e. that their sequences can be traced back to the same ancestral protein sequence via either point mutaitons or no mutations at all).

Demonstration - One Sequence, Two Different Structures - Evolutionary Equivalence without Stuctural Equivalence

Ideally, we would show here an example of naturally-occurring sequences that were clearly related, but whose structures had been determined and found to form globular domains with different folds - but, as yet, we have not found such an example.

However, we do have an example of a sequence - the human Lymphotactin protein (UniProt:XCL1_HUMAN,P47992) - which adopts two different folds depending on conditions. The two sequences are clearly related (they're identical!), but in different contexts, evolutionarily-related/equivalent/homologous residues clearly have different structrues.

This PyMOL session file, which includes data from PDB records 1j8i (local file, PDB, PDBsum) and 2jp1 (local file, PDB, PDBsum), shows the two different forms of the protein in one file.

In the two structures, some regions form equivalent secondary structure elements - however, in at least two cases, the same residues form rather different structures (two such examples are highlighted in cyan and yellow in the PyMOL session file - regions of the same colour are equivalent in the two different structures).

Introducing Pairwise/Multiple Sequence Alignments with JalView

Demonstration

To explore our own ideas about what constitutes a good alignment, here we will carry out 'manual' alignments of some pairs of sequences using JalView. Most of us, when presented with a set of sequences, will have some 'intiuitive'/subconcious feeling about what a good or a bad alignment of them would look like. Trying to build your own 'good' alignments manually should help you understand better what you think makes an alignment 'good' or 'bad'.

This exercise is also useful as it gives you an opportunity to gain familiarity with JalView, which is (in my opinion) a good choice of software for editing sequence alignments - certainly it's what I use.

For instructions on how to use JalView follow this link.

I'll demo for you now:

Loading sequences into JalView - can be done in several different ways
Cutting out sequences from an alignment and saving them in a separate alignment
Changing the way the sequences are displayed - ClustalX colouring scheme, % Identity, etc.
Editing the alignment
Copying sequences in JalView and pasting them to another application

With JalView I will:

open this sequence file (containing four actin sequences)
take out (and save for later) all but two of the sequences
change the way the sequences are displayed - ClustalX colouring scheme, % Identity, etc.
manually adjust the alignment
save the results of the manual alignment

Exercise

Try the same with the sets of sequences provided at the bottom of this exercise - don't try and do all possible pairwise alignments, just pick out a few from each file.

While you're carrying out this exercise, make notes/discuss with your neighbours on the following issues/questions:

How do I decide where to put gaps in the alignment?
Are there regions of the alignment I think are 'better' aligned than others?

If so, what differentiates these different kinds of regions from each other?

Are some sequences harder/slower to align well than others?

If so, what features do the more 'difficult' alignements/sequences have which differentiate them from easier to align sequences?

Sequences:

The sequence sets linked to below are (1-3) three different, non-overlapping regions of the human, rat, and mouse collagen 18s and (4) full length

To check your answers, you could align the sequences automatically using the EBI MUSCLE webserver, and then try comparing these results with the alignments you have built manually.

"Manual" Pairwise Sequence Alignments using JalView - All Files Containing only 2 Sequences

Demonstration

We will demonstrate aligning these two tublin sequences to each other using JalView. The demonstration will involve

Loading the sequences into JalView using both

File->Input Alignment->from File
File->Input Alignment->from Textbox

Changing the residue colouring scheme to:

Percentage Identity
Clustalx

Inserting gaps into a sequence using the mouse
Saving your work

Note that choosing/building an alignment is "simply" about deciding where gaps belong and where they do not belong.

We will also automatically build pairwise alignments of the sequences using three different approaches:

and compare the results of these approaches to the alignment we bulit "manually".

Exercises

Try building pairwise alignments yourself using JalView in a similar way with the pairs of sequences found in the following files.

Try also one (or more) of the automatic pairwise alignment tools listed above.

Compare your "manual" alignment with the automatic alignments and:

identify any differences between them
try to decide which (the manual or one of the automatic alignments) you think is better

While carrying out the manual alignments, write down:

Features that describe a relatively "good" alignment, thinking in terms of

sizes of gaps
numbers of gaps
properties of residues in the same column (the same as each other? different?)

Instructions on how to change a "bad" alignment into a better one
Characteristics of sequences that are more difficult/take more time to align than others

Below are the sequences - if possible, try all of them, as they have been chosen to illustrate a range of different issues/points, and it will hopefully be useful for you to have encountered all of these.

The same two tubulin sequences used in the demonstration above
the coding regions of the cDNA for the same two tubulin sequences - note that to do the automatic alignments here you'll need to choose the nucleotide implementations of the tools
fragments of mouse and rat collagen 18s
ACTB_CERPY and ACTB_TRIVU
SRC_HUMAN and SRK3_SPOLA kinase domains
ARPM1_MOUSE and ACTB_TRIVU
SRC-like Kinases SRC_MOUSE and STK_HYDAT
full length mouse and rat collagen 18s (for those who fancy a challenge)

"Manual" Pairwise Sequence Alignments using JalView

Demonstration

To demonstrate some of the charactistics of sequences that can make them easier or more difficult to align, we will try aligning some sequences ourselves "manually" using the JalView sequence editor.

This exercise is also useful as it gives you an opportunity to gain familiarity with JalView, which is (in my opinion) a good choice of software for editing sequence alignments - certainly it's what I use.

For instructions on how to use JalView follow this link.

I'll demo for you now:

Loading sequences into JalView
Cutting out sequences from an alignment and saving them in a separate alignment
Editing the alignment
Removing all gaps
Removing all empty columns
Copying sequences in JalView and pasting them to another application

With JalView I will:

open this sequence file (containing sequence related to human alpha-2-macroglobulin)
take out (and save for later) all but two of the sequences
manually adjust the alignment
save the results of the manual alignment
compare the results of the alignment to that calculated using smith-waterman and needleman-wunsch algorithms using the EBI EMBOSS website
compare also to the results obtained using BLAST 2 Sequences at the NCBI
do the same with a different pair of alpha-2-macroglobulin-related sequences

I'll also try and make an MSA of all the sequences - and compare that with that obtained using the EBI MUSCLE webserver

Exercise

1. Protein Kinase domains

Load these protein kinase sequences into JalView

Use JalView to make pairwise alignments "manually" as for the above demonstration, and compare the results with those obtained using the smith-waterman, needleman-wunsch, and BLAST 2 Sequence algorithms

Also try to build an MSA using JalView from the same sequences, and compare the result to that obtained using the EBI MUSCLE webserver

2. Collagen 18s

The sequence sets linked to below are (1-3) three different, non-overlapping regions of the human, rat, and mouse collagen 18s and (4) full length

For each set of sequences, use JalView to build a multiple sequence alignment, and compare the results you obtain with those using the EBI MUSCLE webserver.

Are there any obvious features of some of the sequences that make them easier/more difficult to align?

Below are examples of manually-annotated alignments for the different datasets

"Manual" MSA using JalView

Notes

We'll look at making our own multiple sequences alignments "manually" using JalView rather than automatically aligning them with MUSCLE, MAFFT, PRANK etc.

We'll look at two different sets of sequences - one we find easier, the other more difficult, to align manually.

The reason why we find these difficult to align are the same as the reasons why automatic software finds some alignments easier, others more difficult:

overall level of sequence divergence
size of gaps

Here are some links with information on using JalView

Demonstration

Exercise

Making an MSA is Easy?!

The act of taking a set of sequences, inputting them to an automatic MSA tool, and obtaining a result is usually trivial.

Demonstration:

Download this set of proteins in FASTA format that are very similar to the human src (UniProt: SRC_HUMAN, P12931) protein
Copy them into the EBI MUSCLE web server
Press submit
Enjoy the resulting MSA

Exercise:

To see how easy it is to do this, try the demonstration above, trying to carry it out as quickly as possible

Comparing Results from Different Automatic MSA Tools - Short Exercise

Exercise

Using the collagen sequences from the previous exercise to compare the MSAs built automatically using the following list of different MSA software

MUSCLE
CLUSTALW
PRANK
MAFFT
PROBCONS

You'll hopefully see that:

different software gives different results
there are regions where they all agree on the correct alignment

usually the regions where the sequences are more similar to each other

there are other regions where they tend to disagree

usually the regions around the gaps

Comparing Sequence and Structural Alignments

Demonstration

How good are different structural alignment tools compared with purely sequence-based alignments?

We'll make a comparison using the example of two serine protease sequences (shown below) - the two sequences are both taken from a "gold-standard" alignment reference dataset (test set 5ptp_ref1 from BaliBase2)

By examining the UniProt record, I have highlighted in red the active-site residues for 5ptp.

Note that the sequence below describe the proteins found in the PDB (not the UniProt) records - in this case they differ in that the PDB sequence is the same as that in UniProt except for removal of the signal peptide and activation hexapeptide.

>5ptp Bovine Cationic Trypsin TRY1_BOVIN (beta-trypsin form following Lys-23 autocatalytic cleavage)
IVGGYTCGANTVPYQVSLNSGYHFCGGSLINSQWVVSAAHCYKSGIQVRLGEDNINVVEG
NEQFISASKSIVHPSYNSNTLNNDIMLIKLKSAASLNSRVASISLPTSCASAGTQCLISG
WGNTKSSGTSYPDVLKCLKAPILSDSSCKSAYPGQITSNMFCAGYLEGGKDSCQGDSGGP
VVCSGKLQGIVSWGSGCAQKNKPGVYTKVCNYVSWIKQTIASN

>1ton Rat Tonin KLK2_RAT
IVGGYKCEKNSQPWQVAVINEYLCGGVLIDPSWVITAAHCYSNNYQVLLGRNNLFKDEPF
AQRRLVRQSFRHPDYIPLIVTNDTEQPVHDHSNDLMLLHLSEPADITGGVKVIDLPTKEP
KVGSTCLASGWGSTNPSEMVVSHDLQCVNIHLLSNEKCIETYKDNVTDVMLCAGEMEGGK
DTCAGDSGGPLICDGVLQGITSGGATPCAKPKTPAIYAKLIKFTSWIKKVMKENP

Initial PDB files 5ptp 1ton (chain A in both cases)

We will align

structures using

sequences using

Exercise

Carry out a similar exercise yourself with the following amylase sequences.

The three active-site residues in both sequences as described by their UniProt records have been labeled in red in the FASTA format sequences provided below.

>1smd AMY1_HUMAN
GRTSIVHLFEWRWVDIALECERYLAPKGFGGVQVSPPNENVAIHNPFRPWWERYQPVSYK
LCTRSGNEDEFRNMVTRCNNVGVRIYVDAVINHMCGNAVSAGTSSTCGSYFNPGSRDFPA
VPYSGWDFNDGKCKTGSGDIENYNDATQVRDCRLSGLLDLALGKDYVRSKIAEYMNHLID
IGVAGFRIDASKHMWPGDIKAILDKLHNLNSNWFPEGSKPFIYQEVIDLGGEPIKSSDYF
GNGRVTEFKYGAKLGTVIRKWNGEKMSYLKNWGEGWGFMPSDRALVFVDNHDNQRGHGAG
GASILTFWDARLYKMAVGFMLAHPYGFTRVMSSYRWPRYFENGKDVNDWVGPPNDNGVTK
EVTINPDTTCGNDWVCEHRWRQIRNMVNFRNVVDGQPFTNWYDNGSNQVAFGRGNRGFIV
FNNDDWTFSLTLQTGLPAGTYCDVISGDKINGNCTGIKIYVSDDGKAHFSISNSAEDPFI
AIHAESKL

>1bf2 ISOA_PSEAY
DVIYEVHVRGFTEQDTSIPAQYRGTYYGAGLKASYLASLGVTAVEFLPVQETQNDANDVV
PNSDANQNYWGYMTENYFSPDRRYAYNKAAGGPTAEFQAMVQAFHNAGIKVYMDVVYNHT
AEGGTWTSSDPTTATIYSWRGLDNATYYELTSGNQYFYDNTGIGANFNTYNTVAQNLIVD
SLAYWANTMGVDGFRFDLASVLGNSCLNGAYTASAPNCPNGGYNFDAADSNVAINRILRE
FTVRPAAGGSGLDLFAEPWAIGGNSYQLGGFPQGWSEWNGLFRDSLRQAQNELGSMTIYV
TQDANDFSGSSNLFQSSGRSPWNSINFIDVHDGMTLKDVYSCNGANNSQAWPYGPSDGGT
STNYSWDQGMSAGTGAAVDQRRAARTGMAFEMLSAGTPLMQGGDEYLRTLQCNNNAYNLD
SSANWLTYSWTTDQSNFYTFAQRLIAFRKAHPALRPSSWYSGSQLTWYQPSGAVADSNYW
NNTSNYAIAYAINGPSLGDSNSIYVAYNGWSSSVTFTLPAPPSGTQWYRVTDTCDWNDGA
STFVAPGSETLIGGAGTTYGQCGQSLLLLISK

Here are the PDB structure file 1smd 1bf2

As for the serine proteases, align the above (note that pre-calculated results are provided in case the servers are down - try and get these results without consulting these files, however, unless necessary)

pairwise structures using

pairwise sequences using

multiple sequences using PROBCONS (here are the unaligned sequences taken from a BaliBase record 1pamA_ref1 containing both sequences)

execute this from the command line e.g.

probcons 1pamA_ref1.reference.fasta > 1pamA_ref1.probcons.fasta

here is the result

When examining the structural alignments, it may help to label within PyMOL the active site residues for the different chains

here are some basic instructions on how to use PyMOL to:

show the sequences
highlight residues
change colour of residues

Do the sequence- and structure-based alignments successfully align the active-site residues from the two sequences?
Do gaps in the structural alignment tend to be associated with secondary structure elements (alpha helices or beta strands) or with loop regions?

Compatibility of Multiple Pairwise Alignments

Exercise

There are three different pairwise alignments that can be built between these three haemoglobin sequence.

Estimate these using one of these three tools:

Here are the pre-calculated results obtained using EMBOSS's needle program with default settings (which implements the Needleman-Wunsch algorithm):

Try and build an MSA of these sequences (in JalView) so that all the alignments between residues in different sequences described in the three pairwise alignments are present in this MSA.

This exercise should have illustrated that ometimes (often/almost-always if the sequences being compared are very different from each other!) the complete set of all pairwise alignments for a set of sequences will not all be 'compatible' with each other i.e. it is impossible to build a single MSA that contains all of the aligned-column relationships found in all pairwise alignments.

Finding Pre-Calculated Alignments

Demonstration

If you can find a pre-calculated alignment that already contains a set of sequences that is useful for addressing your problem, you may well be able to save lots of time. Using the databases/resources below, we'll look at how you can get access to alignments that you can download yourself and examine locally.

You could try this out using SRC_HUMAN

>SRC_HUMAN|P12931|Proto-oncogene tyrosine-protein kinase Src (EC 2.7.10.2) ENSG00000197122
GSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEP
KLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDW
WLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESE
TTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLC
HRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLK
PGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYL
RLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTA
RQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERG
YRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL

Alternatively, as some of the alignments will be rather large for this well-studied, well-conserved protein, one could try instead with human Sperm-associated antigen 7 (UniProt:SPAG7_HUMAN).

>sp|O75391|SPAG7_HUMAN Sperm-associated antigen 7 ENSG00000091640
MADLLGSILSSMEKPPSLGDQETRRKAREQAARLKKLQEQEKQQKVEFRKRMEKEVSDFI
QDSGQIKKKFQPMNKIERSILHDVVEVAGLTSFSFGEDDDCRYVMIFKKEFAPSDEELDS
YRRGEEWDPQKAEEKRKLKELAQRQEEEAAQQGPVVVSPASDYKDKYSHLIGKGAAKDAA
HMLQANKTYGCVPVANKRDTRSIEEAMNEIRAKKRLRQSGEELPPTS

Ensembl (protein families and/or gene tree)
TreeFam
PFAM
HOMSTRAD
OPTIC (Try using the amniotes - query with Ensembl ID

Exercise

Using these (and perhaps other alignment resources you might find) find a pre-calculated MSA that is will be a good starting place for building an alignment to investigate:

the evolution of FGF genes (e.g. FGF4_HUMAN) in the early vertebrate lineage
variation in the 3D structure of tyrosine kinase domains within the animals

Building Automatic MSAs

Clearly, an important step in obtaining an MSA is the use of automatic MSA software to build an MSA from an inital set of sequences.

Although, as already pointed out, if possible, it is usually a good idea to begin instead with a pre-calculated alignment that already contains many (ideally all!) the sequences needed for your analysis. This can save you a lot of time - not just the time it takes (i) to calculate the alignment, but also the time taken to (ii) select an appropriate set of sequences for the analysis - although the time saved in this way clearly depends on how well-suited the pre-calculated alignment is for your analysis.

Obtaining an Initial Set of Sequences for Automatic Alignment via BLAST

A typical way to build an MSA "from scratch" (i.e. begining with a single sequence of interest, and from there collecting a set of intial sequences from which to build an automatic MSA) is to collect a set of sequences that are evolutionarily related to your initial sequence using BLAST.

Sometimes this is trivial - not always, however, as an initial search may identify:

too many sequences

slow to calculate alignments etc.
provide too much information, making it difficult to get an overview of the alignment

too few sequences

resulting alignment doesn't contain sufficient diversity for your analysis

sequences that are unrelated to the query with highly-significant E-values - due to "broken" BLAST statistics

applications of MSAs are predicated on their including only sequences that are evolutionarily/3D structurally related (as appropriate) - if the statistics are "broken" in this way, then (unless the user notices the problem), any application of the MSA will contain potentially very serious systematic errors

Sequence features that can cause these kinds of problems include:

sequence length

other things being equal (which they aren't... but...), BLAST queries with longer sequences will identify larger numbers of related sequences

conservation/rate of evolution

sequences that evolve slowly, i.e. experience a low rate of substitution, retain significant similarity over greater time periods - thus, other things being equal, searches with more conserved sequences will identify larger numbers of related sequences

taxonomic distribution

in many cases, there are simply no sequence related to the query sequences in many organisms - the greater the number of organisms that possess relatives of a sequence, other things being equal, the greater the number of related sequences that will be detected by the BLAST search

level of interest in the sequence in the scientific comunity

some proteins/molecules/sequences are of particular interest to the scientific community - searches using such sequences as queries tend to produce (very) large numbers of matches. For example:

p53
proteins encoded by the HIV virus such as nef
commonly-used phylogenetic markers such as mitochondrial cytochrome b, this sequence here from a somewhat obscure rodent, the ebony akodont, to illustrate the diversity of organisms for which this gene has been sequence in

structure - BLAST is good at accurately identifying relatedness between globular protein domains - less good when working with intrinsically unstructured protein (IUP) sequence

the substitution marrices (and gap scores) used by BLAST to identify relatedness between the query sequence and sequences in the databases are calculated from alignments of globular protein domains. Thus, they are less good at identifying relatedness between IUP regions, which have rather different patterns of substitutions

Demonstration

Here we will attempt to use BLAST to collect a set of sequences to examine rodent phylogeny.

We will use the mouse cytochrome b protein sequence as a query at the NCBI BLAST server.

As expected, given that this protein is a popular phylogenetic marker, BLAST using default settings do not give us what we need i.e. far too many sequences for us to get an overview of the result.

To find a more appropriate set of sequences, we will try varying several different parameters of the search:

number of matches displayed

as all the 100 matches shown by default are very similar to the query, it seems likely that the 101th etc. matches would also be very similar - so if we want to sample widely across lots of different rodents, we'll probably need to increase the number of matches shown

database

we can choose a smaller database, which is clearly likely to contain fewer sequences than a larger one (SwissProt is smaller than "nr")

filtering by Entrez query

as we're only interested in rodent sequences, we can use an Entrez query to restrict the search space - shortening the list of matches by removing all sequences that are not relevant for our analysis
Here's a link to a description of how to query different database fields using the Entrez system

Exercise

Use BLAST to collect a set of protein sequences related to human p53 protein, with sequences taken from the zebrafish (Danio rerio), mouse (Mus musculus), chicken (Galllus gallus), and african clawed frog (Xenopus laevis). Try to obtain an as-complete-as-possible list of such sequecnes from each organism.

Comparing Automatic MSA Tools

We have included an exercise/demonstration/session on comparing MSAs calculated by several different alignment tools for several reasons:

It provides experience of the user interface to several different tools
It demonstrates the kind of variation you might find in the results obtained by different MSA software
We often want to compare the MSAs built by different tools when attempting to build our own MSAs

In many cases, the MSAs caluclated a range of different tools will all have regions we suspect are mis-aligned. In this case, we need to inspect the resulting MSAs to find the one that looks best - which is the alignment we would then choose as our starting place for further work. Alternatively, we might decide that all the MSA are so bad that we need to build new automatic alignments using different parameters, hoping this improves their quality

Building/Examining MSAs with All Sequences in the Same Order

It may seem trivial - however, anyone who has tried to compare alignments of several similar sequences where the sequences are NOT all in the same order will know how difficult this can make the comparison. Without the same sequence order in all alignments, it can be very difficult to identify regions where the alignments differ.

Demonstration

There are several different ways to obtain obtain alignments with all sequences in the same order. We will demonstrate these using this set of sequences related to human p53

We will align the sequences using four different automatic MSA software packages - in all cases by simply cutting and pasting the unaligned sequences into the webservers provided by the EBI:

These are the alignments obtained using default settings

Note also that CLUSTAL and T-COFFEE have shortened considerably the length of the sequence name/description information provided for each sequence - this could cause trouble if you were comparing alignments in a way which depended on the sequences in different alignments have the same name/description information

Below are the alignments prepared in both .aln and .fasta format, with the order of the sequences determined by running initially MUSCLE over the sequences, and then fixing the order of the sequences when aligning using MAFFT, T_COFFEE, and CLUSTAL

MAFFT

MUSCLE

T-COFFEE

CLUSTAL

Note that by default, as provided above, CLUSTAL and T-COFFEE provide the alignment in .aln/CLUSTAL format, while MAFFT and MUSCLE have returned it in FASTA format.

We'll now look at three different ways of changing the order of sequences in the alignments:

Manually adjust sequence order using an MSA editor:

Disadvantage: time-consuming - in particular, your-time-consuming, rather than just taking a long time for the computer to calculate while you go off and do something else
Advantage: you can organise the sequences exactly the way you want/need to have them

Sort sequence order automatically using the MSA editor JalView

By percentage of pairwise sequence identity (%ID)

Disadvantage: %ID is is measured in the context of the current alignment - thus, if two alignments contatining the same sequences are aligned differently (which may well be the case if the alignments were built using different software) then the sorted order of the sequence might be different
Advantage: more similar sequences are grouped next to each other in the alignment - this makes patterns of conservation/divergence within the alignment easier to spot

Alphabetically (by the first word in the sequence description lines)

Disadvantage: only if the sequences are named in such a way that the identifiers of similar sequences are more similar to each other than they are to more divergent sequences will the sequences be organised so that more similar sequences are adjacent within the alignment (however, if they are named P53_HUMAN, P54_MOUSE, P63_HUMAN, P63_MOUSE or so then this will not be such a problem)
Advantage: as long as the sequence identifiers are the same, the order of the sequences should also be the same for all the alignments

(a) Pre-calculate a quick alignment (using MUSCLE) to order the sequences by %ID (b) remove all gaps from this initial alignment (c) calculate MSAs from this ungapped alignment using the settings that deliver the resulting alignments with the sequences in the same order as provided in the input sequence file

Advantages: (i) sequences are ordered by %ID, which provides the best overview of the conservation/divergence represented in the alignment (ii) ordering of the sequences is independent of any problems associated with the naming of the sequences (i.e. as long as the sequences remain the same, the identifiers can change and the order will remain the same)
Disadvantage: somewhat time-consuming

To specify whether the order of sequences in the output alignment should be that (i) of the order of sequences in the input file or (ii) the order generated by the alignment software (typically or perhaps even always this clusters together sequences in the alignment that are more similar to each other), look around at the options provided by the submissions pages for the MSA software webservers - from which you should choose either (i) input or (ii) aligned

MUSCLE - "OUTPUT ORDER"
MAFFT - "REORDER"
T-COFFEE - "ORDER"
CLUSTAL - "OUTPUT ORDER"

Exercise

Align this set of GTP-binding proteins using several of the different automatic MSA tools listed above.

Sort the order of the sequences in the alignments using each of the different approaches described above.

Which of the ways of sorting the sequences makes it easiest to compare the results of the analyses from the different MSA packages? Why to you prefer that method compared to the others?

Comparing Automatic MSA Tool Results to a Reference Alignment

Demonstration

As you'll do in the exercises below, in this demonstration we'll look at what happens when we re-align a reference alignment from BaliBase (1aboA_ref1) using several different pieces of automatic alignment software.

Exercises

You will be told to take one from the list of BaliBase2 list of reference alignments.

Download the RSF format of the alignments of interest.

Load alignment into CLUSTALX locally, remove all gaps from all sequences (Edit->Select all sequences, Edit->Remove all Gaps) and save in FASTA format.

Alignment->Output Format Options->Output Order change to "Input" (this will preserve the order of the sequences in any alignments done by CLUSTALX to be the same as the order in the input alignment - makes it easier to compare the results of a CLUSTALX alignment with the reference alignment)

Re-align the sequences locally using CLUSTALX (Alignment->Do Complete Alignment) - REMEMBER TO CHOOSE A DIFFERENT NAME FOR THE ALIGNMENT OUTFILE!

Compare the automatic alignment created by CLUSTALX with the reference BaliBase alignment - it will help if you re-arrange the (vertical) order of the sequences in the CLUSTALX alignment.

How many columns are mis-aligned in the CLUSTALX alignment i.e. where residues are aligned in the same column in the CLUSTALX alignment which are not aligned in the same column in the BaliBase alignment? NB - Only consider those columns in the BaliBase alignment which are underlined/in upper case - these are the ones which are assessed as being reliably aligned according to structural comparisons

Do the same for other alignment tools such as:

MAFFT
PROBCONS
MUSCLE
T-COFFEE

MSA analyses from Start to Finish

Decision-making is central to building a good/accurate/appropriate MSA. Here are some examples of the kind of issues that require these decisions:

where to get an initial set of sequences from?
which software to use to create an initial alignment of the sequences?
whether to exclude any of the initial sequences from the analysis?

if so, which ones?

etc.

Your reasons for building the alignment - i.e. the biological question(s) you are aiming to address using it - play a very important role in helping you reach decisions in response to questions such as these. Thus, there is not just one way to build an alignment comparing mammalian COX genes - the alignment you build will be different, depending on what you want to use the alignment for.

Thus, going from the beginning to the end of the process of building an MSA only makes sense in the context of a particular biological questions.

The demonstration/exercise shown below aims to do two things

provide you with a concrete example of building an MSA from scratch
indicate ways in which you might vary this recipe when your sequences/aims are different from those described here (which they almost certainly will be when you carry out similar exercises in your own research)

Given that the other focus of this course is phylogeny, we will take evolutionary questions as the basis for this demonstration and exercise.

We'll go through the analysis step by step "online" using the projector.

After each step, we'll give you a chance to carry out the same actions as I've carried out on your own machines.

However, if you haven't time to carry out these analyses yourself, we've tried to provide files at each stage that you can use to move on to the next step, even if you haven't been able to complete the previous one.

Demonstration - DPW Linear Motif Conservation in the Epsins - an MSA Analysis from Start to Finish

The epsins is a family of proteins involved in regulating receptor-mediated endocytosis, which are characterised by the the Epsin N-terminal homology domain PFAM:ENTH;PF01417, e.g. EPN1_HUMAN or EPN1_RAT.

Most epsins contain functional instances of DPW/AP2alpha linear motifs, which mediate interactions with the alpha-subunit of adaptor protein complex 2 (AP-2 - there are two genes coding for alpha subunits in humans, AP2A1 and AP2A2).

We will now investigate the conservation of these motifs within the epsins by building an MSA of the family.

1. BLAST search at the NCBI using default settings using EPN1_RAT (here as FASTA format file) as a query (which we use rather than the human sequence as the 3D structure of one of the DPW motifs in this sequence has been determined).

This gets quite a long list of sequences - to get a more manageable number (too many sequences is difficult to look at, and processing them [e.g. making an MSA] is slower with more sequences), we'll collect

only the high-scorring ones from the SwissProt database (by re-doing the search, this time against SwissProt)

here's those sequences in fasta format

high-scorring non-mammalian sequences from the initial (default) nr database

here's those sequences in fasta format

2. Build an MSA of the swissprot sequences using the EBI MUSCLE webserver

3. Load the alignment into JalView and look for instances of the DPW motif (regular-expression pattern 'DP[WF]') using JalView's "Find" feature

Here are the sequences aligned by MUSCLE

4. To fill out the picture of conservation a bit better, include also the high-scorring non-mammalian sequences by either

aligning sequences in profile 2 against the pre-aligned mammalian sequences in profile 1 using CLUSTALX
including all sequences together in one file and aligning them all to each other at the same time

Here is the MUSCLE alignment obtained from simply aligning all the combined sequences together at once

5. We'll examine the alignment again using JalView, and will also adjust the alignment manually to try and line up mis-aligned DPW motifs (if there are any...)

Exercise - PTS-1 Linear Motif Conservation in the Catalayses - an MSA Analysis from Start to Finish

Carry out a similar exercise yourself, this time looking at the conservation of PTS-1 peroxisiomal targeting linear motifs in catalases.

1. Search the NCBI nr database using BLAST with the human catalase (UniProt:P04040;CATA_HUMAN) as a query

Use databases of different sizes, and try filtering by E-value threshold and taxonomic group to obtain a set of sequences - at the end of this exercise, I've put together some sets of sequences you might have obtained in this way, but before using them, try adjusting the BLAST searches yourself to obtain a set of sequences which is

relatively small (perhaps only 40 or 50 sequences)
relatively diverse i.e. doesn't contain many sequences which are very similar/nearly identical to each other
with more focus on sequences from "model" organisms

The "Taxonomy Report" from the BLAST output page might be useful when planning your search strategy

2. Automatically align the sequences using the EBI MUSCLE webserver

MUSCLE alignment in FASTA format, with sequence names altered to make it easier to interpret the results

3. Load the sequences into JalView and search for matches to the PTS-1 motif

regular expression is: ' (.[SAPTC][KRH][LMFI]$)|([KRH][SAPTC][NTS][LMFI]$)'

4. Identify some of the taxonomic groups whose catalase sequences do not contain this motif - can you propose some reasons why this might be the case?

BLAST searches

The "Taxonomy Report" from an initial BLAST search (with default settings, which collects only 100 alignemnts) shows hits with high scores only in the animals/metazoans.

To see if there are related sequences in non-animals, I repeated the search using an initial query that filtered out all metazoan sequences, and collected a larger number than the default 100 alignments

This shows the protein to be present in many eukaryotes and prokaryotes (unusual) - and a search restricted to the archaea shows significant hits also in that domain of life - thus we can look to survey the presence/absence of this motif from all three domains of life.

It help to label the sequences so that we can tell which taxonomic group they're from - so I retrieve several different sets of sequences into separate files (one taxonomic group per file), as this makes it easier to re-name them according taxonomy

Archaea

nr database search restricted to "archaea [ORGN]"

All the catalase archaea sequences from the nr database in FASTA format
Archaea catalase sequences in FASTA format after I've done "search and replace" in my editor (">gi" for ">ARCHAEA") so that all the archaea sequences can easily be identified later on

Eubacteria

even a search against only SwissProt yields very many sequences, so I restrict it to just a few well-known organisms
"Bacillus subtilis" [ORGN] OR "deinococcus radiodurans" [ORGN] OR "Escherichia coli" [ORGN] OR "Streptomyces coelicolor" [ORGN]

Eukaryotes

again, too many with an nr database search, so just ran against SwissProt and picked by hand from model organisms

Demonstration - COX1/COX2 Phylogeny - an MSA Analysis from Start to Finish

We want to make inferences about the phylogeny of the prostaglandin G/H synthases (COX1/COX2) gene family in mammals. In particular, we would like to estimate when the duplication event occurred that yielded these two members of the family in humans. Perhaps we are interested in this question due to reading this paper (here's a static version of the Entrez PubMed Abstract for the paper in case the NCBI pages are down) where we read about there being two COX genes in vertebrates/mammals and reference to an "invertebrate precursor" - when did this duplication occur?

This is a typical evolutionary question - using what we see today - (e.g contemporary human COX1/COX2 protein sequences), we try to infer when (and if possible by what mechanism) the differences we see today have occurred.

Identifying Our Initial Sequence of Interest

To begin, we need to find the sequence of at least one of our proteins of interest

A good place to start is searching the UniProt/SwissProt database, which provides high-quality, partially manually-curated sequence and functional annotation data for proteins (as content has relatively high accuracy, and is relatively small, making it faster to search)

We can do this using the EB-eye search tool

We might begin by querying with "human cox1"

This hits many (100+) records in the protein sequence databases

Ideally we would have found just a single record corresponding to our protein of interest - however for several reasons we might hit many records (as is the case here)

Records that mention interactions with cox1 (either as cross-references to other databases, or in the "free text" of the records)
More than one protein with the name cox1

Rather than investigate all these records to find the one we are looking for, we modify the query e.g. "human cox1 prostaglandin"

This matches just one record - reading the annotation we find that, indeed, this is the gene we are looking for

Notice how we need to use our biological knowledge to aid our bioinformatic searching - generally, the more we know about the biology of our sequences, the better we can judge the results of such searches, and the more interesting questions we can pose using our bioinformatic tools

This task (identifying the sequence/record associated with a given gene/protein/biological entity) is one of the main applications of bioinformatic tools (the other is to predict function/structure using comparison to similar genes/proteins/entities with demonstrated function/structure)

This was a relatively easy example - however, if we wanted to find additional information about the cox1 sequence e.g. 3D structure, known protein-protein or protein-chemical interactions etc., we would have to carry out similar kinds of searches of different databases (or, if we're lucky, use cross-references from the UniProt record) All this can take quite a lot of time and effort.

REFLECT Exercise

A short-cut that works in many circumstances is to use the automatic entity-recognition and biological-database integration software REFLECT

you can use the tool by cutting-and-pasting a URL into the bar at the top of the REFLECT webpage
or you can install a Firefox plugin to give you a button on your browser that will REFLECT any page you happen to be viewing

This is what happens when we REFLECT our abstract of interest

Protein/gene names are identified and supplied with pop-up boxes containing information about them
Note the "disambiguation" feature when several proteins have the same name
We can go from these pop-ups to find the sequences for the genes/proteins
Can do this with any HTML page e.g.

Chemicals are highlighted in a different colour and have different pop-ups
Note that, using OnTheFly you can do the same thing with PDFs, Excel spreadsheets, and other document types - useful for those of you obtaining gene-lists from screens/microarray experiments etc.

Anyway - here are the UniProt records of the two human sequences COX1 and COX2

Note that these are NOT related to the mitochondrial protein COX1_HUMAN - these mitochondrial proteins were one of the reasons why our initial search of the UniProt database hit more than one record.

Identifying an Initial Set of Proteins Related to Our Sequence of Interest

We now need to identify a set of proteins that have the same structure/are related to our sequence of interest

We usually do this by beginning either with

a pre-calculated alignment from one of many different web resources/databases - the prefered solution as it can save a lot of time!
the results of a BLAST search

Note that we've already mentioned that different biological questions may well be best addressed by different MSAs. Thus the alignments/sets of sequences we obtain from such resources are unlikely to be the best possible alignments to address our own specific biological question. However, they often provide a good starting place for such an analysis - by carrying out modifications such as those described below, they can be adjusted to suit our specific question

including additional sequences in the alignment
removing sequence from the alignment
adjusting/editing regions to correct misalignments
excluding unneeded regions/columns of the alignment

Identifying pre-aligned sequences related to COX-1/PTGS1

For example, querying TreeFam with the UniProt accession number for COX-1/PTGS (P23219) takes us to the TreeFam record for this family (TF329675) - whose full alignment you can find here.

Identifying Related Sequences using BLAST

To build the MSA as far as possible "from scratch", we'll begin by looking at how we'd use a BLAST search for this.

We begin by choosing

to search against the UniProt/SwissProt (we start with a small database to make the BLAST search quicker)
switching the low-complexity filter on (makes it much less likely that we will obtain alignments with low E-values between unrelated sequences)
using other parameters in their default settings

If this search gives us a good starting point for our analysis, then no need to change the way we did the search.

If we hit many sequences, then downloading them all and building an MSA will take longer than necessary.

To change the search so that we hit fewer sequences we could:

query a smaller database
restrict the search to certain taxa

If we find very few/no sequences, then we may well not have enough sequences to build an MSA appropriate for addressing our question of interest.

To carry out a search to identify more sequences we could:

query a larger/more comprehensive database (perhaps located somewhere other than at the NCBI e.g. one of the JGI genome-specific databases)
adjust the query parameters

change the protein substitution matrix

short sequence? go to a higher BLOSUM/lower PAM matrix
long sequence? go to a lower BLOSUM/higher PAM matrix

try different gap penalties

use a more sensitive search method, typically based on protein profiles

protein domain database e.g. PFAM or SMART
create your own protein profile

PSI-BLAST
HMMER

We examine the results of the BLAST search

there are several sequences in the database that are very similar to our query
there are also many that are likely related (based on annotation and E-values) but relatively dissimilar to the query

We will have around 40 sequences if we collect both of these two groups of sequences - this is a reasonable number to manipulate, so we will take all of them.

We can easily get both GenPept and FASTA format records (and save them locally) for these sequences by re-formating the BLAST output to show only those with low-enough E-values

When I carried this out recently, a threshold of 1e-5 was appropriate
We can get an overview of the different taxa represented in the list by opening the "[Taxonomy reports]" link
here is a link to the sequences in FASTA format

Building an Initial MSA

We build an initial MSA from these sequences using MUSCLE at the EBI

MUSCLE is one of the faster MSA tools available, and is reasonably accurate - thus it is a good choice for carrying out an initial alignment
Here is an example of the output from MUSCLE

We examine the alignment in CLUSTALX

There are clearly two sets of sequences, within which all members of a set are more similar to each other than they are to members of the other set.
Both of our sequences of interest are in the same set - all of which are from mammals (apart from one chicken sequence).

Collecting Additional Sequences of Interest

Perhaps the duplication occurred within the mammals - but we'd be interested to supplement these sequences with some from other chordates, so we'll query a larger database, restricting our search to non-mammalian chordates.

To identify the terms used by the NCBI taxonomy to describe these taxonomic groups we query the taxonomy with the name(s) of organisms that we know belong to groups whose terms we want to use in restricting the BLAST query

here is the page for humans - this gives us the term used in the NCBI Taxonomy for mammals (Mammalia) and chordates (Chordata)

We repeat the BLAST search this time

using a larger sequence database ("Non-redundant protein sequences - nr")
restricting the sequence space to non-mammalian chordates using the "Entrez Query" box "Chordata [ORGN] NOT Mammalia [ORGN]"
here is the result page
we check the taxonomic groups represented in the set of hit sequences - and see a reasonable representation of different groups
to make the next steps quicker, we'll restrict the set of sequences we'll collect using an E-value threshold of 1e-20
if doing this analysis "for real" I would probably have

collected more sequences at this stage i.e. used a higher E-value threshold
extended the range of organisms searched against - perhaps querying against non-mammalian Metazoans as a first step, restricting this more and more if too many sequences were being identified

here are these sequences in FASTA format

Combing Sequence Sets in one Alignment

It's probably easiest just to bundle all the sequences we have collected so far into one file and to submit this to MUSCLE again

combined FASTA file
To be kind to the EBI webservers, rather than all submitting large-ish (84 sequences is not a "small" query) queries to them at the same time, we will use MUSCLE from the command-line

muscle -in combined_Sp_nrNonMammalChordates.fasta -out combined_Sp_nrNonMammalChordates.fasta.muscle

MUSCLE alignment

In this case, we did not spend any time/effort creating a high-quality initial alignment e.g. by correcting mis-alignments, excluding unneeded sequences etc.

Sometimes, however, we want to incorporate additional sequences into an alignment that we have previously edited and manipulated - in which case we will not want to simply combine all (old and new) sequences together and completely re-align them.

Instead, we want to align the new sequences to the initial alignment while preserving the relationships of the sequences in the initial alignment to each other.

We can do this using CLUSTALX (in "profile alignment" mode.)

We can also do this using t_coffee (but only from the command-line)

t_coffee sequences.fasta -profile=prf1.aln -outfile=combined_profiles.aln

sequences.fasta contains any number of sequences - these are assumed to be unaligned
prf1.aln is a sequence alignment
the sequences in sequences.fasta will each be aligned, individually, to prf1.aln, to create a new MSA incorporating all the sequences in sequences.fasta and prf1.aln
it's possible to integrate more than one profile at a time

t_coffee sequences.fasta -profile=prf1.aln,prf2.aln,prf3.aln -outfile=combined_profiles.aln

MUSCLE also allows you to align two alignments with each other (known as profile-profile alignment), although not a set of sequences individually to a pre-existing alignment - this functionality is also only available from the command line

muscle -profile -aln1.fasta -in2 aln2..fasta -out combinedAlignment.fasta
this aligns the alignments in the two files aln1.fasta and aln2.fasta to each other and writes them to combinedAlignment.fasta

Adjusting Sequence Identifiers

As indicated, to make decisions about which sequences to exclude, we need to always consider the biological questions we are trying to address with our alignment.

In this case, we are interested in the timing of the duplication event leading to the two mammalian COX genes - thus we need to change the sequence identifiers so that we can:

use them to recognise the organisms they come from.
as we have multiple sequences from the same organism in our set, we also want to be able to quickly determine which of several sequences from a given organism is associated with a given identifier

When choosing new identifiers for the sequences we need to take into account that:

each sequence in the file should have a unique identifier - many/most MSA software will fail if this is not the case.
some MSA formats only allow a maximum of 10-character identifiers (i.e. PHYLIP format) - therefore we need to make sure that the first 10 characters for every sequence identifier are unique. Thus, we typically place a unique identifier (e.g. a UniProt Accession Number) in the first characters of the identifier
quite a lot of the software that uses MSAs is unable to accept sequence identifiers that include anything other than the alphanumeric characters and underscore (i.e. a-z, A-Z, 0-9 and _) - some software even has problems with lowercase letters

Here is a copy of the MUSCLE alignment created above with the identifiers altered - although note that the policy described above has not been carefully followed closely.

Preparing a "Final" Alignment

Almost always you will want to exclude some of the sequences from such an alignment. This might be for several different reasons:

Many identical/almost identical sequences in the alignment

Keeping them in the alignment adds very little extra information to your analysis, and can slow calculations down considerably
You would usually remove all but one of each such group of sequences
Note, however, that for some applications you may expect and even want many very similar sequences - so this (like everything on this list) needs to be interpreted in the context of the aim of the analysis

Sequences look like they may contain errors

Positions in a sequence where the residue described in the sequence file is different from that which is specified in the organism's genome will introduce systematic errors into any analysis of that sequence that makes use of that residue. Clearly we want to avoid such errors in our analyses

Such errors could have many different sources - for example:

sequencing errors
contig-assembly errors
gene-prediction errors

Likely errors are recognised by finding regions of an alignment where

all/most sequences are conserved
where the likely "wrong" sequence isn't conserved
where the rest of the "wrong" sequence (especially the region flanking the likely error) is conserved in a similar way to other sequences in the alignment
THIS NEED NOT BE AN ERROR! Instead it could be the signature of some atypical evolution occurring in that region of the sequence. However, it's more likely to be an error

Unless you are particularly interested in that sequence, you should discard it
If you ARE particularly interested in that sequence, you will need to do further investigations to try and determine whether or not it is an error

Sequences that contain many gaps

Many analyses (particularly evolutionary analyses) can only make use of columns in which residues are present from all the sequences included in an alignment
For such analyses, as far as possible, we want to exclude sequences that include gaps in columns that we would otherwise expect to use in our analysis
Again, these gaps might be due to errors, or be biologically valid/relevant - as before, the decision on whether to keep a particular sequence or not depends on the aim of the alignment - if the sequence is particularly important for the analysis, we will of course keep it in the alignment, even if it is present in a short fragment

Here is a copy of the above MUSCLE alignment, with altered sequence identifiers, trimmed of many sequences due to reasons described above.

Estimating a Phylogeny from the Alignment

Before estimating the phylogeny, we need to take into account that phylogenetic applications of MSAs are known to be strongly influenced by the inclusion of columns in the analysis where residues in the column are not related by simple substitution events.

Therefore, we use the GBLOCKS server to extract those regions of the alignment that are considered most likely to contain only columns where all residues are related via substitution events.

Here is the HTML page showing which regions GBLOCKS kept and discarded from the alignment
Here is the FASTA format GBLOCKed alignment

Finally, we use a very quick and innaccurate method to estimate phylogeny from this alignment - neighbour-joining as implemented by CLUSTALX - click here for a tree estimated in this way

(We will look at a more accurate approach in tomorrow's session which focuses on phylogenies)

Below is an image of the tree estimated in this way using NJplot

Interestingly, there seem to be two copies of the gene also in the urochodate Ciona intestalis and the cephalochordate Branchiostoma floridae - however, the duplications appear to have occurred independently of each other, and of the duplication that occurred to yield the two vertebrate copies of the gene.

Exercise - Mammalian CDKs 1/2/3 Phylogeny - an MSA Analysis from Start to Finish

Now you have a chance to carry out an similar analysis, this time on your own.

This time, suppose you are interested in a similar question, only this time you want to investigate the timing of the duplications that yielded the three mammalian CDK genes CDK1, CDK2, and CDK3 (for example, see this abstract [Santamaria et al.2007] which mentions all three of the proteins).

Identify the UniProt record(s) associated with at least one of these genes
Use BLAST to identify an appropriate set of proteins to build an alignment from, and download these sequences to a local FASTA-format file
Use MUSCLE to build a multiple sequence alignment of the sequences

if possible, use MUSCLE from the command line e.g.

muscle -in unalignedSequenceFile.fasta -out alignedSequenceFile.fasta

if not, use the MUSCLE server at the EBI

Edit the names of the sequences so that you can tell which organisms each sequence comes from (while keeping the sequence identifiers unique)
Examine the alignment in CLUSTALX (or Jalview)

remove sequences from the analysis to make it more appropriate for the phylogenetic analysis
if necessary, carry out further BLAST searches to obtain additional sequences for the alignment and repeat the MSA and sequence-removal stages

Remove poorly-aligned regions from the alignment using the GBLOCKS server
Quickly estimate a phylogeny from the GBLOCKed alignment, and visualise it using NJplot

Based on this tree, when do you think the three lineages diverged from each other? (i.e. when do you think their ancestral genes underwent duplications?)

If you have time, repeat a similar analysis as described above, this time beginning with the appropriate TreeFam family as the initial source of sequences
(If you can't find the appropriate TreeFam record, here's the link to it).

Demonstration - Phylogeny From Start to Finish with Phylogeny.fr

Another example of an application of MSAs in phylogeny estimation. With this demonstration, the focus is on:

preparing the MSA, not estimating the phylogeny
on keeping comments/asides to a minimum

The Question

We'll try and identify the lamprey sequence (lampreys in wikipedia, UCMP, Tree of Life web project) most closely related to human UniProt link (UniProt:O00622;CYR61_HUMAN, here in FASTA format).

Acquiring initial set of sequences

We will begin with either a set of pre-aligned sequences from a database or unaligned related sequences using BLAST

Pre-aligned sequences from a database

We'll follow the link from the UniProt record for human CYR61 to the HOGENOM database to get our hands on an initial alignment

Unaligned sequences using BLAST

We'll search with the FASTA sequence of CYR61 (here is a local version of the file in case there are problems with UniProt) at the NCBI BLAST webserver

Initial search under default conditions

Output page
Set of sequences obtained using default parameters (unaligned) as of 6.5.2011 (nr database)

too many sequences, not enough taxonomic range, so we repeat using a smaller database, SwissProt

Search against SwissProt with default conditions

Output page
Set of sequences with E-value less than 0.0001
not a good enough representation of taxa of interest - repeat the initial search (against nr database) to have sufficient diversity, but restricting matches to a selected subset of organisms i.e. humans, african clawed frog (Xenopus laevis), zebrafish (Danio rerio), and non-bony-fish (NOT Euteleostomi) deuterostomes

Search against nr, showing 250 matching sequences, with following Entrez query restriction of search space

"homo sapiens" [ORGN] OR "xenopus laevis" [ORGN] OR "Danio rerio" [ORGN] OR "Scyliorhinus torazame" [ORGN] OR (Deuterostomia [ORGN] NOT Euteleostomi [ORGN])
Output page
Set of 52 sequences with E-value less than 1e-51 (chosen so that there aren't too many sequences to align automatically to make exercise go quicker - usually I'd begin with a larger set of sequences

This is a small enough set to begin building an MSA from.

Notice how knowledge of the biology of the system (in this case, the taxonomy/phylogeny of my organisms of interest) drives our decision making during the analysis.

Building an automatic MSA

Using the MUSCLE webserver at the EBI

set of aligned sequences

Refine the MSA

Working with the initial MSA, we usually want to remove some of the sequences (as their unnecessary or would even be problematic for the downstream analysis), and add sequences that are not currently in the alignment.

We might also want to correct possible errors in the alignment, and/or realign the sequences once some have been removed.

Remove unnecessary sequences

To generally reduce the redundancy of the dataset i.e. to avoid having too many similar sequences, we can use an automatic tool:

Decrease redundancy at ExPASy

Produces this alignment in FASTA format that contains fewer sequences than in the input file
Has the problem that we loose the names of the organisms associated with the sequences!

Using JalView (Edit->Reduce Redundancy)

Gives this set of sequences, preserving information on organism names!
Has the problem that we don't have (as far as I can tell!) control over which sequences in a redundant set are kept/removed
However, I still prefer this solution - I just keep an eye on what's removed and if something is removed that I need, I just add it back later

CD-HIT also can be used to reduce redundancy; this tool also preserves the names of the organisms in the output reduced-redundancy alignment

output alignment from CD-HIT

Even if we use such tools as a start to the analysis, we typically also remove additional sequences "by hand"

This alignment has had several additional "unusual" sequences removed

Add missing sequences

There are (well, there were on 5.5.2011) no lamprey sequences in our alignment at the moment!

Obviously we need some to address our biological question of interest.

This link describes how we could go about acquring some - however, to save time, rather than describing how we get them, here is a direct link to an appropriate set of lamprey sequences

We can build an MSA that includes these additional sequences either by:

Realigning all sequences (initial sequences plus additional sequences) from scratch

For example using the MUSCLE webserver at the EBI

Aligning the additional sequences to the initial alignment.

There are not many web tools that allow us to do this, although many automatic MSA packages provide this functionality in the command-line version of their programs e.g. T_COFFEE, MUSCLE, and CLUSTALX.
A local installation of CLUSTALX allows us to do this via a graphical user interface

Reduce the MSA

When using an MSA to estimate phylogeny and/or analyse evolutionary substitution processes, we usually need to remove columns from some of the columns from the alignment before estimating these evolutoinary paramters from our data.
We can do this:

automatically, using GBLOCKS
"manually" using an alignment editor, such as JalView

We have to rename the sequences so that they could be accepted by GBLOCKS to give this file here

Here is the result file from GBLOCKS

Estimating the phylogeny

We can do this using:

CLUSTALX (which yields a rather inaccurate estimate of the phylogeny)
Phylogeny.fr (which should give a somewhat more accurate estimate of the phylogeny)
or do it properly... (and spend a long time doing it)

Examining the phylogeny

We can examine the resulting tree files using software such as

phylowidget (a web tool)
NJplot
Dendroscope

Exercise

Carry out a similar analysis, based on the sequences given here but addressing the question which human sequences in the shark sequences AAL37727 most closely related to.

Patterns of Conservation in Protein Alignments

Either:

Use BLAST and automatic sequence alignment software such as MUSCLE at the EBI, to build MSAs
Obtain sets of pre-aligned sequences for which there are structures, such as HOMSTRAD

Examine the protein structures (which you might obtain from the PDBe or the PDB) using PyMOL

Identify the regions of protein sequence within the alignment that correspond to secondary structure elements (helices, strands, turns) within the protein structure.

The aim is to build/obtain an alignment which highlights patterns of conservation/variation in different kinds of structural element.

You might need to collect and incorporate additional sequences in your alignment and/or remove some sequences from the alignment to make these patterns as clear as possible.

If possible, change the colouring of the PyMOL structure, and save a session file and image of the structure that highlights these patterns.

Example starting points for your analysis:

While this application is perhaps of limited use in your own research, it should serve to highlight, again, the extent to which an awareness of the specific application for which you are building an MSA governs the decisions you make when preparing the alignment.

Conserved Hydrophobic Core

Just to remind you - here we can see clearly that the non-polar residues can be almost completely excluded from the core of a globular protein.

Visualising and Editing Alignments

The aim of this exercise is to give you some practice using alignment viewers/editors to (i) identify and (ii) correct mis-aligned regions.

Download the alignments in the left column of the table below, all of which contain regions that are misaligned (in comparison with a set of reference alignments from BaliBase).

The second column links to an image of the alignment, with a box surrounding a region that contains a misalignment.

Using Jalview and CLUSTALX, try to identify and correct the misalignment.

To see whether you were successful, compare your edited alignments with the reference alignments (provided in the third column in FASTA format, in the fourth as HTML, which highlights the regions BaliBase is confident are correctly aligned). The fifth column links to an image highlighting more closely the region that is misaligned, along with an indication of how the region might be corrected. The final column links to the BaliBase page for the alignment.

Alignment with Misalignment	Misaligned Region Highlighted	Reference BaliBase Alignment (FASTA)	Reference BaliBase Alignment (HTML - Core Regions Highlighted)	Edit Highlighted	BaliBase Page
BB12003.mis.fasta	BB12003.mis.jpg	BB12003.ref.fasta	BB12003.html	BB12003.edit.jpg	BB12003-BaliBase
BB12032.mis.fasta	BB12032.mis.jpg	BB12032.ref.fasta	BB12032.html	BB12032.edit.jpg	BB12032-BaliBase
BB50013.mis.fasta	BB50013.mis.jpg	BB50013.ref.fasta	BB50013.html	BB50013.edit.jpg	BB50013-BaliBase

Back To Common Course Content page