提取大序列文件的子集

最新推荐文章于 2021-05-14 02:00:11 发布

原创最新推荐文章于 2021-05-14 02:00:11 发布 · 448 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#python

python 专栏收录该内容

6 篇文章

订阅专栏

Python3.8

Conda

Python

Python 是一种高级、解释型、通用的编程语言，以其简洁易读的语法而闻名，适用于广泛的应用，包括Web开发、数据分析、人工智能和自动化脚本

有时候需要比较大的序列文件里提取出它的子集，通常是根据序列号来操作的，上代码：

#存入字典
def read_to_dict(in_file):
	sequence = {}
	ac = ''
	seq = ''
	for line in open(in_file):
		if line.startswith('>') and seq != '':
			sequence[ac] = seq
			seq = ''
		if line.startswith('>'):
			ac = line.strip().split()[1] #处理序列key命名，视情况而定
		else:
			seq = seq + line.strip()
	sequence[ac] = seq
	return sequence

#待提取序列号读取入列表
def read_to_list(in_file2):
	list2 = []
	with open(in_file2, 'r', encoding='utf-8') as list1:
		for line in list1:
			list2.append(line.strip())
	return list2

#提取
other_features = read_to_dict('other_genomic.fasta')
list = read_to_list('search_results.txt')
handle = open("selected.txt", "w")
for acc in list:
	if acc in other_features.keys():
		handle.write('>' + acc + '\n' + other_features[acc] + '\n')
	else:
		print(acc)
handel.close()

原文件：

>ABC102 ABC102 SPXID:S000121252, Chr I from 608-776, Genome Release 55-0-1, ABC, “Autonomously Replicating Sequence”
TATACACTTATGTCAATATTACAGAAAAATCCCCACAAAAATCACCTAAACATAAAAATATTCTACTTTT
…

待提取序列：

ABC301
ABC305
ABC306
ABC901
…

输出结果：

>ABC301
TATACACTTATGTCAATATTACAGAAAAATCCCCACAAAAATCACCTAAACATAAAAATATTCTACTTTT
…

BioPython里也有类似的操作，但是输出效果不好，不如上面的灵活。

>>> from Bio import SeqIO
>>> uniprot = SeqIO.index("other_genomic.fasta", "fasta")
>>> handle = open("selected.txt", "w")
>>> for acc in ["ABC301", "ABC305", "ABC306", "ABC901", ... "ABC998"]:
...     handle.write(uniprot.get_raw(acc))
>>> handle.close()