text: user3585: Background
  I have a fasta file with protein sequences. The first sequence in the fasta file is the reference sequence and every subsequent sequence is something I want to align to the reference. The end goal is to identify every unique sequence and be able to compare the regions of variation.
  What I have tried
  I'm looking into different tools such as Clustal Omega and MAFFT, but there is no option to select the reference sequence against which all the other protein sequences will be aligned.
  Questions
  Q1: I want perform one pairwise global alignment of each sequence in the file against the first (reference) sequence. What software can I use to do this?
  Q2: What specific commands can I use to achieve Q1? I'm looking for specific instructions on how to complete this, not just a suggestion of a tool, since I'm not able to find anything from a relatively comprehensive search of available tools.
  
  user2085: Q1
  From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?
  If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.
  Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like bwa mem or minimap2 and to use a bam file.
  Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.
  Q2 - MAFFT/Python
  In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.
  Here is the fasta file from Wikipedia. I called it seqs.fasta:
  >reference
  MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
  EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
  VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
  ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*
  
  >MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
  ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
  FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
  DIDGDGQVNYEEFVQMMTAK*
  
  >gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
  LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
  EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
  LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
  GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
  IENY
  
  Here is a script, mafft_ref.py. Run it with python mafft_ref.py seqs.fasta:
  from Bio import SeqIO
  from Bio.Align.Applications import MafftCommandline
  import tempfile
  
  fasta = sys.argv[1]
  
  reference = 0
  count = 0
  with tempfile.NamedTemporaryFile() as temp:
      for record in SeqIO.parse(fasta, "fasta"):
          if count == 0:
              reference = record
          else:
              SeqIO.write([reference, record], temp.name, "fasta")
              mafft_cline = MafftCommandline(input=temp.name)
              stdout,stderr=mafft_cline()
              print(stdout)
          count += 1
  
  Here is the output of python mafft_ref.py seqs.fasta printed to the terminal. Use python mafft.py prots.pep > output.txt to save the output to output.txt.
  >reference
  ------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
  VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
  ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
  SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
  CLRRDSSKIDTYLKLLNCRIIYNNNC
  >MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
  ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
  AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
  A---------------------------------------ELRHVMTN-----------L
  GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
  -----TAK------------------
  
  >reference
  ------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
  --LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
  GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
  RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
  IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
  PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
  >gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
  LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
  TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
  -FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
  SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
  LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
  ASI----------LYFSIILAFLPIAGXIENY--------------
  
  References:
  Some credit of mafft_ref.py to a biostars answer found here.
  
original_text: 3585: <h3>Background</h3>
  
  <p>I have a fasta file with protein sequences. The first sequence in the fasta file is the reference sequence and every subsequent sequence is something I want to align to the reference. The end goal is to identify every unique sequence and be able to compare the regions of variation.</p>
  
  <h3>What I have tried</h3>
  
  <p>I'm looking into different tools such as Clustal Omega and MAFFT, but there is no option to select the reference sequence against which all the other protein sequences will be aligned.</p>
  
  <h3>Questions</h3>
  
  <p><strong>Q1:</strong> I want perform one pairwise global alignment of each sequence in the file against the first (reference) sequence. <strong>What software can I use to do this?</strong></p>
  
  <p><strong>Q2:</strong> <strong>What specific commands can I use to achieve Q1?</strong> I'm looking for specific instructions on how to complete this, not just a suggestion of a tool, since I'm not able to find anything from a relatively comprehensive search of available tools.</p>
   2085: <h3>Q1</h3>
  
  <p>From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?</p>
  
  <p>If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.</p>
  
  <p>Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like <code>bwa mem</code> or <code>minimap2</code> and to use a bam file.</p>
  
  <p>Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.</p>
  
  <h3>Q2 - MAFFT/Python</h3>
  
  <p>In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.</p>
  
  <p>Here is <a href="https://en.wikipedia.org/wiki/FASTA_format" rel="nofollow noreferrer">the fasta file from Wikipedia</a>. I called it <code>seqs.fasta</code>:</p>
  
  <pre><code>&gt;reference
  MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
  EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
  VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
  ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*
  
  &gt;MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
  ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
  FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
  DIDGDGQVNYEEFVQMMTAK*
  
  &gt;gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
  LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
  EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
  LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
  GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
  IENY
  </code></pre>
  
  <p>Here is a script, <code>mafft_ref.py</code>. Run it with <code>python mafft_ref.py seqs.fasta</code>:</p>
  
  <pre class="lang-py prettyprint-override"><code>from Bio import SeqIO
  from Bio.Align.Applications import MafftCommandline
  import tempfile
  
  fasta = sys.argv[1]
  
  reference = 0
  count = 0
  with tempfile.NamedTemporaryFile() as temp:
      for record in SeqIO.parse(fasta, "fasta"):
          if count == 0:
              reference = record
          else:
              SeqIO.write([reference, record], temp.name, "fasta")
              mafft_cline = MafftCommandline(input=temp.name)
              stdout,stderr=mafft_cline()
              print(stdout)
          count += 1
  </code></pre>
  
  <p>Here is the output of <code>python mafft_ref.py seqs.fasta</code> printed to the terminal. Use <code>python mafft.py prots.pep &gt; output.txt</code> to save the output to <code>output.txt</code>.</p>
  
  <pre><code>&gt;reference
  ------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
  VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
  ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
  SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
  CLRRDSSKIDTYLKLLNCRIIYNNNC
  &gt;MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
  ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
  AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
  A---------------------------------------ELRHVMTN-----------L
  GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
  -----TAK------------------
  
  &gt;reference
  ------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
  --LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
  GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
  RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
  IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
  PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
  &gt;gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
  LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
  TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
  -FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
  SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
  LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
  ASI----------LYFSIILAFLPIAGXIENY--------------
  </code></pre>
  
  <p>References:
  <a href="https://www.biostars.org/p/343605/" rel="nofollow noreferrer">Some credit of <code>mafft_ref.py</code> to a biostars answer found here.</a></p>
  