生信算法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())

1
2
3
4
5
6
7
8
9

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

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 = True
    elif 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()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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()

1
2
3
4
5
6
7
8
9
10

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('>'): 
        # 拼接序列并去除末尾\n
        rna = rna + line.strip()

# 三个不同阅读框，转换为蛋白序列
for frame in range(3):
    # 0,1,2
    protein_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_seq
                protein_seq = protein_seq + codon_table[codon]
        else:
            # 处理非密码子表里的RNA序列，以-符号表示
            protein_seq = protein_seq + '-' 	

    
    # 每行48个氨基酸打印protein_seq
    i = 0
    while i < len(protein_seq):
        print(protein_seq[i:i + 48])
        i = i + 48
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

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] = seq
        seq = ''

    # >开头则获取蛋白序列编号， 否则添加氨基酸序列至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'])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

pdb文件SEQRES：

seq打印结果：

生信算法文章推荐

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

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

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

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

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

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