使用python折叠exons
使用python折叠exons
提问于 2017-03-24 20:33:52
我正在尝试根据包含外显子所属基因信息的另一个文件，从一个文件中折叠外显子。
其中一个文件包含外显子的名称以及它们对应的基因。
举个例子：
Exon1 Gene4
Exon2 Gene5
Exon3 Gene8
Exon4 Gene8另一个文件是包含外显子名称和序列的快速文件，例如：
>Exon1
ACGTCGATTTCGATCA
>Exon2
ACGTCGATTTCGATCAACGTAA
>Exon3
ACGTCGATTTCAGCGGATCA
>Exon4
CCCACGTCGATTTCGATCAACGTAA因此，我的目标是输出一个包含基因名称和相应外显子折叠成一个序列(CDS)的FASTA file。举个例子：
>Gene4
ACGTCGATTTCGATCA
>Gene5
ACGTCGATTTCGATCAACGTAA
>Gene8
ACGTCGATTTCAGCGGATCACCCACGTCGATTTCGATCAACGTAA到目前为止，我已经创建了一个字典，将基因名称作为key，外显子名称作为key列表，如下所示：
exon_to_gene_dict = {}
with open('gene_to_exon', 'r') as file_gene, open('exon.fasta', 'r') as file_exons:
for line in file_gene:
line = line.rstrip().split()
scaffold = line[0]
gene = line[1]
exon_to_gene_dict.setdefault(gene, []).append(scaffold)然后我创建了一个包含外显子id和相应序列的字典：
exon_sequence_dict = {}
with open(exons.fasta,'r') as file_exons:
for line in file_exons:
line = line.rstrip()
if line.startswith('>'):
exon_id = line[1:]
else:
sequence = line
exon_sequence_dict[exon_id] = sequence
gene_to_cds_dict = {}
for exon, seq in exon_sequence_dict.items():
for gene_id, exon1 in gene_to_exon_dict.items():
if exon in exon1:
gene_to_cds_dict.setdefault(gene_id, []).append(seq)
gene_id_to_sequence =[]
for gene_id, sequence in sorted(gene_to_cds_dict.items()):
print(gene_id, file=f_out)
print(''.join(sequence), file=f_out)所以最后我还是想出了上面的方法
该示例的输出为：
Gene4
ACGTCGATTTCGATCA
Gene5
ACGTCGATTTCGATCAACGTAA
Gene8
ACGTCGATTTCAGCGGATCACCCACGTCGATTTCGATCAACGTAA谢谢你的帮助
回答 2
Stack Overflow用户
发布于 2017-03-24 21:43:30
这是一个完整的解决方案，使用BioPython。BioPython会自动处理the issue that Cory Weller raised (如果您的FASTA将序列拆分到多行)。这个解决方案还将来自同一基因的序列连接成CDS。
from Bio import SeqIO
from Bio import Seq
from collections import defaultdict
with open('gene_to_exon.txt') as f:
lines = f.read().splitlines()
exon_to_gene = dict(tuple(line.split()) for line in lines)
exon_records = SeqIO.to_dict(SeqIO.parse('exon.fasta', 'fasta'))
gene_records = defaultdict(list)
for exon_id, record in exon_records.items():
gene_id = exon_to_gene[exon_id]
gene_records[gene_id].append(record)
cds_records = []
for gene_id in gene_records:
gene_record = gene_records[gene_id][0]
if len(gene_records[gene_id]) > 1:
sequence = ''.join([str(gene_record.seq) for gene_record in gene_records[gene_id]])
gene_record.seq = Seq.Seq(sequence)
gene_record.id = gene_id
gene_record.description = gene_id # see https://www.biostars.org/p/156428/
cds_records.append(gene_record)
with open('output.fasta', 'w') as g:
SeqIO.write(cds_records, g, 'fasta')在您的输入文件上，这将给出结果：
>Gene8
ACGTCGATTTCAGCGGATCACCCACGTCGATTTCGATCAACGTAA
>Gene5
ACGTCGATTTCGATCAACGTAA
>Gene4
ACGTCGATTTCGATCAPython :这里的也是一个纯实现。FASTA解析代码来自this code monk blog post (适用于Python3)，作者在其中详细描述了使用itertools.groupby相对于使用for循环的优势。
import itertools
from collections import defaultdict
def isheader(line):
return line.startswith('>')
def aspairs(filehandle):
for header,group in itertools.groupby(filehandle, isheader):
if header:
line = next(group)
ensembl_id = line[1:].split()[0]
else:
sequence = ''.join(line.strip() for line in group)
yield ensembl_id, sequence
with open('exon.fasta') as f:
exon_records = dict(aspairs(f))
with open('gene_to_exon.txt') as f:
lines = f.read().splitlines()
exon_to_gene = dict(tuple(line.split()) for line in lines)
gene_records = defaultdict(str)
for exon_id, seq in exon_records.items():
gene_id = exon_to_gene[exon_id]
gene_records[gene_id] += seq
with open('output.fasta', 'w') as f:
for gene_id, seq in gene_records.items():
f.write('>{}\n{}\n'.format(gene_id, seq))Stack Overflow用户
发布于 2017-03-24 21:10:02
您已经完成了第一部分，现在需要创建一个将外显子链接到序列的字典。你怎么做取决于你的FASTA文件的结构。如果您的FASTA文件在多行上有序列(例如，最大行长为40个字符)，则需要在序列字典中找到将它们连接在一起的方法。否则，如果所有的核苷酸都在一行上，你会更轻松。你应该能够先用">“拆分，然后再用换行符拆分。
with open('exon.fasta','r') as infile:
exon_text = infile.read()
exon_text = exon_text.split(">")[1:]
exon_text = [x.split("\n")[0:2] for x in exon_text]
exonDict = {}
for exon in text:
exonDict[exon[0]] = exon[1]现在您将拥有第二个字典，外显子名称作为关键字，序列作为条目。有了这两个部分，您应该能够遍历和连接所有序列。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/42999643
复制相关文章
相似问题
腾讯云开发者
