蛋白质序列特征提取方法之——PSSM

Overview

我在之前写的一篇博客中谈到整理那些混乱的数据源，发现有pssm fts文件夹中的子文件夹和文件并不清楚来龙去脉，这个问题困扰了我一段时间。最近在研究PSSM算法时，与Chris交流了一下，恍然大悟：这个文件夹中的t3pssm,t4pssm,t6pssm三个子文件夹中的形如t6_12.pssm的文件族，是由t3,t4,t6这三个文件夹中的形如t6_12.fasta的文件族经过BLAST或者PSI-BLAST等算法程序处理过后得到的数据。并且t3pssm,t4pssm,t6pssm三个子文件夹中的每个子文件中包含有两个矩阵，其中第一个即为我们在构建PSSM的步骤中提到的PSSM。这些PSSM将会在下一步特征提取的过程中做为比对(alignment)的模板。
由于并未拿到BLAST和PSI-BLAST程序的源码，所以这里只介绍从PSSM开始到提取特征的步骤。

1.`Python`编程方案

我们的项目在提取特征这一阶段的第六个步骤中，运用的就是这个PSSM方法，这一步只有一个Linux命令

python t34pssm.py T4undrsmp.txt ./t4 ./t4pssm

其中，t34pssm这个Python文件调用了另一个Python文件pssm_ft.py，这两个文件是用来计算提取特征的；T4undrsmp.txt是输入文件，即被处理数据源；./t4和./t4pssm是两个文件夹。
下面我们详细解析这个Python编程方案的每一步的代码和作用。在解析每一步的输出时，我在原代码的基础上引入了一个time包

import time

这个包提供给我们一个很好的函数sleep(),这个函数输入的数基本时间单位为秒。这样我们在遇到大量循环的时候，可以在循环体内加上这么一句：

time.sleep(1)
print #所需要的打印内容

就可以观察每一步循环的输出内容了，睡眠时间可以根据需要调整，十分方便。在终止循环时，可以在命令行界面按下ctrl+c停止进程，这样就能记录下输出内容了。

1.1 解析`t34pssm.py`文件

引入相关的包:

import fileinput
import sys
from os import listdir
from os.path import isfile, join
import re
from pssm_ft import *　#引入pssm_ft.py文件
import time

读取输入参数[1]和[2]：

smplfasta = sys.argv[1]           #T4undrsmp.txt
spfasta = sys.argv[2]             #./t4
check_head = re.compile(r'\>')    #正则表达式匹配符号>

有一点需要说明:由于T4undrsmp.txt这个文件中原有的数据量太大，故输出量太大，耗时太大，我们仅取其第一条蛋白质序列作为例子说明这个算法的运算过程。第一条蛋白质序列如下：

>lpn:lpg0012 hypothetical protein (A)|1
MNTTEHTELGNGLLMDKIEGNPYLRIDEFGVLHLRLQRFGENNLPEPMELEMSAGEIIAMAGDYFTQANW
TMDLDLPKCELFNSPAELGKHLIRKPIDPKEENALITAYNNLAAPDVTRKEIDRIYSINNANYVPFSPTL
NFYAQQLMYYFRVKDYGEMLVRNQTHFTPWSIRVYILGHAIALRYARLSYELKQLATDRNYQSDNPDLST
LKISLQNKNETLSSNTLLDLANRYHAQAYSIELFTFHYYSDHFATGHMSMIGDLRVVLKERFGVWGNILA
NNLHDEVNRVGVYTVRPYDPTPNTTEAPSRARGDGKFDTCLNQFNRLACLNGMTASLKDINQVLSGGVIP
EQKKFGGLVHMPDVDFNSRQHQPLLVLSKGKVYYRNNLSQIHIISPSEYEALRENPEKHGYKELTSKWAA
FKLVTKLRLFPYLYDGSVMPVSDEKLAEIIADEKRRNPQRAPIPTPSCLPESEPTVFDWRTKASWRNNKD
SLDILDGLKKHSILAAKQTHSAIQEEEEVRLNLGV

从T4undrsmp.txt文件中读取undersample fasta格式的蛋白质序列：

smplist = []              #存储蛋白质序列的名称序号信息
smpcnt = 0                #记录读取的蛋白质数量
for line, strin in enumerate(fileinput.input(smplfasta)):
if check_head.match(strin):
    smplist.append(strin.strip())
    #print smplist        #['>lpn:lpg0012 hypothetical protein (A)|1']
    #time.sleep(1)        #休眠一秒以观察输出
    smpcnt += 1

读取1232个fasta格式的文件：

onlyfiles = [ f for f in listdir(spfasta) if isfile(join(spfasta,f)) ]  
#print onlyfiles    
#输出t4_1.fasta等1232个文件的列表,['t4_961.fasta', 't4_456.fasta', 't4_1179.fasta',……, 't4_954.fasta', 't4_567.fasta']
onlyfiles.remove('.DS_Store')        #移除干扰项

将蛋白质名和文件名对应起来：

fastaDict = {}                          #存储“蛋白质名：文件名”的键值对

for fi in onlyfiles:
    cntnt = ''
    for line, strin in enumerate(fileinput.input(spfasta+'/'+fi)):  
        if line == 0:
            cntnt += strin.strip()      #去除前后空格
    if cntnt in fastaDict:              #检查否重复(已经存在)
        print strin                     #打印重复(已经存在)
    fastaDict[cntnt] = fi               #生成键值对
print fastaDict     #{'>lpn:lpg0012 hypothetical protein (A)|1': 't4_1.fasta',……,'>gi|148360344|ref|YP_001251551.1| hypothetical protein LPC_2282 [Legionella pneumophila str. Corby]|-1': 't4_580.fasta'} 一共1232个键值对

读取输入参数[3]：

pssmdir = sys.argv[3]             #PSSM模板文件目录
finalist = []                     #存储PSSM文件名字的列表 
for smp in smplist:
    #print fastaDict[smp]       　#t4_1.fasta
    finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm')       #根据fasta文件的名字找出对应的PSSM文件，输出t4_1.pssm
    #print finalist               #['./t4pssm/t4_1.pssm']

for fi in finalist:            　  #读取PSSM文件列表
    pssm(fi, 't4pssm.txt')        #调用函数文件pssm_ft.py中的pssm(fi,output)函数

调用函数pssm(fi,output)这一步表明真正计算的步骤将要开始，且是t4pssm.txt和PSSM之间的运算，与其他文件无关，更印证了我们最初的想法。

1.2 解析`pssm_ft.py`文件

这个pssm_ft.py文件的作用就是Retrieve PSSM features，也就是检索PSSM特征。将T4undrsmp.txt这个文件中的序列跟PSSM进行比对，得出蛋白质序列的特征向量。
引入相关的包：

import sys
import numpy as np
import math
import re
import fileinput 
import time

numpy这个包是用来处理大型矩阵运算的。
既然pssm(fi,output)这个函数首先被调用，我们就先分析一下它。

def pssm(fi,output):    #pssm(fi, 't4pssm.txt')     fi==./t4pssm/t4_1.pssm
    Amino_vec = "ARNDCQEGHILKMFPSTWYV"    #20个基本氨基酸
    
    PSSM = []            #存储PSSM的特征向量矩阵
    ACC = [ [0.0] * 20 ] * 20
    seq_cn = 0
    for line, strin in enumerate(fileinput.input(fi)):
        if line > 2:
        str_vec = strin.split()[1:22]    #切片截取该行的第１个至第21个字符存入向量str_vec中
        #print str_vec   #['M', '-2', '-3', '-4', '-5', '-3', '-2', '-4', '-4', '-3', '0', '3', '-3', '8', '-1', '-4', '-3', '-2', '-3', '-2', '0']等525个向量
        if len(str_vec) == 0:
            break
        PSSM.append(map(int, str_vec[1:]))    #将每个向量中的氨基酸符号去掉，其他字符转换为int类型，并全部存入PSSM特征向量矩阵    
        seq_cn += 1        #525
        ACC[Amino_vec.index(str_vec[0])] = map(sum, zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])]))    #这句最长，从赋值号(=)右边开始逐步分析

"""
1. Amino_vec.index(str_vec[0])　返回本行氨基酸符号在Amino_vec中的位置
2. ACC[Amino_vec.index(str_vec[0])] 初始值均为20维向量[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
3. map(int, str_vec[1:])　将每个向量中的氨基酸符号去掉，其他字符转换为int类型
4. zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])])　返回一个tuple列表，形如[(-2,0.0),(-3,0.0),……,(-2,0.0),(0,0.0)]
5. map(sum, zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])]))　将tuple列表中各元素转换为sum类型,即求和
6. 循环迭代，最后ACC(20维)为：
[[85.0, -27.0, -42.0, -58.0, -59.0, -14.0, -15.0, -25.0, -46.0, -53.0, -55.0, -20.0, -32.0, -75.0, -32.0, 10.0, -17.0, -88.0, -57.0, -28.0], [-23.0, 90.0, -3.0, -15.0, -80.0, 24.0, -3.0, -52.0, -10.0, -45.0, -42.0, 19.0, -31.0, -63.0, -44.0, -17.0, -10.0, -68.0, -32.0, -48.0], ……　, [-5.0, -33.0, -36.0, -34.0, -36.0, -30.0, -26.0, -58.0, -49.0, 23.0, 11.0, -19.0, -12.0, -36.0, -36.0, -31.0, -5.0, -72.0, -35.0, 42.0]]

"""
    fileinput.close()
    
    ACC_np = np.array(ACC)    #转换成数组
    ACC_np = np.divide(ACC_np, seq_cn)    #将每一个元素除以525
    amcnt_max = np.amax(ACC_np)    #找出最大的元素0.28380952381
    amcnt_min = np.amin(ACC_np)    ##找出最小的元素-0.274285714286
    
    vfunc = np.vectorize(scale) #矢量化scale函数
    """
    numpy中的vectorize函数可以将scale函数转换成可以接受向量参数的函数
    """
    ACC_np = vfunc(ACC_np, amcnt_max,amcnt_min)    #调用scale函数，将矩阵范围控制在[-1,1]之间
    ACC_shp = np.shape(ACC_np)    #(20,20)
    PSSM = np.array(PSSM)
    PSSM_AC = pssm_ac_cal(PSSM, GROUP, seq_cn)      #调用函数pssm_ac_cal(PSSM, g, l)   
"""
pssm_ac_cal(PSSM, g, l)函数，三个参数：1.PSSM矩阵；2.10;3.蛋白质序列长度
返回10×20的矩阵PSSM_AC
"""               
    PSSM_shp = np.shape(PSSM_AC)    #(10, 20)
    file_out = file(output,'a')    #打开文件t4pssm.txt，如果没有则创建一个；ａ表示add,以追加的方式打开
    np.savetxt(file_out, [np.concatenate((np.reshape(ACC_np, (ACC_shp[0] * ACC_shp[1], )), np.reshape(PSSM_AC, (PSSM_shp[0] * PSSM_shp[1], ))))], delimiter=",")    #这句虽然比较长，但是逻辑比较简单清晰
"""
1. reshape函数将矩阵的行、列、维重新调整为200×1的矩阵.
2. concatenate函数将两个矩阵连接起来成为一个新矩阵.
3. savetxt函数将矩阵保存为txt格式文件，分隔符为",".
"""

现在看一下另外两个函数scale和pssm_ac_cal。
scale函数:

def scale(x, max, min):
"""
把矩阵中最大的数和最小的数变为1和-1,其他数按以下比例缩放
"""
    try:
        return -(1-(x-min)/(max-min))+((x-min)/(max-min))
    except Exception:
        print "ZeroDivisionError, max-min:\s", max-min

pssm_ac_cal函数:

def pssm_ac_cal(PSSM, g, l):
    """
    输入PSSM、GROUP、蛋白质序列的长度，返回GROUP×20的PSSM_AC矩阵
    """
    PSSM_AC = np.array([ [0.0] * 20 ] * g)

    for pg in range(g):
        l_g = l - pg - 1
        for pj in range(20):
            sum_jl = 0.0
            for i in range(l):
                sum_jl += PSSM[i][pj]
            sum_jl /= l
        
            pssm_acjg = 0.0
            for i in range(l_g):
                pssm_acjg += (PSSM[i][pj]-sum_jl) * (PSSM[i+pg+1][pj]-sum_jl)
        pssm_acjg /= l_g
            PSSM_AC[pg][pj] = pssm_acjg
    return PSSM_AC

至此，Python编程方案已解析完毕。疏漏之处请与我们联系，一定尽快修改。