# Parsing a GenBank file with multiple gene entries

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
ILVDDGQVRGCVAMNMMEGTFVQINANAVVIATGGGCRAYRFNTNGGIVTGDGLSMAY
RHGVPLRDMEFVQYHPTGLPNTGILMTEGCRGEGGILVNKDGYRYLQDYGLGPETPVG
KPENKYMELGPRDKVSQAFWQEWRKGNTLKTAKGVDVVHLDLRHLGEKYLHERLPFIC
ELAQAYEGVDPAKAPIPVRPVVHYTMGGIEVDQHAETCIKGLFAVGECASSGLHGANR
LGSNSLAELVVFGKVAGEMAAKRAVEATARNQAVIDAQAKDVLERVYALARQEGEESW
SQIRNEMGDSMEEGCGIYRTQESMEKTVAKIAELKERYKRIKVKDSSSVFNTDLLYKI
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.

• 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. – Matteo Ferla Mar 5 '20 at 23:04
• @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 – terdon Jun 14 '20 at 13:52
• 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? – M__ Jun 14 '20 at 18:30

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"  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

• @karel it prints /db_xref="GI:219870283" but not /db_xref="GeneID:7278620". – terdon Jun 14 '20 at 13:50
• @karel see update. Awk to the rescue! – terdon Jun 14 '20 at 13:59