快捷導(dǎo)航

python?服務(wù)器批處理得到PSSM矩陣的問(wèn)題

更新時(shí)間：2022年07月21日 11:37:37 作者：李劃水員

這篇文章主要介紹了python?服務(wù)器批處理得到PSSM矩陣,本文給大家介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或工作具有一定的參考借鑒價(jià)值，需要的朋友可以參考下

1. 在linux上安裝psiblast

最好新建一個(gè)python環(huán)境，因?yàn)槲野l(fā)現(xiàn)conda安裝blast默認(rèn)的是python==3.6.11，可能會(huì)不小心把你的python版本改掉…然后你寫(xiě)好的代碼全die了……

conda create -n blast python==3.6.11
source activate blast
conda install -c bioconda blast

2.下載并編譯用于比對(duì)的大型蛋白質(zhì)數(shù)據(jù)庫(kù)

nr和uniprot是比較通用的數(shù)據(jù)庫(kù)：

ftp://ftp.ncbi.nlm.nih.gov/blast/db/
https://www.uniprot.org/downloads

1)nr是ncbi收集的目前所有微生物的蛋白序列，是用來(lái)計(jì)算氨基酸一般情況下的頻率的，160G

2)uniprot90根據(jù)相似性做了一個(gè)去冗余，所以比nr要小很多，56G

# 以u(píng)niprot90為例
wget ftp://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz # 下載
gzip -d uniref90.fasta.gz # 解壓
makeblastdb -in uniref90.fasta -parse_seqids -hash_index -dbtype prot # 編譯

解析完成后的樣子：

文件是這個(gè)樣子：（只截取了一部分）

3. 獲取PSSM矩陣

我的初始文件是：

P00269.fasta是對(duì)單條蛋白質(zhì)處理，里面的格式是：

testset.fasta是對(duì)蛋白質(zhì)集合批處理，里面的格式是（也可以單獨(dú)蛋白質(zhì)存為.fasta文件，由于blast只能處理單條蛋白糊，把這個(gè)集合知識(shí)歸總的意思，第一步還是要生成單條蛋白質(zhì)的.fasta文件，所以這個(gè)文件看個(gè)人意愿）：

1）單條蛋白質(zhì)序列的處理方法

import os
os.system('psiblast -query dataset/P00269.fasta -db /PSSM/uniref90.fasta -num_iterations 3 -out_ascii_pssm /dataset/P00269.pssm')##這個(gè)蛋白質(zhì)好慢呀

2）批處理獲取的方法

import os
 
file_name='/dataset/testset.fasta'
Protein_id=[]
with open(file_name,'r') as fp:
    i=0
    for line in fp:
 
        if i%2==0:
            # Protein_id.append(line[1:-1])
            id=line[0:-1]
            p=line[1:-1]
            with open ('/dataset/'+str(p)+'.fasta','a') as protein:
                protein.write(id)
                # protein.write()
        if i%2==1:
            seq=line[0:-1]
            with open ('/dataset/'+str(p)+'.fasta','a') as protein:
                protein.write('\n')
                protein.write(seq)
        i=i+1
 
        os.system('psiblast -query '+'/dataset/'+str(p)+'.fasta -db /PSSM/uniref90.fasta -num_iterations 3 -out_ascii_pssm /dataset/'+str(p)+'.pssm')

##PSSM真是太慢了，下面是只生成一個(gè)后的截圖