<p>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.</p>
<ol>
<li>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 <strong>real</strong> sequence will be your representative.</li>
</ol>
<pre><code>from Bio import AlignIO
from Bio import pairwise2
from collections import Counter

aln_file = &quot;MSA.aln&quot;
alignment = AlignIO.read(aln_file, &quot;clustal&quot;) # 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 = &quot;&quot;.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 &gt; highest_score:
        highest_score = pairwise_score
        ref_seq = record
    elif pairwise_score == highest_score:
        # If they are equal, you can prioritize other things
        pass
</code></pre>
<ol start="2">
<li>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 <code>pip install python-Levenshtein</code>.</li>
</ol>
<pre><code>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 &lt; best_distance:
        best_distance = distance 
        ref_seq = record
    elif distance == best_distance:
        # If they are equal, you can prioritize other things
        pass
</code></pre>

CC BY-SA 4.0
2022-09-01T14:43:58.243Z
This is a very good use of `pairwise2`. `Levenshtein distance` in any other walk of life is cool, wrt to genetic distance would have its critiques.
M__
CC BY-SA 4.0
2022-09-02T21:43:56.493Z
I loved your approach it put a smile on my face. Thank you @albertr
Nilay Lale
CC BY-SA 4.0
2022-09-01T14:20:17.733Z
2022-09-01T14:20:17.733Z
albertr
<p>Suggestions in increasing order of relevance:</p>
<ol>
<li><p>The only program that is close is <a href="https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-03956-5" rel="nofollow noreferrer">rootdigger</a>.
This finds the sequence most likely to be the root, as the name suggests. Available via <a href="https://github.com/computations/root_digger" rel="nofollow noreferrer">GitHub</a>, but not via conda.</p>
</li>
<li><p>You could just write your <code>PIM</code> approach as small program? It wouldn't be difficult.</p>
</li>
<li><p>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*).</p>
</li>
<li><p>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).</p>
</li>
</ol>
<p>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.</p>
<p>*, Point 3 and 4 require a phylogram rather than a dendrogram of course.</p>

CC BY-SA 4.0
2022-08-27T11:47:05.220Z
2022-08-27T13:49:05.303Z
2022-08-27T13:49:05.303Z
M__
<p>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:</p>
<ol>
<li>I calculate a Percentage Identity Matrix (PIM) from the MSA file</li>
<li>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.</li>
<li>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.</li>
</ol>
<p>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.</p>

CC BY-SA 4.0
2022-08-26T18:28:02.677Z
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Community
CC BY-SA 4.0
2022-08-27T07:39:33.360Z
Can you give more context on *why* you want a representative sequence and why a consensus sequence is not appropriate? What is your use case?
Chris_Rands
CC BY-SA 4.0
2022-08-30T00:24:33.957Z
I am sorry for late reply, I am in the process of relocating to a different country. We are trying to pick a reference sequence when a reference seq does not exist. Later we will be doing some kind of variant/diversity analysis and we will be needing to map the Highly Convserved Regions to this selected reference.
Nilay Lale
<proteins><multiple-sequence-alignment>
Select representetive sequence within MSA (protein)
# Select representetive sequence within MSA (protein)

(Asked by: username_0 on 2022-08-26)

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.

**Tags:** proteins, multiple-sequence-alignment
> **Comments:**
> **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
> **username_2**: Can you give more context on *why* you want a representative sequence and why a consensus sequence is not appropriate? What is your use case?  
> **username_0**: I am sorry for late reply, I am in the process of relocating to a different country. We are trying to pick a reference sequence when a reference seq does not exist. Later we will be doing some kind of variant/diversity analysis and we will be needing to map the Highly Convserved Regions to this selected reference.  

---

## 2 answers

### Answer by username_3 on 2022-09-01. Score: 2

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.

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


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

> **Comments:**
> **username_4**: This is a very good use of `pairwise2`. `Levenshtein distance` in any other walk of life is cool, wrt to genetic distance would have its critiques.  
> **username_0**: I loved your approach it put a smile on my face. Thank you @albertr  

---

### Answer by username_4 on 2022-08-27. Score: 0

Suggestions in increasing order of relevance:

1. The only program that is close is [rootdigger](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-03956-5).
   This finds the sequence most likely to be the root, as the name suggests. Available via [GitHub](https://github.com/computations/root_digger), 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.


---

