0
$\begingroup$

sorry for the question, I'm trying to extract the proteins sequences from a genbank file.

 gene            complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
 CDS             complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_011081430.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="MARTX multifunctional-autoprocessing
                 repeats-in-toxin holotoxin RtxA"
                 /protein_id="WP_011081430.1"
                 /translation="MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHV
                 TVGSIGATVYTGSGNDTVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGA
                 AGGVSIDHLGNHGDVSYGGAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFA
                 GAGAGNKLDRTWFNRYQGSRGDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKG
                 KEGNHTANLANEDISSANGYHSMGKGGYSLSDLHYSVNAVRSTSETVADIDEYTDQTL
                 FKPATDSGESSGDVRFNGAGGGNVIKSNVTRGNVYFNGGGIANVILHSSQFGHTEFNG
                 GGAANVIVKSGEEGDLTFRGAGLANVLVHQSKQGKMDVYAGGAVNVLVRIGDGQYLAH"

I'm using a code from a previous answer to extract the protein and the locus_tag, but I want some extra features like:

>locus_tag|old_locus_tag|gene|product|complement|protein_id|chromosome|Accession|organism|strain
MGKPFWRSVEYFFTGNYSADDG ..........

the result will be like:

>VV1_RS17390|VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA"|516466..532086|WP_011081430.1|Chromosome-II|NC_004459.3|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDG ..........

this is the code that I'm using, how to add those features, or an other way using R ???

from Bio import SeqIO
file_name = 'CMCP6.gb'
# stores all the CDS entries
all_entries = []
with open(file_name, 'r') as GBFile:
    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)


# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

But the result from the code only add the locus_tag and the sequence:

>VV1_RS17380
MFIKNLSIGKKIAAAFSIIAVINIAFGIFLSTELNTVKSELLNYTEDTLPAMEKVDAVRD
KISYWRRTQFAVFAMNDENQIKQTITRNEGIRREIETELAAYGKSVWPGEEEQTYNRLMS
LWSGYLSTMDKFNDALLAGDKDAAYPILTNSLSTFESIETEVNKLVMILKGAMDSNKNQI
LSSVNGLNTTAVISNIAIFAVMVIMTLLLTRIICGPLNIVVRQANAIARGDLSQDLDRNA
IGNDELGELADATIKMQDDLRQVIDNVIAAVTQLSSAVEEMTQISEMSASGMKDQQMQVT
LVATAMTEMKAAVADVARNTEDSASQAYDANRRTQDGAKETHQMVVSIEEVADIIAKAGD
TVAELEAQSNQINVVVDVIRGIADQTNLLALNAAIEAARAGESGRGFAVVADEVRTLAGR
TQDSTGEITSIIEKLQDLAKQAKHATEDSRTSIGACVEQGNNAQELMGSIEKSIANIADM
GAQIATACGQQDSVAEELSRNIENIHMASQEVAQGSQQTAQACRELTQLSVSLQDVMSRF
KLN
>VV1_RS17385
MNFKKTLLSIAIASASLTPAFSYSAPLLLDNTVHQTSQIAGANAWLEISLGQFKSNIEQF
KSHIAPQTKICAVMKADAYGNGIRGLMPTILEQQIPCVAIASNAEAKLVRESGFEGELIR
VRSASTSEIEQALSLDIEELIGSEQQARELASLAEKYSKTIKVHLALNDGGMGRNGIDMS
TERGPKEAVAIATHPSVAVVGIMTHFPNYNAEDVRTKLKSFNQHAQWLMESAGLKREEIT
LHVANSYTALNVPEAQLDMVRPGGVLYGDLPTNPEYPSIVAFKTRVASLHSLPAGSTVGY
DSTFTTANDAVMANLTVGYSDGYPRKMGNKAQVLINGQRANVVGVASMNTTMVDVSNIKG
VLPGDEVTLFGAQKNQHISVGEMEENAEVIFPELYTIWGTSNPRFYVK
>VV1_RS17390
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND
TVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGAAGGVSIDHLGNHGDVSYG
GAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFAGAGAGNKLDRTWFNRYQGSR
GDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKGKVGDITLQGAGASNRIERTRQA
EDVYAQTRGNIRFEGVGGYNSLYSDVAHGDIHFSGGGAYNTITRKGSGSSFDAQGMEYAK

Thanks so much

$\endgroup$
1
  • $\begingroup$ cds.description = '' $\endgroup$
    – pippo1980
    Nov 2, 2021 at 22:36

1 Answer 1

3
$\begingroup$

ok tried a little bit, not sure is the most elegant/fastest way of doing it

starting from file from : https://www.ncbi.nlm.nih.gov/nuccore/NC_004460.2?report=gbwithparts&log$=seqview {GenBank(Full)}

here my code:

from Bio import SeqIO

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

file_name = 'CMCP6.gb'


# stores all the CDS entries
all_entries = []


gb_record = SeqIO.read(open(file_name,"r"), "genbank")


gb_feature = gb_record.features


accession = gb_record.annotations['accessions'][0]+'.'+str(gb_record.annotations['sequence_version'])

# accession = gb_record.annotations['accessions']
# print(gb_record.annotations['accessions'])

for seq_feature in gb_feature:
    if seq_feature.type=="source":
            organism = str(seq_feature.qualifiers['organism'][0]).rsplit(" ", 1)[0]
            strain = seq_feature.qualifiers['strain'][0]
            chromosome = seq_feature.qualifiers['chromosome'][0]



for seq_feature in gb_feature:
    if seq_feature.type=="CDS":
            if 'translation' in seq_feature.qualifiers:
                if seq_feature.qualifiers['translation'][0] != '':
                    if 'gene' in seq_feature.qualifiers:
                        gene = seq_feature.qualifiers['gene'][0]
                
                    else:
                    
                        gene = 'not defined'
                    
                    if 'product' in seq_feature.qualifiers:
                        product = seq_feature.qualifiers['product'][0]
                    
                    else:

                        product = 'not defined'
                        
                    if 'protein_id' in seq_feature.qualifiers:
                          protein_id = seq_feature.qualifiers['protein_id'][0]
                     
                    else:

                          protein_id = 'not defined'
                         
                         
                    if 'old_locus_tag' in seq_feature.qualifiers:
                          old_locus_tag = seq_feature.qualifiers['old_locus_tag'][0]
                     
                    else:

                          old_locus_tag = 'not defined'
                         
                    if 'locus_tag' in seq_feature.qualifiers:
                          locus_tag = seq_feature.qualifiers['locus_tag'][0]
                     
                    else:

                          locus_tag = 'not defined'
            
                    if seq_feature.location.strand == 1:
                          complement = 'no complement'

                    if seq_feature.location.strand == -1:
                        complement = str(seq_feature.location).strip('[]').split('(-)')[0]
                        
                        
                        
                    pippo = SeqRecord(
                                Seq(seq_feature.qualifiers['translation'][0]),
                                
                                id = seq_feature.qualifiers['locus_tag'][0] ,
                                
                                description = ('|'+str(old_locus_tag)+'|'+str(gene)+'|'+str(product)\
                                                +'|'+str(complement)+'|'+str(protein_id) +'|Chromosome-'+str(chromosome)\
                                                +'|'+str(accession).strip("[]").replace("'", "")+'|'+str(organism)\
                                                    +'|'+str(strain))
                                )
                    all_entries.append(pippo)
                                
# print(all_entries)
# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

output for your example line is :

>VV1_RS17390 |VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA|516465:532086]|WP_011081430.1|Chromosome-II|NC_004460.2|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND.....

I got extra space between locus_tag and | , I dont know how to get rid of it

Its not clear to me what you mean by Accessions, tried two ways (second one is past # in code) and got NC_004460.2 instead of NC_004459.3

$\endgroup$
1
  • $\begingroup$ thanks so much it works perfectly, I really appreciate your help !!!! $\endgroup$
    – abraham
    Nov 5, 2021 at 6:51

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.