<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>

CC BY-SA 4.0
2018-12-20T01:24:53.943Z
2018-12-20T01:47:52.637Z
2018-12-20T01:47:52.637Z
conchoecia
<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>

CC BY-SA 4.0
2018-12-20T08:29:02.487Z
Welcome to the site leeym. How many sequences do you have to compare with your reference sequence?
llrs
<phylogenetics>
Multiple sequence pairwise global alignment against reference sequence
# Multiple sequence pairwise global alignment against reference sequence

(Asked by: username_0 on 2018-12-19)

### 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.

**Tags:** phylogenetics
> **Comments:**
> **username_1**: Welcome to the site leeym. How many sequences do you have to compare with your reference sequence?  

---

## 1 answers

### Answer by username_2 on 2018-12-20. Score: 2 | Accepted answer

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.


---

