3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Welcome to the site leeym. How many sequences do you have to compare with your reference sequence? $\endgroup$
    – llrs
    Commented Dec 20, 2018 at 8:29

1 Answer 1

2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.