How to extract a subset of protein sequences from a big fasta file with python?

Question

I'd like to extract a subset of protein sequences from a .fasta file (swissprot_canonical-isoforms.fasta) in a new file (selected_proteins.fasta) if protein IDs are listed in a txt file (Interested proteins.txt).

Below shows part of protein sequences in the swissprot_canonical-isoforms.fasta. Protein IDs are shown between two "|" in the lines that start with ">". For example, "P04637" is a protein ID.

>sp|P04637|P53_HUMAN Cellular tumor antigen p53 OS=Homo sapiens GN=TP53 PE=1 SV=4
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP
DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK
SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE
RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS
SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP
PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG
GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD
>sp|P04637-2|P53_HUMAN Isoform 2 of Cellular tumor antigen p53 OS=Homo sapiens GN=TP53
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP
DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK
SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE
RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS
SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP
PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQDQTSFQKENC
>sp|P04637-3|P53_HUMAN Isoform 3 of Cellular tumor antigen p53 OS=Homo sapiens GN=TP53
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP
DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK
SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE
RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS
SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP
PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQMLLDLRWCYFLINSS

Here are some protein IDs listed in the Interested proteins.txt

Q6ZWH5
Q8NG66
P51955
P51957
P04629

The final output should look like this (only listed the sequence for Q6ZWH5 as an example) :

>sp|Q6ZWH5|NEK10_HUMAN Serine/threonine-protein kinase Nek10 OS=Homo sapiens GN=NEK10 PE=2 SV=3
MPDQDKKVKTTEKSTDKQQEITIRDYSDLKRLRCLLNVQSSKQQLPAINFDSAQNSMTKS
EPAIRAGGHRARGQWHESTEAVELENFSINYKNERNFSKHPQRKLFQEIFTALVKNRLIS
REWVNRAPSIHFLRVLICLRLLMRDPCYQEILHSLGGIENLAQYMEIVANEYLGYGEEQH
TVDKLVNMTYIFQKLAAVKDQREWVTTSGAHKTLVNLLGARDTNVLLGSLLALASLAESQ
ECREKISELNIVENLLMILHEYDLLSKRLTAELLRLLCAEPQVKEQVKLYEGIPVLLSLL
HSDHLKLLWSIVWILVQVCEDPETSVEIRIWGGIKQLLHILQGDRNFVSDHSSIGSLSSA
NAAGRIQQLHLSEDLSPREIQENTFSLQAACCAALTELVLNDTNAHQVVQENGVYTIAKL
ILPNKQKNAAKSNLLQCYAFRALRFLFSMERNRPLFKRLFPTDLFEIFIDIGHYVRDISA
YEELVSKLNLLVEDELKQIAENIESINQNKAPLKYIGNYAILDHLGSGAFGCVYKVRKHS
GQNLLAMKEVNLHNPAFGKDKKDRDSSVRNIVSELTIIKEQLYHPNIVRYYKTFLENDRL
YIVMELIEGAPLGEHFSSLKEKHHHFTEERLWKIFIQLCLALRYLHKEKRIVHRDLTPNN
IMLGDKDKVTVTDFGLAKQKQENSKLTSVVGTILYSCPEVLKSEPYGEKADVWAVGCILY
QMATLSPPFYSTNMLSLATKIVEAVYEPVPEGIYSEKVTDTISRCLTPDAEARPDIVEVS
SMISDVMMKYLDNLSTSQLSLEKKLERERRRTQRYFMEANRNTVTCHHELAVLSHETFEK
ASLSSSSSGAASLKSELSESADLPPEGFQASYGKDEDRACDEILSDDNFNLENAEKDTYS
EVDDELDISDNSSSSSSSPLKESTFNILKRSFSASGGERQSQTRDFTGGTGSRPRPALLP
LDLLLKVPPHMLRAHIKEIEAELVTGWQSHSLPAVILRNLKDHGPQMGTFLWQASAGIAV
SQRKVRQISDPIQQILIQLHKIIYITQLPPALHHNLKRRVIERFKKSLFSQQSNPCNLKS
EIKKLSQGSPEPIEPNFFTADYHLLHRSSGGNSLSPNDPTGLPTSIELEEGITYEQMQTV
IEEVLEESGYYNFTSNRYHSYPWGTKNHPTKR

Is there any way to do this with python? Any help would be highly appreciated.

Some resources that will help: 1. Use [Biopython](https://biopython.org/) 2. [See example #4 here](http://training.scicomp.jic.ac.uk/docs/python_for_biologists_book/parsing_fasta_files.html#exercises) 3. [Parsing Fasta files with Python](https://stackoverflow.com/a/7655072/1431750) 4. Other [SO questions involving Fasta and Python](https://stackoverflow.com/search?q=fasta+%5Bpython%5D) — aneroid, Jun 19 '20 at 08:01
@aneriod: do not forget [the official Python Absolute Beginners' Tutorial](https://docs.python.org/3/tutorial/) – you don't know OP's level of Python. — Jongware, Jun 19 '20 at 08:27

Michael Mitter · Accepted Answer · 2020-06-19T08:41:54.493

1

You can do this using pyfasta, an interface to the FASTA format from python.

from pyfasta import Fasta

f = Fasta('fasta.fa') # open the file

targets = {"P04637"} # define your target IDs

selection = {}

for key in f:
    candidateKey = key.split("|")[1]
    if candidateKey in targets:
        selection[key] = f[key][:]
        print(key)
        print(selection[key])

Output:

sp|P04637|P53_HUMAN Cellular tumor antigen p53 OS=Homo sapiens GN=TP53 PE=1 SV=4
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAP
TPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAM
AIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS
SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKK
KPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

edited Jun 19 '20 at 08:41

answered Jun 19 '20 at 07:50

Michael Mitter

545
3
9

Thanks for answering my question! It is very useful. But the output format is different to what I want. I don't want a dictionary. Is there a way to covert that to what I listed above? – Chemist learns to code Jun 19 '20 at 08:14
I edited my answer to print the output nicely in addition to storing it in a dictionary. If you want to save this information to a file, you can do that similarly. – Michael Mitter Jun 19 '20 at 08:42

How to extract a subset of protein sequences from a big fasta file with python?

1 Answers1