1. 对小鼠的蛋白序列进行建库
点击查看代码
makeblastdb -in Mus_musculus.GRCm39.pep.all.fa -dbtype prot -parse_seqids -out ./mm.blast.db/index
2. 将自己的物种的蛋白序列进行比对
点击查看代码
blastp -query pep.fa -db ./mm.blast.db/index -evalue 1e-6 -outfmt 6 -num_threads 6 -out out1.aglin.txt
### outfmt可以提取小鼠对应的蛋白序列
blastp -query pep.fa -db ./mm.blast.db/index -evalue 1e-6 -outfmt '6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore' -num_threads 12 -out out1.aglin.txt
3. 根据out1.aglin.txt中结果提取evalue值最小的记录
点击查看代码
awk '!a[$1]++{print}' out1.aglin.txt > uniq.blast
4. 根据提取到最优的小鼠蛋白名提取序列
点击查看代码
awk -F '\t' '{print $2,$12}' OFS="\n" out1.aglin.txt|sed 's/ENSMUSP/>ENSMUSP/g' >> query_2.fa
awk 'NR==FNR{a[$2];next} $2 in a {print $0}' pep.fa out2.aglin.txt | head