命令
makeblastdb -in CAZyDB.07202017.fa -dbtype prot -out CAZydb
nohup blastp -db ~/database/pfam/Pfam-A-duplicate.fasta -query ../assembly.faa -num_threads 8 -out pfam.out -outfmt "6" -evalue 1e-5 > pfam.log 2>&1 &
pfam nr trembl很耗时间
1. 关于每个gene id只保留最优hit
1.1 方法一
sort -k1,1 -k12,12gr -k11,11g -k3,3gr blastout.txt | sort -u -k1,1 --merge > bestHits
1.2 方法二
blast加参数-max_target_seqs 1
1.3 结果对比
phi.out 13213
phi_1.out :blast调参数 1579 （仍有少量第一行重复）
phi_bestHits: sort 1453 (去重复更严格）
2. blast输出结果过滤
import pandas as pd
prote_df= pd.read_csv('blast',delimiter='\t',names=('qseqid','sseqid','pident','length','mismatch','gapopen','qstart','qend','sstart','send','evalue','bitscore','qlen','slen'))
blast_flit=prote_df[(prote_df['pident']>60)&(prote_df.length/prote_df.qlen>0.80)]
blast_flit.to_csv('gene.csv',index=False)
