以下代碼均為個(gè)人原創(chuàng)，如有疑問，歡迎交流。新浪微博：拾毅者

本節(jié)內(nèi)容：

• pssm生成
• pssm簡化
• 標(biāo)準(zhǔn)的pssm構(gòu)建
• 滑動pssm生成

在基于蛋白質(zhì)序列的相關(guān)預(yù)測中。使用PSSM打分矩陣會得將預(yù)測效果大大提高，同一時(shí)候，假設(shè)使用滑動的PSSM，效果又會進(jìn)一步提高。這里主要以分享代碼為主，以下介紹下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 用法。輸入一個(gè)序列，加上相關(guān)參數(shù)，直接輸出PSSM文件

代碼

#一個(gè)命令函數(shù)，依據(jù)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))

上面是運(yùn)行的命令，封裝成函數(shù)，以下為實(shí)際代碼：

#對每一個(gè)序列進(jìn)行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_nameoutput_file = outdir + '/' + chain_name + '.out'pssm_file = outdir + '/' + chain_name + '.pssm' command_pssm(input_file, output_file,pssm_file)temp_file.closecontent = ''chain_name = eachline[1:5] + eachline[6:7]content += ''.join(eachline)#print content#print chain_nameif len(content):temp_file = open(outdir + '/fasta/' + chain_name,'w')temp_file.write(content)input_file = outdir + '/fasta/' + chain_nameoutput_file = outdir + '/' + chain_name + '.out'pssm_file = outdir + '/' + chain_name + '.pssm'command_pssm(input_file, output_file,pssm_file) temp_file.closeinputfile.close()

測試用例：

''' #生成pssm文件，迭代次數(shù)為3 proseq = '/ifs/home/liudiwei/experiment/step2/data/protein.seq' outdir = '/ifs/home/liudiwei/experiment/step2/pssm' pssm(proseq,outdir) '''

PSSM輸出例子：

2.簡化PSSM數(shù)據(jù)

通常我們須要的僅僅是前面的20列

以下通過代碼來實(shí)現(xiàn)上面的功能：

#格式化pssm每行數(shù)據(jù) def formateachline(eachline):col = eachline[0:5].strip()col += '\t' + eachline[5:8].strip() begin = 9 end = begin +3for i in range(20): begin = beginend = begin + 3col += '\t' + eachline[begin:end].strip()begin = endcol += '\n'return col

簡化pssm。僅僅要得到前面的20個(gè)氨基酸的打分值

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 = 0for eachline in inputpssm:count +=1if count <= 3:continueif not len(eachline.strip()):breakoneline = 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.得到標(biāo)準(zhǔn)的PSSM

通過上面抽取出來的PSSM，以下通過代碼來獲得一個(gè)滑動的PSSM

#標(biāo)準(zhǔn)的pssm,直接依據(jù)標(biāo)準(zhǔn)的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 countlinetemp = 0for eachline in inputfile:if linetemp < linenum-window_size/2:linetemp += 1continueif 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 -=1else:breaklinetemp += 1while countline:for i in range(20):last.append('\t0')countline -=1content += first + second + third + lastoutfile.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.依據(jù)滑動窗體求出滑動的PSSM

#依據(jù)窗體大小,計(jì)算出滑動后的20個(gè)氨基酸打分值 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 countlinetemp = 0for eachline in inputfile:if linetemp < linenum-window_size/2:linetemp += 1continueif 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 += tlineelse: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不同，直接調(diào)用standardPSSM函數(shù)

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

最后得到的是一個(gè)滑動的PSSM矩陣，特征的維數(shù)隨窗體的大小逐漸增減。

---

本欄目python代碼分享持續(xù)更新中，歡迎關(guān)注：Dream_Angel_Z博客 新浪微博：拾毅者

總結(jié)

以上是生活随笔為你收集整理的PSSM特征-从生成到处理的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。