171110 blast

命令

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输出结果过滤

link

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)

3.结果进一步注释