python 动态规划算法 蛋白质对比序列 蛋白质序列比对blast

主要内容

1 背景

2 在线blast

3 本地blast

3.1 老版本blast

3.2 新版本blast

背景

序列比对（Sequence Alignment）的基本问题是比较两个或两个以上序列的相似性。

blast作为一种序列相似性比对工具，是生物信息分析最常用的一款软件，必须掌握。不管是做两序列相似性的简单比对，还是引物特异性、序列的来源等个性化分析，都会用到blast比对。许多看似高大上的基因分析，都可归类于序列间的比较，因此blast是生信分析中基础性的工具。

如果你有一堆测序回来的序列，想要看看它们是来自于哪个物种的，或者想在数据库中搜索对应的同源序列，使用ncbi的在线blast，数据库很全，速度很快，马上就能知道结果。

1. blastp：蛋白序列与蛋白库作比对，直接比对蛋白序列的同源性。
2. blastx：核酸序列与蛋白库作比对，将核酸序列先翻译成蛋白序列，再将其与蛋白库作比对。
3. blastn：核酸序列与核酸库的比对，直接比对核酸序列的同源性。
4. tblastn：蛋白序列对核算库的比对，现将核酸库翻译成蛋白库，再将蛋白序列与翻译后的蛋白库进行比对。
5. tblastx：核酸与核酸数据库在蛋白质水平比较

同源性（homology）VS 一致性（identity）
同源性是来描述物种之间的进化关系的，所以在同源性的表达中只能用“有”或者“无”，对于有同源性的物种可以描述为“部分同源”或者“完全同源”。

有些小伙伴们会说“序列A和序列B之间有85%的同源性”，这种说法是不正确的，A和B之间要么有同源性，要么没有同源性，可以这样说：序列A和序列B之间有85%的identity，A和B之间有同源性。

在线blast

blast网站 https://blast.ncbi.nlm.nih.gov/Blast.cgi

Nucleotide BLAST

核酸序列到核酸库中的一种查询，库中存在的每条已知序列都将同所查序列作一对一的核酸序列比对。

Protein BLAST

蛋白序列到蛋白库中的一种查询，库中存在的每条已知序列将逐一地同每条所查序列作一对一的序列比对。

BLASTX

核酸序列到蛋白库中的一种查询，先将核酸序列翻译成蛋白序列，再对翻译成的每一条序列作一对一的蛋白序列比对。

TBLASTN

蛋白序列到核酸库的一种查询，与BLASTX相反，它是将库中的核酸序列翻译成蛋白质序列，再对所查序列作蛋白与蛋白的比对。

我们点击“Nucleotide BLAST”后就到了blastn的初始界面，输入以下序列

GACGCGGCCGTCGAGGGCGCTCAGGTGGACTTCGACACCGCCAGCCGCATCGAGAGCCGCTACTTCACCCAGCTGGTCACCGGCCAGGTCGCCAAGAACATGATCCAGGCGTTCTTCTTCGACCTCCAGCACATCAACGGCGGCGGCTCCCGCCCCGAGGGCATCGAGCCGGTCAAGATCAACAAGATCGGTGTGCTCGGCGCGGGCATGATGGGCGCCGGCATCGCCTACGTCTCGGCCAAGGC

数据库可以有很多选择，默认是nr/nt库

点击“BLAST”即可进行比对，20秒之内就会出现如下结果：

发现这是一段来 Mycobacterium avium 的序列。

其中，Expect(E值)和Identities（一致性）是评价blast结果的标准。E值接近零或者为零时，说明比上的序列很接近；一致性：匹配上的碱基数占总序列长的百分数。

如果比对其他序列出现如下界面，说明没比上，可以试着选择其他参数

可以选择 Somewhat similar sequences (blastn) 进行再次比对

MEGABLAST : 采用贪婪式算法，多用于比较相似性比较高（相似性在>95%）的序列，灵敏度高，速度快。

Discontiguous MEGABLAST : 灵敏度更高，用于更精确的序列的比对。主要用于跨物种之间的同源比对。

BlastN : 用于比对相似性较差的序列，，相似度较低的序列也可以查找到，比对结果最多，速度最慢。 但允许更短序列的比对（如短到7个碱基的序列），例如做短的引物的比较可以使用这个选项。

如果想比较手头上的序列A和序列B之间的相似性，可以在blastn的初始界面点击
Align two or more sequences，就可以分别放入序列进行比对了。

本地blast

老版本blast

一些老的服务器上面自带的只有老版本的blast

很简单的两个步骤：formatdb 建库和 blast 比对

首先建库，把database.fasta这个序列文件变成blast专用的库，-p选项中的T是代表蛋白库，F是代表核酸库

1. !formatdb -i database.fasta -p F -o T

然后需要用-p来选择比对程序，核酸比核酸是blastn

1. !blastall -p blastn -d database.fasta -i query.fasta -o query.fasta.bls -v 1 -b 1 -e 1e-10 -F F
2. !head -80 query.fasta.bls
3. BLASTN 2.2.26 [Sep-21-2011]
4. Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
5. Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
6. "Gapped BLAST and PSI-BLAST: a new generation of protein database search
7. programs", Nucleic Acids Res. 25:3389-3402.
8. Query= query
9. (1026 letters)
10. Database: database.fasta
11. 503 sequences; 3,420,533 total letters
12. Searching..................................................done
13. Score E
14. Sequences producing significant alignments: (bits) Value
15. gnl|test|test_1 2034 0.0
16. >gnl|test|test_1
17. Length = 883888
18. Score = 2034 bits (1026), Expect = 0.0
19. Identities = 1026/1026 (100%)
20. Strand = Plus / Plus
21. Query: 1 agatctgacggcgggttgatcgccatgcctcgtttccttcttttgttaaatgcaccataa 60
22. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
23. Sbjct: 8776 agatctgacggcgggttgatcgccatgcctcgtttccttcttttgttaaatgcaccataa 8835
24. Query: 61 taggtcaaacttgactgaaacttgtgggcgaattcagctccgcttttgcgactcatccga 120
25. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
26. Sbjct: 8836 taggtcaaacttgactgaaacttgtgggcgaattcagctccgcttttgcgactcatccga 8895
27. Query: 121 tgttttgcggctgacaaattcttcgatcaggctggcaagctcgcgccgctcacggtcagc 180
28. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
29. Sbjct: 8896 tgttttgcggctgacaaattcttcgatcaggctggcaagctcgcgccgctcacggtcagc 8955
30. Query: 181 ctgcgcatcgttcagccggtcttgcagaagttccacgcgccgttgctctttcaaaagtgc 240
31. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
32. Sbjct: 8956 ctgcgcatcgttcagccggtcttgcagaagttccacgcgccgttgctctttcaaaagtgc 9015
33. Query: 241 ttttcgcgccaattcgagccgtgattcgacctgttgcatttcaaccgcattggtctccag 300
34. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
35. Sbjct: 9016 ttttcgcgccaattcgagccgtgattcgacctgttgcatttcaaccgcattggtctccag 9075
36. Query: 301 ccaccggatgacgagcgagggatcactatcgagcgcattcgcgtcatagcgccggtcctg 360
37. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
38. Sbjct: 9076 ccaccggatgacgagcgagggatcactatcgagcgcattcgcgtcatagcgccggtcctg 9135
39. Query: 361 catggcaatcagcgccgcgcgctcctcatcaagcgccgactttcgcgcatttgcgatggc 420
40. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
41. Sbjct: 9136 catggcaatcagcgccgcgcgctcctcatcaagcgccgactttcgcgcatttgcgatggc 9195
42. Query: 421 aagtgctgcttcgtgccgcgcacagccaagcttttgcagttccaccagtttcctcagatt 480
43. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
44. Sbjct: 9196 aagtgctgcttcgtgccgcgcacagccaagcttttgcagttccaccagtttcctcagatt 9255
45. Query: 481 gcgtgttccagcctgcgtcatgtcgccctcagaaaagttcgcgcaactgggttgcggttt 540
46. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
47. Sbjct: 9256 gcgtgttccagcctgcgtcatgtcgccctcagaaaagttcgcgcaactgggttgcggttt 9315
48. Query: 541 catgcgtgaaaaagcgcagaagctccggcaggatgaaatagacaaccagcagccctcctc 600
49. ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
50. Sbjct: 9316 catgcgtgaaaaagcgcagaagctccggcaggatgaaatagacaaccagcagccctcctc 9375

我们输出格式的m8，这个格式很重要，我们还可以设置-a控制cpu数量，和-e控制阈值等

1. -p：blastn/blastp/blastx/tblastn/tblastx 分别对应不同的数据比对方式
2. -e：指定一个实数，过滤掉期望值大于这个数的比对结果，默认为10，建议蛋白比对设置为1E-5，核酸比对设为1E-10
3. -F：用来屏蔽简单重复和低复杂度序列( T/F) ，默认为T。可提高比对的精确度，但对于引物特异性检测等分析，建议设置为F，不然获得的identity可能会有问题。
4. -m：设定输出格式，-m 0~6展示了subjects间的比对结果，-m8~9以表格形式展示比对结果，默认0。
5. -v：输出中每一个query的比对列表最多显示subject个数
6. -b：每个query 最多显示与多少条 subject 的比对条形图
7. !blastall -p blastn -d database.fasta -i query.fasta -o query.fasta.bls -a 8 -v 1 -b 1 -e 1e-10 -F F -m 8
8. !head -80 query.fasta.bls
9. query gnl|test|test_1 100.00 1026 0 0 1 1026 8776 9801 0.0 2034

新版本blast

也是两个步骤：makeblastdb 建库和 blast 比对

1. !makeblastdb -in database.fasta -dbtype nucl -parse_seqids -out databasename
2. Building a new DB, current time: 03/01/2019 16:55:18
3. New DB name: databasename
4. New DB title: database.fasta
5. Sequence type: Nucleotide
6. Keep Linkouts: T
7. Keep MBits: T
8. Maximum file size: 1000000000B
9. Adding sequences from FASTA; added 503 sequences in 0.149536 seconds.
10. 参数说明:
11. -in：待格式化的序列文件
12. -dbtype：数据库类型，prot或nucl
13. -out：数据库名

blastn比对

---

---

1. !blastn -query query.fasta -out query.fasta.bls -db databasename -outfmt 6 -evalue 1e-10 -num_threads 8
2. 参数说明:
3. -query： 输入文件名
4. -out：输出文件名
5. -db：格式化了的数据库名
6. -outfmt：输出文件格式，总共有12种格式，6是tabular格式对应BLAST的m8格式
7. -evalue：设置输出结果的e-value值
8. -num_threads：线程数
查看结果
1. !more query.fasta.bls
2. query gnl|test|test_1 100.00 1026 0 0 1 1026 8776
3. 9801 0.0 1895