Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-09-01T14:43:58.243Z
Answers.Text: 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.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-09-02T21:43:56.493Z
Answers.Text: I loved your approach it put a smile on my face. Thank you @albertr
Answers.UserDisplayName: Nilay Lale
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-09-01T14:20:17.733Z
Answers.LastActivityDate: 2022-09-01T14:20:17.733Z
Answers.UserDisplayName: albertr
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-27T11:47:05.220Z
Answers.LastActivityDate: 2022-08-27T13:49:05.303Z
Answers.LastEditDate: 2022-08-27T13:49:05.303Z
Answers.UserDisplayName: M__
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-26T18:28:02.677Z
Comments.Text: 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.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-27T07:39:33.360Z
Comments.Text: 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?
Comments.UserDisplayName: Chris_Rands
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-08-30T00:24:33.957Z
Comments.Text: 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.
Comments.UserDisplayName: Nilay Lale
Tags: <proteins><multiple-sequence-alignment>
Title: Select representetive sequence within MSA (protein)
ThreadText: # 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.
  
  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
  
  ```
  
  > **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 @username_3  
  
  ---
  
  ### 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.
  
  
  ---
  
  