In [6]:
from sys import path
from os.path import join as joinpath
from os.path import normpath
path.append(normpath(".."))
from pyprot.data.fasta import getSequencesFromFasta
from pyprot.align.align import Align, Aligned
from pyprot.align.score import PSSM
A profile is a way to represent multiple aligned sequences and their similarities or patterns.
A protein domain is part of a protein that has its own structure and function, independent of the rest of the protein. Despite their use being similar across the proteins they can be found in, the same domains do not necessarily have the same amino acid sequence. However, they do have enough common amino acids to be identified and located by the means of sequence alignment. Our goal in this project is to locate a given domain within a given sequence, should it exist one or more times within that sequence. One such domain is the WW domain, known to be used by multiple species and sometimes included multiple times in a single protein. In this project we will use sequences from the WW domain that belong to human proteins to test our implementations : these have been saved to the file to-be-aligned.fasta. These domain sequences and many others can be found on the SMART database.
Multiple Sequence Alignment (or MSA) allows us to align multiple sequences together, and therefore observe how well some amino acids are preserved at some positions. This is of great interest in the study of domains, since it reveals which positions are of greater importance when trying to align new sequences to known domains. There are of course other applications, as well as numerous methods that achieve MSA, however we won't go into them since we won't be implementing any of these in this project. Two online tools were used to align all the sequences from to-be-aligned.fasta together : MUSCLE and CLUSTAL Omega. The resulting sequences can be respectively found in files msaresults-MUSCLE.fasta and msaresults-CLUSTALO.fasta.
A Position Specific Scoring Matrix (or PSSM) is a kind of profile that uses a matrix, where each column contains information about the same columns of the aforementioned sequences, and each row matches one amino acid. In this manner, each cell points to a specific column of all observed sequences, as well as one amino acid : the value it contains is the frequency of this amino acid within these columns. PSSMs can be used to align sequences with, the result indicating whether that sequence contains a subsequence that's similar to all of the sequences represented by the PSSM.
One way to represent a profile such as a PSSM is a table containing each value from the matrix. A WebLogo is a more human-readable representation, where amino acids that appear in each column are represented sequencially and with sizes depending on their frequency within that column. Here is an example of a WebLogo, gathered from this site:
The eponym class PSSM contains the frequency matrix as well as metadata that will allow us to use it efficiently as a scoring system when aligning sequences with it (more on that later). Here is the class, as well as a dictionary of absolute amino acid frequencies within the UniProt database (also used to provide scores) :
Suppose we want to align a sequence to a PSSM: it will be the PSSM that determines the scores for the alignment -based on the frequencies gathered from previously aligned sequences-. Sure, PSSM has "Scoring Matrix" in it, but we haven't talked about scores yet (that is, unless you've read the code above). That's when you should think : "oh, but I know what a scoring matrix is, it gives out scores for each amino acid pair, that's what this is". Well, when we were aligning two sequences together, that was the case since we made the assumption that the positions of amino acids didn't matter. However the goal now is to align one sequence with a load of other aligned sequences -represented by a PSSM- that often don't have the same amino acid at the same position. The frequency and range of amino acids can vary for each column of the PSSM, therefore the column value is required (instead of the "other" amino acid) to provide a score.
So how is that score calculated you ask ? Well, exactly the same way we did for the Blosum matrices: with the log-odds ratio $\frac{q_{a,c}}{p_a}$, where $q_{a,c}$ is the evolutionary probability of amino acid $a$ being located at column $c$, and $p_a$ is the random probability of amino acid $a$ (which doesn't depend on the column). What differs from the previous project is the way these terms are calculated:
This is going to be easy: in order to align a sequence with a PSSM, we can re-use the algorithm written for the first project to compute the alignment. The code has been adapted to better manage and represent MSAs, but the alignment is basically a local+suboptimal sequence alignment, with no affine gap penalty, where scores as well are provided by the PSSM. Note that gap penalties can be made to depend on the column and therefore are provided by the PSSM class, however in this case we use a constant gap penalty throug the whole alignment.
Here is the HMM Logo (similar to a WebLogo) of the WW domain, provided by this site (note that columns with no "main" amino acid are not represented):
Once aligned with the two mentioned online tools, we can represent our human proteins-only WW domains with WebLogos directly from this site. Here are the results of the CLUSTAL Omega alignment:
and the MUSCLE alignemnt:
Despite the ranges being slightly different between the WW domain and our subset aligned "by hand", the similarities speak for themselves: the same amino acids are of great importance to the domain, only "smaller" amino acids differ. Conserved positions are mainly 3,6,18,19,54,57.
Protein D6C652 (Transcriptional coactivator YAP1-A) domains according to UniProt:
This tells us that protein contains two WW domains, starting at positions 141 and 199. Let's see if these match our results for its alignment against our PSSMs for the two MSAs of the WW domain (note that by default, the function multiAlign only returns one result per optimal or suboptimal lookup, and only goes to suboptimal depth of 3):
In [7]:
musclePath = normpath("resources/fasta/msaresults-MUSCLE.fasta")
clustalPath = normpath("resources/fasta/msaresults-OMEGA.fasta")
testPath = normpath("resources/fasta/test.fasta")
for aligned in (musclePath, clustalPath):
print("\n >>> Creating PSSM from file {} ...".format(aligned), end="")
pssm = PSSM("WW domain")
for seq in getSequencesFromFasta(aligned):
pssm.add(seq)
pssm.setGapPenalty(4)
print(" done\n\n")
al = Align(pssm)
toalign = [p for p in getSequencesFromFasta(testPath)][0] #First protein
for aligned in al.multiAlign(toalign):
print(aligned)
Protein P46935 (E3 ubiquitin-protein ligase NEDD4) domains according to UniProt:
Let's see how these correlate to our own results:
In [8]:
for aligned in (musclePath, clustalPath):
print("\n >>> Creating PSSM from file {} ...".format(aligned), end="")
pssm = PSSM("WW domain")
for seq in getSequencesFromFasta(aligned):
pssm.add(seq)
pssm.setGapPenalty(4)
print(" done\n\n")
al = Align(pssm)
toalign = [p for p in getSequencesFromFasta(testPath)][1] #Second protein
for aligned in al.multiAlign(toalign):
print(aligned)
The PSSM provides the required information to find local alignments within other sequences. By looking at the WebLogos, we can see that both online tools have aligned our sequences from the WW domain in a very similar fashion and resemble the HMM Logo of the WW domain, yet the MUSCLE alignment provides better results than the OMEGA alignemnt.
All of the WW domains from both aligned sequences have been found, however the algorithm can probably be tweaked to have better start and end conditions for the alignment, since it sometimes returns only part of the domain. We can also notice that, when a complete domain is found, its start and end index within the sequence differ very slighlty (one or two amino acids) with the UniProt alignment; this could be due to tiny differences in the PSSM values, to the arbitrary selection of results from the backtracking (only the first is displayed for each lookup), or even to the constant and arbitrary gap penalty (different gap penalties per column have been implemented but not used).
Considering that the two tested sequences come from different organisms than the ones used for creating the PSSM, the results are positive overall.