3
$\begingroup$

I need the name of the protein in /product="protein_name" using bash commands. Beware, there is a lot of whitespace between lines.

FEATURES             Location/Qualifiers
     source          1..1266
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /isolate="Salem"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
                     /country="Australia: Taronga Zoo, Sydney"
     gene            1..1266
                     /gene="UROS"
                     /note="Derived by automated computational analysis using
                     gene prediction method: Gnomon. Supporting evidence
                     includes similarity to: 10 Proteins, and 100% coverage of
                     the annotated genomic feature by RNAseq alignments,
                     including 5 samples with support for all annotated
                     introns"
                     /db_xref="GeneID:100933602"
     CDS             60..908
                     /gene="UROS"
                     /codon_start=1
                     /product="uroporphyrinogen-III synthase isoform X1"
                     /protein_id="XP_003755055.2"
                     /db_xref="GeneID:100933602"
                     /translation="MRIPTSFSRGAWWKMKVLLLKDAKDDDSGLDPYLQELGSCGFEA
                     TLIPVLAFEFLSLPSFSEKLSHPEAYGGLVFTSPRAVEALQLSLERAGQLEVWKKSLK
                     EKWRLKPAYVVGQATASLVKKVGLAAEGEKCGNGEKLAEYICSRETPATLPLLFPCGT
                     LKGEVLPKTLKQKGILLESLTVYQKVQHPRMQDSLKSYFSTQGVPAGVTFFSPSGLTY
                     SLPHIQELSGASFDRIKFAAIGPTTARALASAGIPVSCTAEKPCPQDLAAGLQKALQP
                     SGPSSS"

Right now I have this:

while read p;
do cat "${p}.gb" | tr "\"" " "| egrep -v "/db_xref= taxon"| egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk '{print($2)}END{print("")}'
done <XM_codes.txt

Which results in this:

GeneID:100933602
UROS
uroporphyrinogen-111
XP_003755055.2

My output should be:

GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
$\endgroup$
0
$\begingroup$

If you have GNU grep (default on Linux), you can do this with a single grep command using Perl Compatible Regular Expressions (PCREs):

$ grep -oP 'GeneID:\d+|/(gene|product|protein_id)="\K[^"]+' file | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

The -o flag tells grep to only print the matching portion of the line and the -P enables PCREs which let us use d+ for "one or more digits" and \K which means "discard whatever was matched before the \K".

Note that I'm not using cat here. The command from your question is what is known as a UUoC (useless use of cat) and is a surefire way of annoying any *nix geeks around ;).


If you don't have GNU grep, you can do something like this:

$ awk -F= '/GeneID:|(gene|product|protein_id)=/{print $2}' file | 
    tr -d '"' | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2
$\endgroup$
  • $\begingroup$ Thank you so much! Also, what to do if the "uroporphyrinogen-III synthase isoform X1" string was on more than one line? $\endgroup$ – Aron van Beelen Sep 17 '18 at 16:13
  • $\begingroup$ @AronvanBeelen ah, that is a very, very different situation. Is that even a valid genbank file? How would it look like? You would also need to remove the spaces on the second line, right? That is approaching the level of complexity where it is better to write a proper parser instead of hacking something together on the command line. You might want to ask a new question, showing an example that accurately represents your real data. $\endgroup$ – terdon Sep 17 '18 at 16:18
  • $\begingroup$ I've asked a new question, I hope you can help me! :) $\endgroup$ – Aron van Beelen Sep 17 '18 at 16:37
3
$\begingroup$

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t"). In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter.

 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

results:

GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
$\endgroup$
0
$\begingroup$

The solutions here are more robust than mine since I rely on the structure of the file, but if you want a solution that is really easy to understand:

$grep -A 5 "CDS" file | grep -v "CDS" | grep -v "codon_start" | cut -f 2 -d \"
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
GeneID:100933602

I take CDS line and 5 lines bellow, then I kick out CDS and codon_start (lines we are not interested in) and then I just bring what is between quote using cut. I also report the values in the order as they appear in the original file.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.