Select representetive sequence within MSA (protein)

Question

I am trying to pick the best representative sequence from an amino acid MSA. I do not want to use the consensus, as it generates an artificial one. Here is my approach:

1. I calculate a Percentage Identity Matrix (PIM) from the MSA file
2. I check the highest value in each column and rank the respective sequences. What I mean by this, for example, the values in the first column would be the Percentage Identity of the first seq compared the rest of the dataset. If I pick the highest value in there, that would give me the sequences that are highly identical to the first seq.
3. I do this for the every column and rank the sequences. The top ranked sequence would yield me the best representetive sequence within my MSA data set.

The problem with this approach is that its computationally expensive to calculate the PIM as its doing a character match, comparing each sequence's position with all of the dataset. Is there any easier way to achieve this? Im not sure, if we could achieve something similar by looking at the tree. I think we can find the lowest common ancestor by the tree but I do not think that would yield a representetive sequence.

Answer 1

I have two suggestions that I used in the past (both using Python/BioPython though). They are arguably the same, but obtained through different paths.

1. Calculate the (artificial) consensus sequence of your alignment. Then, compare each sequence in your MSA to the consensus in a pairwise alignment, keeping the score. The highest scoring real sequence will be your representative.

from Bio import AlignIO
from Bio import pairwise2
from collections import Counter

aln_file = "MSA.aln"
alignment = AlignIO.read(aln_file, "clustal") # whatever format you have

for i in range(0,alignment.get_alignment_length()):  # iterate position by position
    aa_position = []

    # Get all residues in that position
    for record in alignment:
        aa_position.append(record.seq[i])

     # Use the most common residue as representative
    consensus_list.append(Counter(aa).most_common()[0][0])
consensus_seq = "".join(consensus_list)

ref_seq = None
highest_score = 0
for record in alignment:
    # globalxx finds the best global alignment between both sequences.
    # The first x indicates that identical characters are given 1 point, 0 otherwise
    # Second x indicates no gap penalties
    pairwise_score = pairwise2.align.globalxx(record.seq,
                                              consensus_seq,
                                              score_only=True
                                              )
    if pairwise_score > highest_score:
        highest_score = pairwise_score
        ref_seq = record
    elif pairwise_score == highest_score:
        # If they are equal, you can prioritize other things
        pass

2. You can again obtain the consensus sequence, and then calculate the Levenshtein (edit) distance between each record and the consensus. The sequence with the lowest distance will be the most representative. The easiest way to calculate it is via installing the python-Levenshtein module with pip install python-Levenshtein.

from Levenshtein import distance as edit_distance

# First, calculate the consensus as described above

ref_seq = None
best_distance = 999999999999
for record in alignment:
    distance = edit_distance(record.seq, consensus_seq)
    if distance < best_distance:
        best_distance = distance 
        ref_seq = record
    elif distance == best_distance:
        # If they are equal, you can prioritize other things
        pass

Answer 2

Suggestions in increasing order of relevance:

1. The only program that is close is rootdigger. This finds the sequence most likely to be the root, as the name suggests. Available via GitHub, but not via conda.

2. You could just write your PIM approach as small program? It wouldn't be difficult.

3. Alternatively, you could just select the longest, shortest or median branch length and use that as a representative sequence. This wouldn't be computationally expensive at all (you'd need the correct root however and a rooted tree*).

4. The shortest branch length from the root is a reasonable solution. You can do it via a GUI (if its e.g. a Bayesian tree, there's one program that does this), write you own code, or else manually via optical analysis (easy to spot the shortest branch in a rooted tree).

The shortest branch-length is the closest sequence to the ancestral sequence for the in-groups (non-root). You can obtain a reconstruction of the ancestral sequence - it is computationally (very) expensive and isn't a real sequence.

*, Point 3 and 4 require a phylogram rather than a dendrogram of course.