以下代碼均為個人原創，如有疑問，歡迎交流。新浪微博：拾毅者
本節內容：
- pssm生成
- pssm簡化
- 標准的pssm構建
- 滑動pssm生成
在基於蛋白質序列的相關預測中。使用PSSM打分矩陣會得將預測效果大大提高，同一時候，假設使用滑動的PSSM，效果又會進一步提高。這里主要以分享代碼為主，以下介紹下PSSM從生成到處理的全過程。
1.PSSM的生成
PSSM的生成有多種方式，這里使用的psiblast軟件。ncbi-blast-2.2.28+/bin/psiblast。下載地址：http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastNews#1 用法。輸入一個序列，加上相關參數，直接輸出PSSM文件
代碼
#一個命令函數，依據pdb文件。使用blast生成pssm文件
def command_pssm(content, output_file,pssm_file):
os.system('/ifs/share/lib/blast/ncbi-blast-2.2.28+/bin/psiblast \ -query %s \ -db /ifs/data/database/blast_data/nr \ -num_iterations 3 \ -out %s \ -out_ascii_pssm %s &' %(content, output_file,pssm_file))上面是運行的命令，封裝成函數，以下為實際代碼：
#對每一個序列進行PSSM打分
def pssm(proseq,outdir):
inputfile = open(proseq,'r')
content = ''
input_file = ''
output_file = ''
pssm_file = ''
chain_name = []
for eachline in inputfile:
if '>' in eachline:
if len(content):
temp_file = open(outdir + '/fasta/' + chain_name,'w')
temp_file.write(content)
input_file = outdir + '/fasta/' + chain_name
output_file = outdir + '/' + chain_name + '.out'
pssm_file = outdir + '/' + chain_name + '.pssm'
command_pssm(input_file, output_file,pssm_file)
temp_file.close
content = ''
chain_name = eachline[1:5] + eachline[6:7]
content += ''.join(eachline)
#print content
#print chain_name
if len(content):
temp_file = open(outdir + '/fasta/' + chain_name,'w')
temp_file.write(content)
input_file = outdir + '/fasta/' + chain_name
output_file = outdir + '/' + chain_name + '.out'
pssm_file = outdir + '/' + chain_name + '.pssm'
command_pssm(input_file, output_file,pssm_file)
temp_file.close
inputfile.close() 測試用例：
''' #生成pssm文件，迭代次數為3 proseq = '/ifs/home/liudiwei/experiment/step2/data/protein.seq' outdir = '/ifs/home/liudiwei/experiment/step2/pssm' pssm(proseq,outdir) '''PSSM輸出例子：
2.簡化PSSM數據
通常我們須要的僅僅是前面的20列
以下通過代碼來實現上面的功能：
#格式化pssm每行數據
def formateachline(eachline):
col = eachline[0:5].strip()
col += '\t' + eachline[5:8].strip()
begin = 9
end = begin +3
for i in range(20):
begin = begin
end = begin + 3
col += '\t' + eachline[begin:end].strip()
begin = end
col += '\n'
return col簡化pssm。僅僅要得到前面的20個氨基酸的打分值
def simplifypssm(pssmdir,newdir):
listfile = os.listdir(pssmdir)
for eachfile in listfile:
with open(pssmdir + '/' + eachfile,'r') as inputpssm:
with open(newdir + '/' + eachfile,'w') as outfile:
count = 0
for eachline in inputpssm:
count +=1
if count <= 3:
continue
if not len(eachline.strip()):
break
oneline = formateachline(eachline)
outfile.write(''.join(oneline))
''' Test example pssmdir = '/ifs/home/liudiwei/experiment/step2/pssm/oldpssm' newdir = '/ifs/home/liudiwei/experiment/step2/pssm/newpssm' simplifypssm(pssmdir, newdir) '''3.得到標准的PSSM
通過上面抽取出來的PSSM，以下通過代碼來獲得一個滑動的PSSM
#標准的pssm,直接依據標准的pssm滑動
def standardPSSM(window_size,pssmdir,outdir):
listfile = os.listdir(pssmdir)
for eachfile in listfile:
outfile = open(outdir + '/' + eachfile, 'w')
with open(pssmdir + '/' + eachfile, 'r') as inputf:
inputfile = inputf.readlines()
for linenum in range(len(inputfile)):
content = []
first = [];second = [];third=[];last=[]
if linenum < window_size/2:
for i in range((window_size/2 - linenum)*20):
second.append('\t0')
if window_size/2 - linenum > 0:
countline = window_size - (window_size/2 - linenum)
else:
countline = window_size #get needed line count
linetemp = 0
for eachline in inputfile:
if linetemp < linenum-window_size/2:
linetemp += 1
continue
if linetemp == linenum:
thisline = eachline.split('\t')
for j in range(0,2):
if j>0:
first.append('\t')
first.append(thisline[j].strip())
if countline > 0:
oneline = eachline.split('\t')
for j in range(2,len(oneline)):
third.append('\t' + oneline[j].strip())
countline -=1
else:
break
linetemp += 1
while countline:
for i in range(20):
last.append('\t0')
countline -=1
content += first + second + third + last
outfile.write(''.join(content) + '\n')
outfile.close()
'''Test example pssmdir = '/ifs/home/liudiwei/experiment/step2/pssm/newpssm' newdir = '/ifs/home/liudiwei/experiment/step2/pssm/standardpssm' window_size = 5 standardPSSM(window_size,pssmdir, newdir) '''4.依據滑動窗體求出滑動的PSSM
#依據窗體大小,計算出滑動后的20個氨基酸打分值
def computedPSSM(window_size,pssmdir,outdir):
listfile = os.listdir(pssmdir)
for eachfile in listfile:
outfile = open(outdir + '/' + eachfile, 'w')
with open(pssmdir + '/' + eachfile, 'r') as inputf:
inputfile = inputf.readlines()
for linenum in range(len(inputfile)):
content = []
first = [];second = []
if window_size/2 - linenum > 0:
countline = window_size - (window_size/2 - linenum)
else:
countline = window_size #get needed line count
linetemp = 0
for eachline in inputfile:
if linetemp < linenum-window_size/2:
linetemp += 1
continue
if linetemp == linenum:
thisline = eachline.split('\t')
for j in range(0,2):
if j>0:first.append('\t')
first.append(thisline[j].strip())
if countline > 0:
oneline = eachline.split('\t')[2:len(eachline)]
tline = []
for i in range(len(oneline)):
tline.append(int(oneline[i]))
if len(second)==0:
second += tline
else:
second = list(map(lambda x: x[0]+x[1], zip(second, tline)))
countline -=1
else:
break
linetemp += 1
format_second = []
for i in range(len(second)):
format_second.append('\t' + str(second[i]))
content += first + format_second
outfile.write(''.join(content) + '\n')
outfile.close()
''' pssmdir = '/ifs/home/liudiwei/experiment/step2/pssm/newpssm' newdir = '/ifs/home/liudiwei/experiment/step2/pssm/computedpssm' window_size = 5 computedPSSM(window_size,pssmdir, newdir) '''平滑的PSSM，僅僅是pssmdir不同，直接調用standardPSSM函數
def smoothedPSSM(window_size,pssmdir,outdir):
standardPSSM(window_size,pssmdir, outdir)
'''Test example pssmdir = '/ifs/home/liudiwei/experiment/step2/pssm/computedpssm' newdir = '/ifs/home/liudiwei/experiment/step2/pssm/smoothedpssm' window_size = 5 smoothedPSSM(window_size,pssmdir,newdir) '''最后得到的是一個滑動的PSSM矩陣，特征的維數隨窗體的大小逐漸增減。