Subset a protein fasta based on sequence length?

Question

I have a strange protein fasta file. There are several entries for the same gene and I need to extract the entry with the longest width.

>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 a
lpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..8944
8))] [gbkey=CDS]
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFGSTGFIGRYVANKLGKIGTQL
ILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKALKYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAG
VEKFIHLSALNCDGPNEAIFSRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWGYKRYDMKYDMFFKLKVSLTQ
KFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLRRMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017771299.1_10 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017771299.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SAVVCGFPRDEGGCMSNVDRWYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQ
KILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHR
TEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017772069.1_11 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X2] [protein_id=XP_017772069.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109546..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD
DGPLSKILCEEVQENEDTQPNMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMSLVSMIDPSPDWIVGVSGLEL
CLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHPIRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEK
NCVAQVDVSEEDGGVEGDKCEMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERECGPSDRRPLQESKECEGRAW
SNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDTEDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDN
SVTENYFGRIVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGFSTKHRFIRVHPSNGGKRCPQKILKKR
KCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCVMTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPC
SLLACPE
>lcl|NW_017095468.1_prot_XP_017772815.1_12 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X3] [protein_id=XP_017772815.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,109573..109812,109882..110114)] [gbkey=CDS]
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNIFLRSNGEYQAKNKFKDFLLL
VEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVMQTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTD

I've tried using BioStrings but I'm unfamiliar with how to index this type of file. Any help appreciated.

Can you give us an example that includes a duplicated gene? How is the "gene" defined? Should we look for the RefSeq protein ID (XP_123), for the gene=LOC string? All of the above? Will the entire header line be identical or only parts of it? — terdon, Jul 09 '19 at 17:08
I've added more of the file to my question. I think the unique/isoform bits are the XP_ (which has information about isoforms) ... NW_ bits are the genes. — R-MASHup, Jul 09 '19 at 19:22
NW_017095468 is a whole genome shotgun sequence, it will have many genes. The XP bits are RefSeq protein accessions. These can be different for different protein products of the same gene. It looks like the solutions will need to look at either the gene= or the GeneID: values. — terdon, Jul 09 '19 at 19:30
Perhaps more important is for R-MASHup to explain what is "strange" about these sequences and why. It is well known that a gene may encode multiple proteins, so the fact that there are multiple entries per gene doesn't surprise me. — Daniel Standage, Jul 09 '19 at 19:47
Yes. And choosing the longest is kind of arbitrary, but all this depends on what you're doing. For instance, I once needed to do something very similar because I was looking at various features of certain genes and I wanted to see if protein product length correlated with some other characteristic, so I chose the longest isoform of each gene's products. But, @R-MASHup, if you think the longest will be somehow special (most highly expressed or more representative or anything like that), you will be surprised. — terdon, Jul 09 '19 at 20:47
Well it's a nonmodel organism, trying to find orthologs, and so need to reduce the number of entries somehow (esp as it's computationally so slow) ... how would you pick? — R-MASHup, Jul 09 '19 at 20:58
For orthologs, I'd probably not pick since it's impossible to predict which protein will have ended up shared with the other organism. I'd just run a tblastn or better exonerate on using all proteins as queries. But yeah, picking the longest might make sense if you don't have a decent machine. Also, if my answer isn't working for you, you probably have ID lines that are different from what you've shown here. My solution just takes the 2nd field of the line ([gene=LOCNNNN] in your example), if that changes, it will not work. — terdon, Jul 11 '19 at 08:24

score 1 · Accepted Answer · answered Jul 09 '19 at 19:46

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do:

FastaToTbl file.fa | 
    awk '{
            len=length($NF); 
            if(len > l[$2]){
                l[$2]=len; 
                seq[$2]=$0; 
            }
          } 
          END{
            for(s in seq){
                print seq[s]
            }
        }' | TblToFasta

Running this on your example gives:

$ FastaToTbl file.fa | awk '{len=length($NF); if(len > l[$2]){l[$2]=len; seq[$2]=$0; }} END{for(s in seq){print seq[s]}}' | TblToFasta 
>lcl|NW_017095468.1_prot_XP_017770525.1_9 [gene=LOC108557919] [db_xref=GeneID:108557919] [protein=spondin-1 isoform X1] [protein_id=XP_017770525.1] [location=join(95471..95648,106560..106774,106834..107025,107093..107299,107562..107816,107872..108474,108526..108799,108885..109058,109127..109237,109573..109812,109882..110114)] [gbkey=CDS] 
MRLKVAFLWLVSSISWIGEALRCDRTPEGTFSPRTRADGRFVIEVSGNPDTYVPGEQYNI
FLRSNGEYQAKNKFKDFLLLVEHEPSEKILGEVHNPSVGTLQLLGDMLMKFSEKCRNAVM
QTNSLPKSEVQVLWVAPPSGSGCVAIRATVVESKEFWYTDDGPLSKILCEEVQENEDTQP
NMLRQCCACDEAKYEVTFEGLWSRNTHPKDFPSNGWLTRFSDIIGASHTFDYTFWNYGEI
ASNGLRQLAENGNTRMLESELKAKSEHIRTIIKARGISYPNITGKTFAVFRVDKRHHLMS
LVSMIDPSPDWIVGVSGLELCLRNCSWVESRVLNLYPWDAGTDDGPTYISANQPSMPPHP
IRRIKSNSPNDPRSPFYDPTGTEMKPLARLYLSRQRLYEKNCVAQVDVSEEDGGVEGDKC
EMEEWSEWSKCTVTCGRGFKYKQRAYKNPASNFVCNKPLTKRASCVAILEHCSNQQRPQE
ADPSCSLTGWGNWSSCTAPCGPGWKTRSRRYKNRKAAKRCAAGNENPEPLEQNLECMERE
CGPSDRRPLQESKECEGRAWSNWSPCSSTCGKGIKVRRRMAYRSLWGRSPARYNRGLFDT
EDTSRDDDDDGSDEDPCMNLDEKVECINDDVPVCEDTVDNSAVVCGFPRDEGGCMSNVDR
WYFDVIKGNCDIFSYSGCQGNKNNFKTLERCENVCDSYKKELLANRTAYKRQLGVTVSGV
LSYNLHHMQNDDADNCVPGSQTRQDQDKKIIQEPIGEVVDCQMSEWTNWSGCNATCGRGF
STKHRFIRVHPSNGGKRCPQKILKKRKCKIPCPGDYTKRDPMLPTWGTANSLEHVQIDCV
MTGWSAWSPCSRSCGPNAVQQRTRGILLPPSGRGEPCLHRTEERPCSLLACPE
>lcl|NW_017095468.1_prot_XP_017786561.1_8 [gene=LOC108569215] [db_xref=GeneID:108569215] [protein=NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial] [protein_id=XP_017786561.1] [location=complement(join(88030..88986,89044..89243,89412..89448))] [gbkey=CDS] 
MAAVIFTGAQLLKQQSGLVGIAYIRVNNYSSDAKVYNLASLKRGTGGRSSFNGIVATVFG
STGFIGRYVANKLGKIGTQLILPYRGNNYETMRLKLCGDLGQVLFQPYYLKDLESIDKAL
KYSNVVINLVGRDWETKNFSFQDVHVKGARDLARAAKKAGVEKFIHLSALNCDGPNEAIF
SRQGSKFLSSKWEGEQAVLEEFPEATIFRPSDVYGQEDRFLRYYAHNWRRQGQYMPMWKN
GEATIKQPVHVSDLAAGIVAAIKDPEAAGKVYQAVGPRRYQLNELVDWFYRVMRKDSEWG
YKRYDMKYDMFFKLKVSLTQKFSPAFPVGNLHWERLEREFVTDRVNQAIPTLEDLGVNLR
RMEDQVPWELKPFTYGLYFGTDSEEPVVEPKPPKYVS

Subset a protein fasta based on sequence length?

1 Answers1