1

I'm new to python and biopython, so please bear with me if i ask sth really stupid or absurd ;P

So I'm working on a group projet of school, i have been asked to write a genbank file which must contains :

  • for each contigs : name, circular or not, number of proteins, CG%, taxonomic catagories, genome
  • for each protein : orthology group, strand, coordinate, nucleotic sequence and protein sequence and at last taxonomic annotation.

Here is my code 

################################################################################
################################################################################
from collections import defaultdict
import re
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from Bio.SeqFeature import SeqFeature, FeatureLocation
################################################################################
################################################################################

################################################################################
dict_idContig_genome=defaultdict(str) # done
dict_idContig_Circ=defaultdict(str) # done
dict_idContig_nbProteine=defaultdict(int) # done
dict_idContig_pourcentageCG=defaultdict(str) # done
################################################################################

with open ("contigs_circ.fna","r") as f :
    for li in f :
        li=li.rstrip('\n')
        if li.startswith('>'):
            li=li.lstrip('>')
            circ=re.search('circ',li)
            ls=li.split()
            idContig=ls[0]
            dict_idContig_nbProteine[idContig]=0
            if circ :
                dict_idContig_Circ[idContig]="Circular"
            else :
                dict_idContig_Circ[idContig]="Not circular"
        else :
            dict_idContig_genome[idContig]=li

for idContig in dict_idContig_genome :
    CG=re.findall('C|G',dict_idContig_genome[idContig],re.I)
    nb_CG=len(CG)
    nb_nucloTotal=len(dict_idContig_genome[idContig])
    percentage_CG=str(round(nb_CG/nb_nucloTotal*100,2))+"%"
    dict_idContig_pourcentageCG[idContig]=percentage_CG

################################################################################
dict_idGene_brinCodant=defaultdict(int) #
dict_idGene_coordonne=defaultdict(list) #
dict_idGene_seqNucleo=defaultdict(str) #
################################################################################

with open ('contigs_circ_prot.fna','r') as f :
    for li in f :
        li=li.rstrip('\n')
        if li.startswith('>') :
            li=li.lstrip('>')
            ls=li.split()
            idContig_nProt_partiality=ls[0]
            idGene=idContig_nProt_partiality.split('_')[0]+'_'+idContig_nProt_partiality.split('_')[1]
            idContig=idContig_nProt_partiality.split('_')[0]
            dict_idContig_nbProteine[idContig]+=1
            start=int(ls[1])
            end=int(ls[2])
            brin=int(ls[3])
            dict_idGene_coordonne[idGene]=[start,end]
            dict_idGene_brinCodant[idGene]=brin
        else :
            dict_idGene_seqNucleo[idGene]=li

################################################################################
dict_idProt_seqProt=defaultdict(str)
dict_idProt_brinCodant=defaultdict(int)
dict_idProt_coordonne=defaultdict(list)
################################################################################

with open ('contigs_circ_prot.faa','r') as f :
    for li in f :
        li=li.rstrip('\n')
        if li.startswith('>') :
            li=li.lstrip('>')
            ls=li.split()
            idContig_nProt_partiality=ls[0]
            idProt=idContig_nProt_partiality.split('_')[0]+'_'+idContig_nProt_partiality.split('_')[1]
            partiality=idContig_nProt_partiality.split('_')[2]
            start=int(ls[1])
            end=int(ls[2])
            brin=int(ls[3])
            dict_idProt_coordonne[idProt]=[start,end]
            dict_idProt_brinCodant[idProt]=brin
        else :
            dict_idProt_seqProt[idProt]=li

################################################################################
################################################################################

for idContig in dict_idContig_genome :
    seqNucleo_string = dict_idContig_genome[idContig]
    seqNucleo_objet = Seq(seqNucleo_string, IUPAC.unambiguous_dna)
    record = SeqRecord(seqNucleo_objet, id=idContig, name="<unknown name>", description=dict_idContig_Circ[idContig],annotations=None, letter_annotations=None)
    for idGene in dict_idGene_seqNucleo :
        if idGene.split('_')[0] == idContig :
            feature_gene=SeqFeature(FeatureLocation(start=dict_idGene_coordonne[idGene][0],end=dict_idGene_coordonne[idGene][1]),type="gene",strand=dict_idGene_brinCodant[idGene])
            record.features.append(feature_gene)
            feature_protein=SeqFeature(FeatureLocation(start=dict_idProt_coordonne[idGene][0],end=dict_idProt_coordonne[idGene][1]),type="protein",strand=dict_idProt_brinCodant[idGene])
            record.features.append(feature_protein)
    output_file = open ('test.gb','a')
    SeqIO.write(record,output_file,'genbank')

the output is like this 

LOCUS       c7                     10849 bp    DNA              UNK 01-JAN-1980
DEFINITION  Not circular.
ACCESSION   c7
VERSION     c7
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
     gene            4..1046
     protein         4..1046
     gene            complement(1048..5231)
     protein         complement(1048..5231)
     gene            complement(5368..7175)
     protein         complement(5368..7175)
     gene            7233..7777
     protein         7233..7777
     gene            7762..8420
     protein         7762..8420
     gene            complement(8429..9072)
     protein         complement(8429..9072)
     gene            9112..10790
     protein         9112..10790
ORIGIN
        1 ataaggtaaa aaaaaccgac aggttcggta aaacggaaaa aacctcggta aaagtgtgcg
       61 taaaggttat tttcaataac tttattgata tgtcagctga aaaaagaccc cgtttgtttg
      121 gtgagattta tagtatagga aaagaaaccc ttccaccacc accaggtgtt agtaaagaag
      181 cgtttaaaag aatgcttgag ttttgtgaaa gtggtgttcc tggaacgaaa ccctactact
      241 acagatgggc actttggttg tatagagtac ttgcagaacg aggtagtaaa aaatctgatt
      301 ttgaggcttt tatgaagagg ctacaatata ttcttgatgg tgttttaggt atttatgctg
      361 atgcgattta tccaagctct atttcttctc catatctttt ggaggatatt ccttcaatta
      421 gaaggcaagt tgatgagttc atgagacgtg aaggttgggg gtctttacta gagtttcttc
      481 attcatattc agatgatttt gaggatcctg atttcattaa tagaggggat tttactatac
      541 gtttttcttt cctgttgcct tttgttttac ttgcgttgca gaataattat aattttgata
      601 tgattcctcc tggttcgtta agaagtaagt tcgttatgta tgagcctctt gttatagaaa
      661 atgctcataa gttgtacatg aaaagtttgt ttggtgagga cccaacagga aaacgtaagt
      721 tacatgttcc tacatgggct tatccttggc accgtaaagg ccctgttcct gaagcaagat
      781 tagatgccga agacgccttt ataaatgata caaaacaaat ttcagaatat agagactttg
      841 atctttggca gaattatgaa gagcttttac aacgaaataa ggcagaagca tttaaaagat
      901 tgcttgccga tgatcctaat gctgttactt atatatgggc tggtacacat gcaacagggt

Can someone please tell me how to add feature sequence using biopython function ? And how can i delete some of the unwanted section such as VERSION. Also i have noticed that biopython create the genbank file with a date, but it's incorrect and i would like to get rid of it.

Thanks in adanve for all your help.

mewu3
  • 11
  • 1

0 Answers0