下载数据
#将数据改名为Arabidopsis_thaliana.fna
mv NC_003070.fna Arabidopsis_thaliana.fna
#我下载的是拟兰介的1号染色体，取其前100000行作为我的测试数据
head -n 100000 Arabidopsis_thaliana.fna > sample.01.fna
下载blast
cd ~
#下载本地blast
nohup wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.6.0+-x64-linux.tar.gz &
#解压，重命名,环境变量的配置
sudo tar zxvf tar zxvf ncbi-blast-2.6.0+-x64-linux.tar.gz -C /opt/biosoft/ #解压
sudo mv /opt/biosoft/ncbi-blast-2.6.0+ /opt/biosoft/ncbi-blast #重命名
echo 'PATH=/opt/biosoft/ncbi-blast/bin:$PATH' >> ~/.bashrc #加入环境变量
source ~/.bashrc
blastn -version
cd ~ && vim .ncbirc
#.ncbirc的文件中输入
[BLAST]
BLASTDB=/home/train/blast/db
现在才是重点：
Blastn要求格式化好的数据库，从NCBI ftp站点下载下来的fasta格式的数据库不能直接使用，需要使用makeblastdb程序来格式化，命令如下：
../ncbi/bin/makeblastdb -in Arabidopsis_thaliana.fna -parse_seqids -hash_index -dbtype nucl
-in参数后面接将要格式化的数据库，-parse_seqids, -hash_index两个参数一般都带上，主要是为blastdbcmd取子序列时使用，-dbtype nucl告诉程序这是核酸数据库。
检查数据库是否正确
../ncbi/bin/blastdbcheck -dir ./
运行blastn程序，
../ncbi/bin/blastn -task blastn -query sample.01.fna -db Arabidopsis_thaliana.fna out test.out
-task blastn：任务类型，前面已有描述；
-query NGB.rna：查询序列，必须是fasta格式；
-db human.rna.fna：格式化好的查询数据库；
-out NGB_blast.txt：结果输出文件。
查看结果
less test.out|more
real 18m58.015s
user 18m53.860s
sys 0m2.861s
以下是比对结果的前几行：
BLASTN 2.2.27+
Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A.
Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J.
Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of
protein database search programs", Nucleic Acids Res. 25:3389-3402.
Database: Arabidopsis_thaliana.fna
1 sequences; 30,427,671 total letters
Query= gi|240254421|ref|NC_003070.9| Arabidopsis thaliana chromosome 1,
complete sequence
Length=6999930
Score E
Sequences producing significant alignments: (Bits) Value
ref|NC_003070.9| Arabidopsis thaliana chromosome 1, complete seq... 1.262e+07 0.0
>ref|NC_003070.9| Arabidopsis thaliana chromosome 1, complete sequence
Length=30427671
Score = 1.262e+07 bits (13999860), Expect = 0.0
Identities = 6999930/6999930 (100%), Gaps = 0/6999930 (0%)
Strand=Plus/Plus
Query 1 ccctaaaccctaaaccctaaaccctaaacctctgaatccttaatccctaaatccctaaat 60
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1 CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAAT 60
sample.01.fna中一共7099929个碱基，
数据库中一共30,427,671个碱基，
比对花费时间：
real 18m58.015s
user 18m53.860s
sys 0m2.861s
BLAST中的E值(E-value)是什么意思？
BLAST是指Basic Local Alignment Search Tool，是生物信息学中的一种序列比对算法，用于寻找蛋白质或核酸的相似序列。
下面是一个BLAST结果，
| Sequences producing significant alignments: | Score (S) | E |
| gi|83574104|Moth_2374|sporulation – prote… | 202 | 2e-53 |
| gi|83573446|Moth_1696|Sporulation – prote… | 112 | 1e-26 |
| gi|83571874|Moth_0087|sporulation – prote… | 95 | 3e-21 |
| gi|83573435|Moth_1685|Substrate-binding -… | 27 | 1.0 |
后面有两个值，一个是S值，一个E值。可以发现，结果是依据S值的高低来显示的。
S值表示两序列的同源性，分值越高表明它们之间相似的程度越大。
E值就是S值可靠性的评价。它表明在随机的情况下，其它序列与目标序列相似度要大于S值的可能性。所以它的分值越低越好。
E值的计算：
E=K*m*n*(e-lambda*S)
其中，K和lambda与数据库和算法有关，是个常量；m代表目标序列的长度，n代表数据库的大小，S就是前面提到的S值。
通常来讲，我们认为E值小于10-5就是比较可性的S值结果。我们可以想象，相同的数据库，E=0.001时如果有1000条都有机会比现在这个S值要高的话，那么把E设置为10-6时可能就会只得到一条结果，就是S值最可靠的那个。
但是E值也不是万能的。它在以下几个情况下有局限性：
1. 当目标序列过小时，E值会偏大，因为无法得到较高的S值。
2. 当两序列同源性虽然高，但有较大的gap（空隙）时，S值会下降。这个时候gap scores就非常有用。
3. 有些序列的非功能区有较低的随机性时，可能会造成两序列较高的同源性。
E值总结：
E值适合于有一定长度，而且复杂度不能太低的序列。
当E值小于10-5时，表明两序列有较高的同源性，而不是因为计算错误。
当E值小于10-6时，表明两序列的同源性非常高，几乎没有必要再做确认。