linux下用psiblast批量生成pssm矩陣

目錄

• 在linux上安裝psiblast
• 下載下傳并編譯用于比對的大型蛋白質資料庫
• 生成pssm
• 批量生成方法
• 文章中看到的手動生成pssm的方法

在linux上安裝psiblast

最好建立一個python環境，因為我發現conda安裝blast預設的是python==3.6.11，可能會不小心把你的python版本改掉…然後你寫好的代碼全die了……

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

下載下傳并編譯用于比對的大型蛋白質資料庫

ftp://ftp.ncbi.nlm.nih.gov/blast/db/

https://www.uniprot.org/downloads

nr和uniprot是比較通用的資料庫：

• nr是ncbi收集的目前所有微生物的蛋白序列，是用來計算氨基酸一般情況下的頻率的，160G
• uniprot90根據相似性做了一個去備援，是以比nr要小很多，56G

# 以uniprot90為例
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 # 編譯
           

溫馨提示：

1. nr有點太大了，我在makeblastdb的時候由于存儲空間不夠，Failed了。當時還有700G左右的空間……是以我選擇了uniprot90，成功了~

   

2. swissprot那個資料庫雖然很小，但是估計sequence太少了，導緻我一些序列會被忽略，就是根本不進行pssm的計算，啥也得不到……

   我的猜測是不對的，因為uniprot90也不想，哪個大神指點一下吧~

生成pssm

psiblast -query LYTMKGPDWNTELYN.fasta -db uniref90.fasta -num_iterations 3 -out_ascii_pssm LYTMKGPDWNTELYN.pssm
           

• -query：需要生成pssm的fasta檔案
• -db：用于比對的大資料庫的名稱
• -out_ascii_pssm：生成的pssm檔案名

批量生成方法

用python寫個.sh腳本，然後用bash .sh指令

def sh_text(sequences):
    for item in tqdm(sequences):
        item = str.upper(item) # 因為我的fasta裡有小寫字母……
        fasta_file = item + '.fasta'
        pssm_file = item + '.pssm'
        s = 'psiblast -query ' + fasta_file + ' -db uniref90.fasta -num_iterations 3 -out_ascii_pssm ' +  pssm_file + '\n'
        with open('pssm_blast.sh', 'a') as f:
            f.write(s)
           

文章中看到的手動生成pssm的方法

文章名：PSSMHCpan_ a novel PSSM-based software for predicting class I peptide-HLA binding affinity