How can I pull out specific protein fastas from one file using information from the protein header?

1

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

example of my .faa file

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:

list of gene IDs

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

• 85 views

$ more prot.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

$ more id
A00100
A20100

# linearize fasta (code from @Pierre) 

$ awk '/^>/ {printf("%s%st",(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%st",(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


Login
before adding your answer.

Traffic: 2374 users visited in the last hour



Source link