How can I pull out specific protein fastas from one file using information from the protein header?
1
0
Entering edit mode
5 weeks ago
ges29 ▴ 10

I have a file containing all the protein fastas from my genome of interest (downloaded from NCBI):

I've been conducting some analysis and have a list of genes I'm interested in; however, I don't have the accession numbers, I have a list of partial gene IDs assigned by the people who did the annotation:

Is there a way I can search the protein fastas file for a list of proteins using information from the protein headers?

I've tried grep -f genes_of_interest.txt protein_fastas_file.faa > output.fa

in an attempt to at least pull out the headers so I can get the accession numbers but all it did was return every single protein fasta (ie. exactly the same file as protein_fastas_file.faa). I'm assuming this is because the actual sequence part doesn't actually exist on a new line or something?

Thanks in advance for any help anyone can give!

protein grep accession fasta • 408 views
0
Entering edit mode

Input files are taken from GenoMax post

input:

$cat test.fa >QBMSHHS Hypotherical protein METSCH_A00100 someting MPPVGGKKAKKGILERLNAGEIVIGDGGFVFALEKRGYVKAGPWTPEAAVEHPEAVRQLHREFLRAGSNV MQTFTFYASEDKLENRGNYVLEKISGQEVNEAACDIARQVADEGDALVAGGVSQTPSYLSCKSETEVKKV >QBMSMMS Hypotherical protein METSCH_A00101 someting FLQQLEVFMKKNVDFLIAEYFEHVEEAVWAVETLIASGKPVAATMCIGPEGDLHGVPPGECAVRLVKAGA SIIGVNCHFDPTISLKTVKLMKEGLEAARLKAHLMSQPLAYHTPDCNKQGFIDLPEFPFGLEPRVATRWD >QQDSHHS Hypotherical protein METSCH_A20100 someting IQKYAREAYNLGVRYIGGCCGFEPYHIRAIAEELAPERGFLPPASEKHGSWGSGLDMHTKPWVRARARKE YWENLRIASGRPYNPSMSKPDGWGVTKGTAELMQQKEATTEQQLKELFEKQKFKSQ$ cat ids.txt

A00100


ouput:

$seqkit -w 0 grep -nr -f ids.txt test.fa >QBMSHHS Hypotherical protein METSCH_A00100 someting MPPVGGKKAKKGILERLNAGEIVIGDGGFVFALEKRGYVKAGPWTPEAAVEHPEAVRQLHREFLRAGSNVMQTFTFYASEDKLENRGNYVLEKISGQEVNEAACDIARQVADEGDALVAGGVSQTPSYLSCKSETEVKKV  Download seqkit from here: https://bioinf.shenwei.me/seqkit/download/ ADD REPLY 2 Entering edit mode 5 weeks ago GenoMax 102k $ more prot.fa
>QBMSHHS Hypotherical protein METSCH_A00100 someting
MPPVGGKKAKKGILERLNAGEIVIGDGGFVFALEKRGYVKAGPWTPEAAVEHPEAVRQLHREFLRAGSNV
>QBMSMMS Hypotherical protein METSCH_A00101 someting
FLQQLEVFMKKNVDFLIAEYFEHVEEAVWAVETLIASGKPVAATMCIGPEGDLHGVPPGECAVRLVKAGA
SIIGVNCHFDPTISLKTVKLMKEGLEAARLKAHLMSQPLAYHTPDCNKQGFIDLPEFPFGLEPRVATRWD
>QQDSHHS Hypotherical protein METSCH_A20100 someting
IQKYAREAYNLGVRYIGGCCGFEPYHIRAIAEELAPERGFLPPASEKHGSWGSGLDMHTKPWVRARARKE
YWENLRIASGRPYNPSMSKPDGWGVTKGTAELMQQKEATTEQQLKELFEKQKFKSQ

$more id A00100 A20100 # linearize fasta (code from @Pierre)$ awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' prot.fa | grep -f id | tr "\t" "\n" > new.fa

$more new.fa >QBMSHHS Hypotherical protein METSCH_A00100 someting MPPVGGKKAKKGILERLNAGEIVIGDGGFVFALEKRGYVKAGPWTPEAAVEHPEAVRQLHREFLRAGSNVMQTFTFYASEDKLENRGNYVLEKISGQEVNEAACDIARQVADEGDALVAGGVSQTPSYLSCKSETEVKKV >QQDSHHS Hypotherical protein METSCH_A20100 someting IQKYAREAYNLGVRYIGGCCGFEPYHIRAIAEELAPERGFLPPASEKHGSWGSGLDMHTKPWVRARARKEYWENLRIASGRPYNPSMSKPDGWGVTKGTAELMQQKEATTEQQLKELFEKQKFKSQ  Add fold -w 70 if you want the sequence lines folded every 70 aa. $ awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' prot.fa | grep -f id | tr "\t" "\n"| fold -w 70 > new.fa

0
Entering edit mode

I tried awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' prot_one.fa | grep -f id | tr "\t" "\n" > new.fa

but this has resulted in the same issue as before, the output file (new.fa) has all genes in it.

I checked this with grep -c "^>" new.fa - which returns 5800 genes.

but my list of genes I'm interested in only has 1550

0
Entering edit mode

I showed a working example modeled after what your screen shot had (BTW you should copy paste data rather than screen-shots). It should work. Not sure what is going on.

0
Entering edit mode

Apologies, I did try to do that with my original question but whenever I copied in the data the ">" symbol for the next gene causes half of it to turn into a quote rather than code. As below:

 >QBM85393.1 hypothetical protein METSCH_A00100 [Metschnikowia aff. pulcherrima]
MIQAHESLQFHQGFTQQMVKYLHKIVVEDFENFPGFHVLPKMHKSPPSSRPIVPTHSWGTMRISKVLSFMLRPQLEHLPW
LGTLTKFVMENNYFTSGKRLFRQIQGTAMGTNMAPEFANIYMGVAEKQKLTEVVNMGSCLYLRYIDDVLIIGTKAELDRV
AEQEHFLDLPLEVNWEKPAKSLPFLDLQLSISEHSIVYELYRKPINTYQYLEWDSDHPRATKKGFVIGEIVRIFRNCSER
ATALKHVNFFHQKIRARGYPPRLTNTWIKAAIDSIQERELGHRQAPPAENSALKIFRLKVKYSPSWEIFDAKVLRKRIEQ
LGIRRHAKRARSGRNRRLRPRKRHVRRTVVRRRKRARGNTDPKAPKRLRYQPTLDSLRGVHTS
>QBM85394.1 hypothetical protein METSCH_A00110 [Metschnikowia aff. pulcherrima]
KDSNLNPVGVNFFPAYIPPNTLMYHSREEPNIPELFEWLAMDYEFSYSFARFTRGKSRGGPHKPPKDPGDPGQPGIPGGR
SLRKDPPEKNKGLSLDEKTQKRDHLGAGPSYLFTFRNKTPLNRLILLDGASAAKSSTGEMDQQMILSQQKDVEGYVDEEI
AANNICNWGKLFGLQGVIRLEVGFEIILCDFHKDIELVSNVTLSNVTKMVGFPYEPSEPTTDLEKSRTKLIDRWQAMSEF
EWMQAGGRVNNGDGRILIDYANMVTPLNKTWVNSDPYQRRINKLSEPLKSQIITQLENTILKGVDPFYKTDWPLVIGAIV
EKFSPMLVGVNATFTVFEHEANINLSSALKNLSDNLSSSTFNFVRRYTDGNAHEAIQGTCARRDAVEDYVHNTFAISTTL
CAVPGWGPGPMGWSFGGRKSLHFENGSYRINNELECINVDQLRLR

0
Entering edit mode

Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.

Thank you!

0
Entering edit mode

As I was saying above, I tried that but when I pasted it in, the formatting messed up. Apologies, I'm relatively new to the site and don't know how to fix the issue.

0
Entering edit mode

This sequence works for me in same was as my fake example. Are you using one ID per line in file where you are specifying ID's you want to keep? What OS are you working on?

0
Entering edit mode

Given that your worked example is working fine, I suspect there's something wrong with the files I've got.

A00120
A00140
A00170
A00280


This is directly copied from my file. As above, I don't understand the formatting for this forum but should I assume that the values aren't separated by a new line but rather a tab?

EDIT: Sorry, forgot to also say I'm on MacOS Catalina 10.15.7

EDIT2: Ok, nope, the lines are definitely separated by a newline

1
Entering edit mode

Ok I got it working. I had two trailing newlines at the end of the id file which meant grep -f was returning everything rather than lines that matched my search terms.

Thanks GenoMax!

0
Entering edit mode

Just tested this on macOS big sur (v.11.3.1). Code works for me.