3
$\begingroup$

I want to compute the optimal alignment of two amino acid sequences as per the following definition from a patent:

"The percentage of identity between two peptidic or nucleotidic sequences is a function of the number of amino acids or nucleotide residues that are identical in the two sequences when an alignment of these two sequences has been generated. Identical residues are defined as residues that are the same in the two sequences in a given position of the alignment. The percentage of sequence identity, as used herein, is calculated from the optimal alignment by taking the number of residues identical between two sequences dividing it by the total number of residues in the shortest sequence and multiplying by 100. The optimal alignment is the alignment in which the percentage of identity is the highest possible. Gaps may be introduced into one or both sequences in one or more positions of the alignment to obtain the optimal alignment. These gaps are then taken into account as non-identical residues for the calculation of the percentage of sequence identity."

The Needleman and Wunsch implementation at NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) mostly works but not exactly. Thanks to @David & @WYSIWIG from a related SE Thread that suggested this (Computing Percent Identity between DNA / Amino Acid Sequence)

I want to know if there is a way to fix the mismatch.

e.g. My test case is:

Seq1: ABDE
Seq2: AAAAAAAAAAABCDE

The NCBI implementation yields the following alignment which has only 3 identical residues:

enter image description here

But shouldn't an optimal alignment with 4 identical residues be possible like so:

Seq1: ----------AB-DE
Seq2: AAAAAAAAAAABCDE

Thoughts? Any way to tweak the implementation to give the result I want? Alternatively any other algorithm that can be coerced into getting this alignment? BLAST or a variant?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
3
$\begingroup$

Needleman-Wunsch does an end-to-end (global) alignment (BLAST uses Smith-Waterman). Needle from the EMBOSS toolkit performs Needleman-Wunsch alignment. It will report the highest scoring alignment. I am not sure which alignment it reports when there are two of them with equal scores (I don't think it is random).

I just tried your case: replaced B with W as the former does not denote any specific amino acid (it is ambiguous). It does give:

     1 ----------AW-DE      4
                 || ||
     1 AAAAAAAAAAAWCDE     15

Note that you can change this behaviour by altering the gap-open and gap-extend penalties. You can also change the end-gap penalties (gaps in the beginning or end of the alignment; not in the middle) in Needle.

In this case gap-open=10, gap-extend=0.5, end gap penalty=false and matrix=BLOSUM62

For performing a local alignment you can use Smith-Waterman. It just aligns the highest scoring region and does not do an end-to-end alignment. You can change gap-open and gap-extend penalties in Smith-Wateman and BLAST too, but these algorithms do not begin or end an alignment with a gap.

This is the output of Smith-Waterman with gap-open=1 and gap-extend=0.5 and BLOSUM62 as scoring matrix.

     1  AW-D  3
        || |
    11  AWCD 14

For more information see What is the difference between local and global sequence alignments?

Improve this answer
$\endgroup$
3
  • $\begingroup$ Dang. What luck that I came across the wrong codes online at this link: fao.org/docrep/004/y2775e/y2775e0e.htm It says B denotes asparagine or aspartic acid. It is the top link at Google even. And the FAO seemed like an authoritative site. My bad! $\endgroup$
    – curious_cat
    Jun 22, 2016 at 17:33
  • $\begingroup$ @curious_cat Oh yes I didn't realize. Seems Needle accepts ambiguous codes like B and Z but many algorithms probably do not (though I am not too sure about this). Anyway, is your doubt clear? $\endgroup$
    – WYSIWYG
    Jun 22, 2016 at 17:43
  • $\begingroup$ Thanks. Yes, I think it is clear now. I am getting close to the patent optimal alignments by setting all penalties to zero. Maybe I will also provide a custom matrix now to make it a 1-0 alignment instead of the scoring function matrix that is natural propensity based like BLOSUM. $\endgroup$
    – curious_cat
    Jun 22, 2016 at 18:36
1
$\begingroup$

That is a very bad test case.

The problem is that the sequences are too short and involve a long repetition. This means that the default gap penalties and the gap-length penalties are not applicable. They are designed to work with longer sequences, where the penalty of inserting a gap can be offset by an increase in matches. In any case you can get bad alignments at the ends of DNA sequences where the stretch of DNA that is aligned by introducing a gap is ‘off-screen’.

Even though the algorithm may be watertight, its implement in an alignment program makes certain assumptions where values for scores and penalties are involved. You have to know what these are and when they are applicable. There are even situations where you need to produce your own matrix to force alignment between pairs that you know from other information must align (active sites, regulatory motifs). This is quite valid because the program doesn't know about biology — you do.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ An addition (on the same lines as providing your own matrix): PSI-BLAST does a calculation of a position specific scoring matrix using the homolog hits. With each iteration the hits are supposed to get better. This is a sort of machine learning; not really like putting your own matrix but it does take some biology into account. $\endgroup$
    – WYSIWYG
    Jun 22, 2016 at 20:29
  • $\begingroup$ As you say, this is a bit different. It's standard in MSA (multiple sequence alignment) — or at least Clustal does the same thing. But the key point is to think about the algorithm and the assumptions — difficult for a beginner obviously. $\endgroup$
    – David
    Jun 22, 2016 at 21:31
1
$\begingroup$

Your guess is made without knowing protein substitution probabilities.

Optimality of alignment in BLAST is a metrics of the scoring function. The scoring function depends on the word size (length of the seed that initiates an alignment), rewards and penalties for matches and mismatches - gap cost and substitution matrix. In general, BLAST uses BLOSUM and PAM- these matrices are based on evolutionary data. More about scoring parameters you could read in the BLAST Help.

You could fine tune these parameters in "Algorithm parameters" in the very bottom of the analysis page. Knowing parameters you could decide does your wild guess alignment make sense or not.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .