Добавить функцию seauence в файл genbank с помощью biopython - PullRequest
Новая
       64

Добавить функцию seauence в файл genbank с помощью biopython

0 голосов
/ 28 февраля 2020

Я новичок в python и био python, поэтому, пожалуйста, потерпите меня, если я спрошу, что это действительно глупо или абсурдно; P

Так что я работаю над групповым проектом школы, меня попросили написать файл genbank, который должен содержать:

  • для каждого контига: имя, циркуляр или нет, количество белков, CG%, таксономия c категории, геном
  • для каждого белка: ортологическая группа, цепь, координата, нуклеотидная c последовательность и последовательность белка и, наконец, таксономическая c аннотация.

Вот мой код 

################################################################################
################################################################################
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')

вывод такой: 

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

Может кто-нибудь подскажите, пожалуйста, как добавить последовательность функций с использованием bio python функция? И как я могу удалить некоторые из нежелательных разделов, таких как VERSION. Также я заметил, что bio python создает файл genbank с датой, но это неверно, и я хотел бы от него избавиться.

Спасибо в adanve за всю вашу помощь.

...