This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to extract certain CDS from a GenBank file in linux terminal

Hi, I have a large GenBank file that contains multiple records. I got this file after performing BLAST against my query sequence. This .gb file contains all the results (full record in GenBank format for each accession) from the blast search. I want to extract a certain CDS from each record in that .gb file. As additional complexity, this CDS is a join of two regions in the genome. Is there any automated process by using code in the Linux terminal, so that I can extract my target CDS from all records and save them in a separate file? Thanks in advance.

gene sequence genome
• 7,553 views • link
updated by AP • AP • written by Anisur Rahman • Anisur Rahman

Please post example gb file and expected output.

• 0 views • link
written by cpad0112 cpad0112

Here I don't find any option to upload the GenBank file. But it is very long file, so do I just copy-paste it here?

• 0 views • link
written by Anisur Rahman • Anisur Rahman

Post one or two records as input and expected output for those one or two records.

• 0 views • link
written by cpad0112 cpad0112

Okay, according to your suggestion I cut off some portion of these genbank records. Here, there are two genbank record in a single file:

LOCUS       MK733608              136382 bp    DNA     linear   VRL 13-MAY-2020
DEFINITION  Human gammaherpesvirus 8 strain UNC_KICS010, complete genome.
ACCESSION   MK733608
VERSION     MK733608.1.
FEATURES             Location/Qualifiers
     source          1..136382
                     /organism="Human gammaherpesvirus 8"
                     /mol_type="genomic DNA"
                     /strain="UNC_KICS010"
     gene            1082..17050
                     /gene="ORF4"
                     /note="HHV8GK18_gp02"
     polyA_site      1113..1113
                     /gene="K1"
     CDS             1142..2794
                     /gene="ORF4"
                     /note="complement control protein; membrane protein;
                     contains four SCR domains"
                     /codon_start=1
                     /product="KCP"
                     /protein_id="QFU18707.1"
                     /translation="MAFLRQTLWILWTFTMVIGQDNEKCSQKTLIGYRLKMSRDGDIA
                     VGETVELRCRSGYTTYARNITATCLQGGTWSEPTATCNKKSCPNPGEIQNGKVIFHGG
                     QDALKYGANISYVCNEGYFLVGREYVRYCMIGASGQMAWSSSPPFCEKEKCHRPKIEN
                     GDFKPDKDYYEYNDAVHFECNEGYTLVGPHSIACAVNNTWTSNMPTCELAGCKFPSVT"
ORIGIN
        1 TACTAATTTT CAAAGGCGGG GTTCTGCCAG GCATAGTCTT TTTTTCTGGC GGCCCTTGTG
       61 TATACCTGTC TTTCAGACCT TGTTGGACAT CCTGTACAAT CAAGATGTTC CTGTATGTTG
      121 TCTGCAGTCT GGCGGTTTGC TTTCGAGGAC TATTAAGCCT TTCTCTGTTA TCGTCTCCAA
      181 CTTTGTGCCC TGGAGTGATT TCAACGCCTT ACAATTTGAC CTGTCCGTCT AATACATCCT
      241 TGCCAACATC CTGGTATTGC AACGATACTC GGCTTTTACG AGTGACGCAG GGAACATTGA
//
LOCUS       MK733606              136317 bp    DNA     linear   VRL 13-MAY-2020
DEFINITION  Human gammaherpesvirus 8 strain UNC_KICS009, complete genome.
ACCESSION   MK733606
VERSION     MK733606.1
FEATURES             Location/Qualifiers
     source          1..136317
                     /organism="Human gammaherpesvirus 8"
                     /mol_type="genomic DNA"
                     /strain="UNC_KICS009"
                     /gene="ORF4"
                     /note="HHV8GK18_gp02"
     polyA_site      1113..1113
                     /gene="K1"
     CDS             1142..2794
                     /gene="ORF4"
                     /note="complement control protein; membrane protein;
                     contains four SCR domains"
                     /codon_start=1
                     /product="KCP"
                     /protein_id="QFU18535.1"
                     /translation="MAFLRQTLWILWTFTMVIGQDNEKCSQKTLIGYRLKMSRDGDIA
                     VGETVELRCRSGYTTYARNITATCLQGGTWSEPTATCNKKSCPNPGEIQNGKVIFHGG
                     QDALKYGANISYVCNEGYFLVGREYVRYCMIGASGQMAWSSSPPFCEKEKCHRPKIEN
                     GDFKPDKDYYEYNDAVHFECNEGYTLVGPHSIACAVNNTWTSNMPTCELAGCKFPSVT
                     HGYPIQGFSLTYKHKQSVTFACNDGFVLRGSPTITCNVTEWDPPLPKCVLEDIDDPNN"
ORIGIN
        1 TACTAATTTT CAAAGGCGGG GTTCTGCCAG GCATAGTCTT TTTTTCTGGC GGCCCTTGTG
       61 TATACCTGTC TTTCAGACCT TGTTGGACAT CCTGTACAAT CAAGATGTTC CTGTATGTTG
      121 TCTGCAGTCT GGCGGTTTGC TTTCGAGGAC TATTAAGCCT TTCTCTGTTA TCGTCTCCAA
      181 CTTTGTGCCC TGGAGTGATT TCAACGCCTT ACAATTTGAC CTGTCTGTCT AATGCATCCT
      241 TGCCAATATC CTGGTATTGC AACAATACTC GGCTTTTGCG ACTGACGAAG AACACACTCC
      301 TTGCTGTCGC CATTGCCTGC AATTTTACTT GTGTGGAACA ATCTGGGCAT CGACAGAGC
//

And from this file I want to extract a cds suppose ORF4 (as fasta format) from every record. The output would be as follow (for character limitation I deleted some sequence from the fasta files):

>MK733608.1:1142-2794 Human gammaherpesvirus 8 strain UNC_KICS010, complete genome
ATGGCCTTTTTAAGACAAACACTGTGGATTTTATGGACATTTACCATGGTCATTGGCCAGGACAATGAAA
AGTGTTCCCAAAAAACCTTAATTGGATATAGACTTAAAATGTCTCGTGACGGTGACATTGCAGTTGGAGA
AACAGTGGAATTACGTTGTAGATCTGGATACACTACTTATGCCCGCAATATAACAGCAACATGTTTACAA
GGTGGGACGTGGTCTGAACCAACGGCAACATGTAACAAAAAGTCCTGTCCAAACCCAGGTGAAATACAAA
ATGGAAAGGTTATATTTCATGGTGGACAAGATGCCTTAAAATATGGGGCAAACATTTCATATGTTTGTAA
>MK733606.1:1142-2794 Human gammaherpesvirus 8 strain UNC_KICS009, complete genome
ATGGCCTTTTTAAGACAAACACTGTGGATTTTATGGACATTTACCATGGTCATTGGCCAGGACAATGAAA
AGTGTTCCCAAAAAACCTTAATTGGATATAGACTTAAAATGTCTCGTGACGGTGACATTGCAGTTGGAGA
AACAGTGGAATTACGTTGTAGATCTGGATACACTACTTATGCCCGCAATATAACAGCAACATGTTTACAA
GGTGGGACGTGGTCTGAACCAACGGCAACATGTAACAAAAAGTCCTGTCCAAACCCAGGTGAAATACAAA
• 0 views • link
written by Anisur Rahman • Anisur Rahman

try this code:

from Bio import SeqIO

    for rec in SeqIO.parse("test.gb", "genbank"):
        for feature in rec.features:
            if feature.type == "CDS":
                seq = feature.qualifiers["translation"]
                id = ">{}:{} {}".format'('rec.id, feature.location, rec.description)
                print("{}\n{}".format(id, *seq))
  1. test.gb = input genbank
  2. Please mind tabs in python (3)
  3. Current code prints the sequences to the screen. Write it to a file once you are okay with the code.
  4. Replace '(' with (
• 0 views • link
written by cpad0112 cpad0112

Thank you for your code. I run the code according to your instruction. But, it generating all amino acid sequences for every CDS. But my target is to get DNA sequence of a specific CDS (lets say, ORF4).

• 0 views • link
written by Anisur Rahman • Anisur Rahman

For nucleotide sequence, please try this:

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

_Sequences = []

for rec in SeqIO.parse("test.gb", "genbank"):
    id = "{}:{}{}".format'('rec.id, rec.features[2].location, rec.description)
    sequence = SeqRecord(rec.seq, id=id, description="")
    _Sequences.append(sequence)

SeqIO.write(_Sequences, "test.fa", 'fasta')

Replace '(' with (

• 0 views • link
written by cpad0112 cpad0112

Hey, I am trying this code but here I getting the translated sequence of all the genes. But I want a translated sequence of a particular gene only from the GB file. Can you help me out with this?

• 0 views • link
written by Gen • Gen

Hey Gen did you figure out your problem I am having the same issue

• 0 views • link
written by AP • AP

1 answer

See the BIRCH tutorial CREATING DATASETS OF RELATED SEQUENCES at http://home.cc.umanitoba.ca/~psgendb/birchhomedir/BIRCHDEV/public_html/tutorials/bioLegato/dataset/dataset.html

• 0 views • link
written by brian.fristensky brian.fristensky

Log in to answer this question.

More posts like this