我是 python 和 biopython 的新手，所以如果我问的东西真的很愚蠢或荒谬，请多多包涵；P
所以我正在研究一个学校的小组项目，我被要求编写一个必须包含的 genbank 文件：
- 对于每个重叠群：名称、环状与否、蛋白质数量、CG%、分类学类别、基因组
- 对于每个蛋白质：直系群、链、坐标、核序列和蛋白质序列，最后是分类注释。
这是我的代码
################################################################################
################################################################################
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
有人可以告诉我如何使用 biopython 函数添加特征序列吗？以及如何删除一些不需要的部分，例如 VERSION。我也注意到 biopython 创建带有日期的 genbank 文件，但它不正确，我想摆脱它。
感谢您的所有帮助。