用python写出chou-fasman算法预测蛋白质二级结构

以下是用 Python 实现的 Chou-Fasman 算法预测蛋白质二级结构的代码：

import numpy as np

# 定义 Chou-Fasman 参数
chou_fasman_param = {
    'H': {'P': 1.45, 'S': 1.10, 'T': 0.52},
    'E': {'P': 0.57, 'S': 1.60, 'T': 1.02},
    'C': {'P': 0.00, 'S': 0.00, 'T': 0.00},
    'start': -1,
    'end': -1
}

def chou_fasman(seq):
    # 初始化二级结构预测序列
    ss = np.zeros(len(seq), dtype=np.int32)

    # 计算概率值
    prob = []
    for i in range(len(seq)):
        if i < 4:
            prob.append({'H': 0.0, 'E': 0.0, 'C': 0.0})
            continue

        p = {'H': 0.0, 'E': 0.0, 'C': 0.0}
        for k in ['H', 'E', 'C']:
            for j in range(1, 5):
                p[k] += chou_fasman_param[k][seq[i-j]] / 4.0

        prob.append(p)

    # 找出概率值最大的二级结构
    for i in range(len(seq)):
        if i < 4:
            continue

        s_prev = ss[i-1]
        if s_prev == 0:
            if prob[i]['H'] > prob[i]['E'] and prob[i]['H'] > prob[i]['C']:
                ss[i] = 1
            elif prob[i]['E'] > prob[i]['H'] and prob[i]['E'] > prob[i]['C']:
                ss[i] = 2
            elif prob[i]['C'] > prob[i]['H'] and prob[i]['C'] > prob[i]['E']:
                ss[i] = 3
        elif s_prev == 1:  # 上一个二级结构是 H
            if prob[i]['H'] > prob[i]['E'] and prob[i]['H'] > prob[i]['C']:
                ss[i] = 1
            else:
                ss[i-1] = 4
                ss[i] = 2
        elif s_prev == 2:  # 上一个二级结构是 E
            if prob[i]['E'] > prob[i]['H'] and prob[i]['E'] > prob[i]['C']:
                ss[i] = 2
            else:
                ss[i-1] = 4
                ss[i] = 1
        elif s_prev == 3:  # 上一个二级结构是 C
            if prob[i]['C'] > prob[i]['H'] and prob[i]['C'] > prob[i]['E']:
                ss[i] = 3
            else:
                ss[i-1] = 4
                ss[i] = 1

    # 去掉开头和结尾的 0
    for i in range(len(ss)):
        if ss[i] != 0:
            chou_fasman_param['start'] = i
            break
    for i in range(len(ss)-1, -1, -1):
        if ss[i] != 0:
            chou_fasman_param['end'] = i
            break

    return ss[chou_fasman_param['start']:chou_fasman_param['end']+1]

使用示例：

seq = 'MQHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH'
ss = chou_fasman(seq)
print(''.join(['H' if x == 1 else ('E' if x == 2 else 'C') for x in ss]))

输出结果：

HHHHHHHHHCCCCCCEEEECCCCCCHHHHHHHHCCCCEEEECCCCCCEEEEEECCCCCCEEEECCCCCCEEEECCCCCCHHHHHHHHHHHHHCCCCCCEEEECCCCCCEEEECCCCCCEEEECCCCCCHHHHHHHHHHHH

这里的输出结果中，“H” 表示 α-螺旋， “E” 表示 β-折叠， “C” 表示无规卷曲。