生信算法7 - 核酸序列Fasta和蛋白PDB文件读写与检索

python 3.9实现以下算法。

1. 简单的写文件和读文件

# 写
file1 = open('count.txt','w')
file1.write('this is a test')
file1.close()# 读
file2 = open('my_file')
print(file2.read())

2. 将列表内容写入文本文件

# 生成100-500数字列表
data = [i * 100 for i in range(1, 6)]
print(data)context = []
for value in data:# 将内容追加至context列表context.append(str(value) + '\n')# 写入文件
open('results.txt', 'w').writelines(out)# 文件内容
# 100
# 200
# 300
# 400
# 500

3. 将NCBI Genbank文件转换为fasta文件

Genbank包含了所有已知的核酸和蛋白质序列，以及发表的期刊及生物学注释等信息。

AY810830.gb文件下载：
https://www.ncbi.nlm.nih.gov/nuccore/AY810830

genbank_file = open("AY810830.gb")
fasta_file = open("AY810830.fasta","w")flag = False
# 遍历文件每行
for line in genbank_file:# 写入ACCESSION编号if line.startswith('ACCESSION'):accession = line.split()[1].strip()fasta_file.write('>' + accession + '\n')# 存在ORIGIN，则存在fasta序列if line.startswith('ORIGIN'):flag = Trueelif flag:fields = line.split()if fields != []:print(seq)# fasta序列seq = ''.join(fields[1:])fasta_file.write(seq.upper() + '\n')genbank_file.close()
fasta_file.close()

4. 提取fasta序列header信息

fasta_file = open('AY810830.fasta','r')
out_file = open('AY810830.header','w')for line in fasta_file:# > 开头为fasta序列header信息if line.startswith('>'):out_file.write(line)out_file.close()

5. 转换RNA fasta序列为氨基酸序列

# 定义密码子表字典
codon_table = {'GCU':'A', 'GCC':'A', 'GCA':'A', 'GCG':'A', 'CGU':'R', 'CGC':'R',   'CGA':'R', 'CGG':'R', 'AGA':'R', 'AGG':'R', 'UCU':'S', 'UCC':'S','UCA':'S', 'UCG':'S', 'AGU':'S', 'AGC':'S', 'AUU':'I', 'AUC':'I','AUA':'I', 'UUA':'L', 'UUG':'L', 'CUU':'L', 'CUC':'L', 'CUA':'L','CUG':'L', 'GGU':'G', 'GGC':'G', 'GGA':'G', 'GGG':'G', 'GUU':'V','GUC':'V', 'GUA':'V', 'GUG':'V', 'ACU':'T', 'ACC':'T', 'ACA':'T','ACG':'T', 'CCU':'P', 'CCC':'P', 'CCA':'P', 'CCG':'P', 'AAU':'N','AAC':'N', 'GAU':'D', 'GAC':'D', 'UGU':'C', 'UGC':'C', 'CAA':'Q','CAG':'Q', 'GAA':'E', 'GAG':'E', 'CAU':'H', 'CAC':'H', 'AAA':'K','AAG':'K', 'UUU':'F', 'UUC':'F', 'UAU':'Y', 'UAC':'Y', 'AUG':'M','UGG':'W','UAG':'STOP', 'UGA':'STOP', 'UAA':'STOP'}# 读取RNA fasta文件
rna = ''
for line in open('A06662-RNA.fasta'):# 过滤>开头行if not line.startswith('>'): # 拼接序列并去除末尾\nrna = rna + line.strip()# 三个不同阅读框，转换为蛋白序列
for frame in range(3):# 0,1,2protein_seq = '' print('Reading frame ' + str(frame + 1))for i in range(frame, len(rna), 3):codon = rna[i:(i + 3)]if codon in codon_table:# 判断是否为终止密码子if codon_table[codon] == 'STOP':# *符号表示终止密码子protein_seq = protein_seq + '*'else: # 不是终止密码子则添加转换后的氨基酸名称至protein_seqprotein_seq = protein_seq + codon_table[codon]else:# 处理非密码子表里的RNA序列，以-符号表示protein_seq = protein_seq + '-' 	# 每行48个氨基酸打印protein_seqi = 0while i < len(protein_seq):print(protein_seq[i:i + 48])i = i + 48

6. 将fasta序列转换为字典

P62258氨基酸序列下载：
https://www.ncbi.nlm.nih.gov/protein/P62258

sequences = {}
ac = ''
seq = ''# 遍历fasta文件
for line in open("P62258.fasta"):# header信息保存至字典if line.startswith('>') and seq != '':sequences[ac] = seqseq = ''# >开头则获取蛋白序列编号， 否则添加氨基酸序列至seq变量if line.startswith('>'):ac = line.split('|')[1]else:seq = seq + line.strip()sequences[ac] = seq# 打印全部key
print(sequences.keys())# 打印指定Key字典氨基酸序列
print(sequences['P62258.1'])

7. 从pdb文件提取氨基酸序列

PDB数据库是一个数据中心，主要包含：原子坐标，蛋白质结构的其他信息和除蛋白以外生物大分子的信息。pdb文件可以从该数据下载。

牛β-胰蛋白酶 pdb文件下载:
https://www.rcsb.org/structure/1TLD

# 氨基酸简写字典
aa_codes = {'ALA':'A', 'CYS':'C', 'ASP':'D', 'GLU':'E','PHE':'F', 'GLY':'G', 'HIS':'H', 'LYS':'K','ILE':'I', 'LEU':'L', 'MET':'M', 'ASN':'N','PRO':'P', 'GLN':'Q', 'ARG':'R', 'SER':'S','THR':'T', 'VAL':'V', 'TYR':'Y', 'TRP':'W'}seq = ''
# 遍历.pdb文件
for line in open("1tld.pdb"):# SEQRES开头行为氨基酸序列if line.startswith("SEQRES"):line_split = line.split()print(line_split)# 拼接氨基酸序列for aa_code in line_split[4:]:seq = seq + aa_codes[aa_code]# 打印拼接结果， 每行63个氨基酸
i = 0
print(">1tld")
while i < len(seq):print(seq[i:(i + 64)])i = i + 64

pdb文件SEQRES：

seq打印结果：

生信算法文章推荐

生信算法1 - DNA测序算法实践之序列操作

生信算法2 - DNA测序算法实践之序列统计

生信算法3 - 基于k-mer算法获取序列比对索引

生信算法4 - 获取overlap序列索引和序列的算法

生信算法5 - 序列比对之全局比对算法

生信算法6 - 比对reads碱基数量统计及百分比统计