个性化本地blast

这是这个教程的延续：如何批量获取基因的转录组序列
在进行本地blast之前需要制作自己的比对数据库，注意所提供的fasta的header长度不能超过50字符，因此从genecode中提取出来的fasta序列得改下header的长度，下面是我用的方法。最好把需要检测的fasta文件的名字也改了。
更改fasta Header教程
blast参考教程

###Blast
ml blast
cd ${data_dir}
# rename the mouse fasta file because the db just allow header length less than 50
#${mouse_fa_dir}是我的fasta文件目录，请按照自己的更改
awk '{ if ($0~/^>/) { n=split($0, a, "|"); gsub(/_/," ", a[1]); printf("%s|%s\n", a[1], substr(a[2], 2)); } else { print $0; } }' ${mouse_fa_dir}/mouse.lncRNA.fa>${mouse_fa_dir}/mouse.lncRNA.renmae.fa
blast_db_dir=/your_dir/BLAST_DB
makeblastdb -in ${mouse_fa_dir}/mouse.lncRNA.renmae.fa -dbtype nucl -out ${blast_db_dir} -parse_seqids

###blast human to mouse
ml blast
data_dir=/your_dir/sequence
human_blast_result=/your_dir/human/blast_results
blast_db_dir=/your_dir/BLAST_DB
mkdir -p ${human_blast_result}

cat ${data_dir}/human_genelist.txt | while read line
do
blastn -db ${blast_db_dir} -query ${human_fa_dir}/${line}.fa  -out ${human_blast_result}/${line}.blast2mouse.tsv
done

#blastn -db ${blast_db_dir} -query ${mouse_fa_dir}/mouse.lncRNA.renmae.fa  -out ${human_blast_result}/test.blast2mouse.tsv

blast_results_analysis=/your_dir/human/blast_results_analysis
#check the hit genes
cd  ${human_blast_result}
grep -H "Identities" *.tsv > ${blast_results_analysis}/human_blast2mouse_summary.tsv