0
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ 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? $\endgroup$
    – terdon ♦
    Commented Jul 9, 2019 at 17:08
  • $\begingroup$ 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. $\endgroup$
    – R-MASHup
    Commented Jul 9, 2019 at 19:22
  • $\begingroup$ 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. $\endgroup$
    – terdon ♦
    Commented Jul 9, 2019 at 19:30
  • 1
    $\begingroup$ 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. $\endgroup$
    – Daniel Standage
    Commented Jul 9, 2019 at 19:47
  • $\begingroup$ 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. $\endgroup$
    – terdon ♦
    Commented Jul 9, 2019 at 20:47

1 Answer 1

Reset to default
1
$\begingroup$

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
Improve this answer
$\endgroup$
0

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.