如何修复将一个FASTA文件拆分为多个文件的程序中的悖论/偶然循环?

2024-05-29 04:12:04 发布

您现在位置:Python中文网/ 问答频道 /正文

我创建了一个文件,它包含一个GFF3文件和一个FASTA基因组文件,并将GFF3文件中的每个外显子写入一个新文件,将每个外显子基因组写入一个新的CDS文件。问题是,外显子写入一个单独的文件,但CDS文件只有一个头,没有CDS序列。你知道吗

我怀疑问题始于第一线

   if NMID != previous:

当我试着看外显子的父基因是否与前一个外显子相同时。如果不是,则在打开新的CDS文件之前,将写入并关闭当前的CDS文件。你知道吗

有没有关于修改从开始的if else条款的建议

   if NMID != previous:

我想这是因为我在程序的范围内设置了值previous。你知道吗

如何在将一个文件拆分为多个文件的程序中修复这个悖论/偶然循环?你知道吗

输入文件:

ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.49_FB2013_01/fasta/dmel-all-chromosome-r5.49.fasta.gz这是FASTA文件

ftp://ftp.ensembl.org/pub/release-77/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP5.77.gtf.gz这是gff3文件。尽管从技术上讲它是gtf,但这是我用作gff3文件的文件,因为它包含了外显子。你知道吗


我有一个gff3文件和一个FASTA基因组文件

gff3文件如下:

20 protein2genome exon 12005 12107 . - . ID=Blah_exon:1;Name=Blah:1;Parent=Blah

20 protein2genome exon 12108 12200 . - . ID=Blah_exon:2;Name=Blah:1;Parent=Blah

20 protein2genome exon 12005 12107 . - . ID=Blah_exon:3;Name=Blah:1;Parent=Blah

20 protein2genome exon 12342 12542 . - . ID=Blah_exon:4;Name=Blah:1;Parent=Blah

20 protein2genome exon 13005 13107 . - . ID=ABC_exon:1;Name=ABC:1;Parent=ABC

对于每个外显子,我需要编写一个名为gene\u name\u exon的新文件_编号.fa你知道吗

文件中的标题如下:

gene_name exon_number chromosome_name exon_start exon_end

然后是它染色体上的核苷酸碱基

所以第一个外显子文件看起来像:

Blah 1 20 12005 12107 ACGGGCCTCAAAGCCGCGACACGACGGCTGTCGGCCGGTAACAGTAACCCCGGAGTGAACTCCTAT

每个外显子将是一个新的文件。你知道吗

另外,我需要为每个基因创建一个包含所有名为gene\u name的外显子的CDS文件_光盘.fa你知道吗

文件中的标题如下:

gene_name number_of_exons chromosome_name cds_start cds_end cds_total_length

所以第一个CDS文件看起来像:

Blah 4 20 12005 12542 ACGGGCCTCAAAGCCGCGACACGACGGCTGTCGGCCGGTAACAGTAACCCCGGAGTGAACTCCTAT ACGGGCCTCAAAGCCGCGACACGACGGCTGTCGGCCGGTAACAGTAACCCCGGAGTGAACTCCTAT

注:核苷酸碱基只是所有外显子的核苷酸碱基粘贴在一起。你知道吗

最后,在一个包含几百个基因和几千个外显子的GFF3文件中,每个基因将有几个外显子文件和一个CDS文件。你知道吗

对于外显子,我需要这个:

gene_name exon_number chromosome_name exon_start exon_end

对于CD,我需要这个:

gene_name number_of_exons chromosome_name cds_start cds_end cds_total_length

cds\u总长度应为cds\u end-cds\u start

from collections import defaultdict
from Bio import SeqIO

Exon_file = "exon.gff3"
FASTA_file = "melanogaster.fasta"

handle = open(FASTA_file, "rU")
FASTA_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta"))
handle.close()

Exon_dict = {}
with open(Exon_file) as f:
    for line in f:
        Exon_dict[line.split()[-1].rsplit(";")[0].replace(":","")[3:]] = line.split()

previous = ""
CDS_sequence = ""
CDS_number_of_exons = int(0)
CDS_start = int(0)
CDS_end = int(0)
CDS_file = open (str("testfolder/" + "DUMMY" +"_CDS" + ".fas"),"w+")

for key, value in Exon_dict.iteritems():
    for k,v in FASTA_dict.iteritems():
        chromosome = value[0]
        direction = value[6]
        NMID = value[-1].split()[0].rsplit(';')[0][3:].rsplit("_")[0]
        if chromosome == k:
            start = value[3]
            end = value[4]
            exon_file = open (str("testfolder/" + key + ".fas"),"w+")
            exon_file.write(">" + key + " " + k + " " + start + " " + end + "\n")
            if direction == "+":
                exon_file.write(str(v.seq[int(int(start)-1):(int(end))]))
                CDS_sequence.__add__(str(v.seq[int(int(start)-1):(int(end))]))
            elif direction == "-":
                exon_file.write(str(v.seq[int(int(start)-1):(int(end))].reverse_complement()))
                CDS_sequence.__add__(str(v.seq[int(int(start)-1):(int(end))].reverse_complement()))

        if NMID != previous:
            CDS_file.write(">" + str(NMID) + " " + str(CDS_number_of_exons) + " " + str(chromosome) + " " + str(CDS_start) + " " + str(CDS_end) + " " + str(int(CDS_end)-int(CDS_start)) + "\n")
            CDS_file.write(CDS_sequence)
            CDS_file.close()
            CDS_file = open (str("testfolder/" + NMID +"_CDS" + ".fas"),"w+")
            previous = NMID
        else:
            CDS_number_of_exons += 1
            CDS_end += int(end)
            if start > CDS_start:
                pass
            else:
                CDS_start = start

Tags: 文件nameifstartfastafileintend

热门问题