PSSM矩阵（位置权重矩阵）

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

PSSM矩阵（位置权重矩阵）

1. 定义

-参考维基百科 -同时参考这篇优秀的博客： www.nohup.cc/article/112…

可以反映出每个位置上不同碱基出现的频率，矩阵每一列表示相应位置上碱基出现的频率。 构造PSSM的第一步：通过计算每个位置上每个碱基出现的次数来创建一个基本频率矩阵（PFM) 第二步：标准化，用每个位置的原核苷酸计数除以序列数。构建位置频率矩阵 给定一个长度为l的序列集合X (N), 第三步：构建位置比重矩阵

2. 获得PSSM矩阵 参考了以下博客：

-www.nohup.cc/article/115… -blog.csdn.net/weixin_4340… -www.nohup.cc/article/110… -PSSM特征-从生成到处理 -如何生成PSSM矩阵 -blog.csdn.net/cpc78422148… -www.cnblogs.com/cong3Z/p/12…

• 第一步：下载blast，使用psiblast
• 第二步：下载数据库，下载swissprot，网址为：

ftp.ncbi.nlm.nih.gov/blast/db/或者… 进去后选择：swissprot， 下载链接为：ftp.ncbi.nlm.nih.gov/blast/db/sw… ftp://ftp.ncbi.nlm.nih.gov/blast/db/swissprot.tar.gz 注意下载到blast/bin路径下 解压缩：tar zxvf filename.tar

• 第三步：批量处理序列：

i = 0
fw = open('/blast/bin/0.txt', 'w')
for line in open('filename.fasta', 'r'):
    fw.write(line)
    i += 1
    if i % 2 == 0:
        fw.close()
        fw = open(str(i) + '.txt', 'w')
fw.close()

• 第四步：格式化数据库（可以不用做这一步）：

makeblastdb -in swissprot -parse_seqids -hash_index -dbtype prot(nulc)

• 第五步：批量对序列进行处理，获得其pssm矩阵：

import os
os.chdir('/blast/bin/sequence/')
for i in range(0,num(sequence),2):
    os.system("/blast/bin/psiblast -query /blast/bin/sequence/"+ str(i) + ".txt"+" -db /blast/swissprot -evalue 0.001 -num_iterations 3 -num_threads 38 -out /blast/bin/pssm_file/"+str(i)+".pssm")

注意-db前的空格，否则会报错“too many positional arguments” 第六步：处理得到的PSSM矩阵

pssm = []
with open('out.pssm') as f:
            lines = f.readlines()[3:-6]
            pssm = np.array([line.split()[2:22] for line in lines], dtype=int)```