How to get the product of a CDS

Question

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

Answer 1

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

Answer 2

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

Answer 3

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.