1
$\begingroup$

Basically a GenBank file consists of gene entries (announced by 'gene') followed by its corresponding 'CDS' entry (only one per gene) like the two shown here below. I would like to extract part of the data from the input file shown below according to the following rules and print it in the terminal.

  • There are two blocks of gene data shown below. I want to extract part of both blocks.

  • The extracted text for each block starts with a line that contains spaces at the beginning of the line followed by gene

  • The extracted text for each block ends with a line that contains /db_xref="GeneID

  • The perl and awk tags are just suggestions. An answer can use a different program(s).

Input file (source: biotech):

     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
                     /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"
                     /translation="MQTVNVDVAIVGAGGGGLRAAIAAAEANPNLKIALISKVYPMRS
                     HTVAAEGGAAAVAKEEDSYDKHFHDTVAGGDWLCEQDVVEYFVEHSPVEMTQLERWGC
                     PWSRKADGDVNVRRFGGMKIERTWFAADKTGFHLLHTLFQTSIKYPQIIRFDEHFVVD
                     ILVDDGQVRGCVAMNMMEGTFVQINANAVVIATGGGCRAYRFNTNGGIVTGDGLSMAY
                     RHGVPLRDMEFVQYHPTGLPNTGILMTEGCRGEGGILVNKDGYRYLQDYGLGPETPVG
                     KPENKYMELGPRDKVSQAFWQEWRKGNTLKTAKGVDVVHLDLRHLGEKYLHERLPFIC
                     ELAQAYEGVDPAKAPIPVRPVVHYTMGGIEVDQHAETCIKGLFAVGECASSGLHGANR
                     LGSNSLAELVVFGKVAGEMAAKRAVEATARNQAVIDAQAKDVLERVYALARQEGEESW
                     SQIRNEMGDSMEEGCGIYRTQESMEKTVAKIAELKERYKRIKVKDSSSVFNTDLLYKI
                     ELGYILDVAQSISSSAVERKESRGAHQRLDYVERDDVNYLKHTLAFYNADGTPTIKYS
                     DVKITKSQPAKRVYGAEAEAQEAAAKKE"

Desired output:

     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

I tried using pcregrep --multiline .*'START-SEARCH-TERM.*(\n|.)*END-SEARCH-TERM' path/to/SOURCE-FILE. I installed pcregrep (grep utility that uses Perl-style regexps) in Ubuntu with sudo apt install pcregrep. I also installed Biopython with sudo apt install python3-biopython and ran the Simple GenBank parsing example from Biopython Tutorial and Cookbook.

Improve this question
$\endgroup$
3
  • $\begingroup$ I would strongly suggest simply using biopython, bioruby or biojulia etc. These libraries are really good for extracting data from genbank files. You tagged perl, bioperl can definitely do it, but Perl is a dying language and bioperl is not straightforward to install IIRC from 10 years ago —but worth a try. $\endgroup$
    – Matteo Ferla
    Commented Mar 5, 2020 at 23:04
  • $\begingroup$ @MatteoFerla take that back! Just because young whippersnappers today don't appreciate the power and beauty of Perl does not make it a dying language! Them's fighting words! :P $\endgroup$
    – terdon ♦
    Commented Jun 14, 2020 at 13:52
  • $\begingroup$ Yeah agreed, code is code. Two things will continue Perl in any age, regex and Perl one liners (definitely stylish). Have you ever heard of a Python one-lliner? $\endgroup$
    – M__ ♦
    Commented Jun 14, 2020 at 18:30

2 Answers 2

Reset to default
1
$\begingroup$

You can simply use grep for this purpose as shown below,

grep -v /translation bio.txt | grep -B100000000 /db_xref= > output_file.txt

Just make sure that you keep the number with B bigger than the number of lines of your file.

If you print the contents of the above file you get your desired output as given below,

:~$ cat output_file.txt 
gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"
Improve this answer
$\endgroup$
1
$\begingroup$

You can do this very easily with awk:

$ awk '/^ {5}\w/{a=1} /\/translation/{a=0}a' file.gb 
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

The idea here is to set a to 1 if this line starts with 5 spaces followed by a word character. Then, we set a back to 0 if this line matches /translation. The default action for awk when an expression evaluates to true (not 0) is to print, therefore the final a will cause all lines read while a is not 0 to be printed, effectively removing everything after each /translation line.

A more easily understandable version of the same code would be:

awk '{ 
        if(/^ {5}\w/){
            a=1
        }
        if(/\/translation/){
            a=0
        }
        if(a==1){
            print
        }
    }' file.gb
Improve this answer
$\endgroup$
2
  • $\begingroup$ @karel it prints /db_xref="GI:219870283" but not /db_xref="GeneID:7278620". $\endgroup$
    – terdon ♦
    Commented Jun 14, 2020 at 13:50
  • $\begingroup$ @karel see update. Awk to the rescue! $\endgroup$
    – terdon ♦
    Commented Jun 14, 2020 at 13:59

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.