blast和dot生物信息工具

blast主要方法使用总结（Linux）

2023-01-15  本文已影响0人  马尔代夫Maldives

BLAST的主要理念

blast是区段比对，对于给定的两个序列，blast会把具有相似性的片段（hit）找出来，显示的是hit的信息。

一、 linux 安装BLAST

linux下的blast是一系列工具集合，包括：blastn, blastp, blastx, tblastn, tblastx等。
采用conda安装，可能出现网络连接超时等问题，重复尝试：

# 安装blast
$>>conda install -c bioconda blast

# blast安装perl模块的方法
$>>conda install perl-digest-md5

二、 BLAST使用方法

第（1）步：makeblastdb命令建立检索所需“数据库”

blast数据库有两种类型：核酸数据库（关键字：nucl）和氨基酸数据库（关键字：prot）。
blast创建数据库的命令：makeblastbd（该命令的参数可以通过“-help”查看，这很重要！）

$>>makeblastdb -help              
USAGE
  makeblastdb [-h] [-help] [-in input_file] [-input_type type]
    -dbtype molecule_type [-title database_title] [-parse_seqids]
    [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
    [-mask_desc mask_algo_descriptions] [-gi_mask]
    [-gi_mask_name gi_based_mask_names] [-out database_name]
    [-blastdb_version version] [-max_file_sz number_of_bytes]
    [-logfile File_Name] [-taxid TaxID] [-taxid_map TaxIDMapFile] [-version]

DESCRIPTION
   Application to create BLAST databases, version 2.9.0+

REQUIRED ARGUMENTS
 -dbtype <String, `nucl', `prot'>
   Molecule type of target db
……

示例：

以拟南芥氨基酸数据为例（拟南芥基因组数据方法类似），下载地址：ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
下载后解压得到氨基酸序列文件（就是一个普通的fasta格式文件）：Arabidopsis_thaliana.TAIR10.pep.all.fa

1）构建拟南芥氨基酸的blast数据库命令：

$>>makeblastdb -in  Arabidopsis_thaliana.TAIR10.pep.all.fa -input_type fasta -dbtype prot -out TAIR10 -parse_seqids

Building a new DB, current time: 01/13/2023 09:34:33
New DB name:   /data/home/chengaoxiang/blast_data_cgx/TAIR10
New DB title:  Arabidopsis_thaliana.TAIR10.pep.all.fa
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 48321 sequences in 1.11516 seconds.

数据库建成后，生成了下图中的几个文件：


生成数据库文件.png

2）查看blast数据库信息：

第（2）步：选择blast中不同的工具进行计算

根据不同的需求，比如说你的待比对序列是氨基酸还是核苷酸，搜索数据库是核甘酸还是氨基酸，选择合适的blast工具。不同需求的对应关系可以见下图(来自biostars handbook)。
（不同工具的应用范围虽然不同，但是基本参数都是一致的，学会blastn，基本上其他诸如blastp，blastx也都会了。）


根据待比对序列类型和数据库数据类型选择对应方法.png

（延伸知识：一个DNA顺序可能有3种阅读框，但只有一种具有编码作用，称为开放阅读框（open reading frame or ORF）。有的阅读框因终止密码出现频繁，故不能生成蛋白，这种阅读框称为封闭阅读框（block reading frame）。若一个顺序所有的三个阅读框都是封闭的，则它无编码蛋白的功能。一个翻译成蛋白的顺序有一个阅读框，开始于AUG起始密码子，通过一系列有义密码子，直到终止密码子结束。通常3个阅读框中总有封闭阅读框的存在。
六种阅读框架：一段DNA或RNA序列有多种不同读取方式，因此可能同时存在许多不同的开放阅读框架。开放阅读框包含一段可以编码蛋白的碱基序列，不能被终止子打断。 当一个新基因被识别，其DNA序列被解读，人们仍旧无法搞清相应的蛋白序列是什么。这是因为在没有其它信息的前提下，DNA序列可以按六种框架阅读和翻译（核酸有两条链，每条链三种，对应三种不同的起始密码子）。）

示例：

以人血红蛋白alpha亚基（HBA_HUMAN.FASTA，内含1个氨基酸序列）作为查询序列（query sequence）（所谓查询序列就是我们要用于在数据库中比对的序列，目的在于从数据库中找到与查询序列相似的序列），如下所示，：

$>>cat HBA_HUMAN.FASTA 
>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1 
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR

然后在前面建的TAIR10数据库中进行比对，根据上表，应采用blastp方法进行对比，比对时有很多参数可以选择，下面挑重点的几项说明：

1）最简洁比对命令

$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 | less

比对结果：
列出了所有被筛出来的氨基酸序列（6个），同时后面分别列出了每个氨基酸的详细信息。

BLASTP 2.9.0+
……
                                                                      Score        E
Sequences producing significant alignments:                          (Bits)     Value

AT4G37710.1 pep chromosome:TAIR10:4:17718118:17718742:1 gene:AT4G...  31.2       0.17 
AT4G37710.2 pep chromosome:TAIR10:4:17718128:17718771:1 gene:AT4G...  31.2       0.17 
AT2G03880.1 pep chromosome:TAIR10:2:1181512:1183703:1 gene:AT2G03...  29.6       0.99 
AT4G16420.2 pep chromosome:TAIR10:4:9262621:9266157:-1 gene:AT4G1...  26.9       9.7  
AT4G16420.3 pep chromosome:TAIR10:4:9262621:9266086:-1 gene:AT4G1...  26.9       9.7  
AT4G16420.1 pep chromosome:TAIR10:4:9262621:9266157:-1 gene:AT4G1...  26.9       9.7  

>AT4G37710.1 pep chromosome:TAIR10:4:17718118:17718742:1 gene:AT4G37710 transcript:AT4G37710.1 
gene_biotype:protein_coding transcript_biotype:protein_coding 
gene_symbol:VQ29 description:VQ motif-containing 
protein 29 [Source:UniProtKB/Swiss-Prot;Acc:Q9SZG3]
……

2）将比对结果保存在文件中命令（关键参数-out）：

$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -out HBA_SW.TXT

3）选择输出结果格式命令（关键参数-outfmt）：

$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -outfmt 7 | less

执行结果：

# BLASTP 2.9.0+
# Query: sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
# Database: /data/home/chengaoxiang/blast_data_cgx/TAIR10
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 6 hits found
sp|P69905|HBA_HUMAN AT4G37710.1 48.276  29  15  0   78  106 45  73  0.17    31.2
sp|P69905|HBA_HUMAN AT4G37710.2 48.276  29  15  0   78  106 45  73  0.17    31.2
sp|P69905|HBA_HUMAN AT2G03880.1 31.579  38  26  0   76  113 41  78  0.99    29.6
sp|P69905|HBA_HUMAN AT4G16420.2 27.419  62  39  2   15  76  125 180 9.7 26.9
sp|P69905|HBA_HUMAN AT4G16420.3 27.419  62  39  2   15  76  125 180 9.7 26.9
sp|P69905|HBA_HUMAN AT4G16420.1 27.419  62  39  2   15  76  125 180 9.7 26.9
# BLAST processed 1 queries

上述输出格式每列的意义：
query acc.ver（查询序列）、subject acc.ver（结果序列）、% identity（相似度）、alignment length（比对长度）、mismatches（错配数）、gap opens（起始空位数）、q. start（查询序列起始位点）、q. end（查询序列终止位点）、s. start（结果序列起始位点）、s. end（结果序列终止位点）、evalue（期望值）、bit score（得分）。

4）★设置evalue值（非常重要，前面几个示例没有设置该参数，实际中得设置）

将3）中的命令增加“-evalue 0.2”，表示只输出evalue小于0.2的匹配结果。从输出结果可以看出，只输出了evalue=0.17的两条序列。

$>>blastp -query /data/home/chengaoxiang/blast_test/HBA_HUMAN.FASTA -db /data/home/chengaoxiang/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 | less

输出结果：

# BLASTP 2.9.0+
# Query: sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
# Database: /data/home/chengaoxiang/blast_data_cgx/TAIR10
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 2 hits found
sp|P69905|HBA_HUMAN     AT4G37710.1     48.276  29      15      0       78      106     45      73      0.17    31.2
sp|P69905|HBA_HUMAN     AT4G37710.2     48.276  29      15      0       78      106     45      73      0.17    31.2
# BLAST processed 1 queries

blast对比对结果的评价主要有两个参数：evalue（E值）和bit score（S值）

E值的意义：
在一次blast分析过程中，假设query序列与subject序列（数据库序列）命中了N个hit（命中片段），这N个hit一定可信吗？不一定，因为即使随机产生一组序列，也有可能和query序列存在hit。基于这点，生成一个大小与原数据库相同的随机数据库，然后计算query与该随机数据库之间的hit（假设得到了M个hit）。最后判断M中有多少个hit与N中的hit相同或更好，假设有H个。H应该越小越好，H越小说明N中的hit都是通过blast方法得到的，而不是随机得到的。E值就是这种随机性的度量（公式详见后面），与这里的H等价，因此E也是越小越好，它越小，说明被随机匹配到的可能性越小。这就是E的意义。

E经验值：
E<1e-5：匹配结果基本可以确保是同源关系，这是一个非常高的置信度；
1e-5<E<0.01：匹配结果可以视为同源性关系；
0.01<E<10：匹配结果基本没有同源性，但可能存在暂时的远程同源关系，需要额外的证据来确认；
10<EE：没有同源性，或极远的同源性，但已经低于当前方法所能达到的检测极限了。

E值和S值的计算方式：
S=？？？？？？？？？？？？？（待增加）：比特分数是序列数据库的所需大小，其中当前匹配可能只是偶然发现。比特分数是log 2缩放和标准化的原始分数。每增加一倍，所需的数据库大小（2的"位分数"次方）就会增加一倍。
E=K*m*n*(e-λ*S)：（K和λ与数据库和算法有关，是常量；m是目标序列长度，n是数据库大小，S就是前面提到的S值。）【根据公式，对于某个给定的query序列，其E值随数据库增大而增大（因为参数K和λ），但该query序列和原数据库序列中某个序列的进化关系其实并未改变。这意味着，随着数据库的增大，相同query和数据库序列对之间的匹配可信度会降低，即随着数据库的增大，原来被检测出的有效结果可能会被淹没。因此，需要一种替代E值计算的方法。】

E值的局限性：

  1. 当目标序列过短时，由于无法得到较大的S值，导致E值偏大。
  2. 当两序列同源性虽然高，但有较大的gap（空隙）时，S值会下降。这个时候gap scores就非常有用。
  3. 有些序列的非功能区有较低的随机性时，可能会造成两序列较高的同源性。

5）设置字长“-word_size”

$>>blastp -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -word_size 2 | less

6）设置打分矩阵（-matrix）

$>>blastp -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -matrix PAM250 | less

7）psiblast（循环搜索命令） : 敏感度更高的蛋白序列与蛋白序列之间的比对

Position-Specific Iterated (PSI)-BLAST，是一种更加高灵敏的Blastp程序，对于发现远亲物种的相似蛋白或某个蛋白家族的新成员非常有效。当你使用标准的Blastp比对失败时，或比对的结果仅仅是一些假基因或推测的基因序列时（”hypothetical protein” or “similar to…”），你可以选择PSI-BLAST重新试试。

$>>psiblast -query xx/HBA_HUMAN.FASTA -db xx/blast_data_cgx/TAIR10 -evalue 0.2 -outfmt 7 -num_iterations 2 | less

8）双序列比对，指定起始位点

以癌胚抗原搜索结构域为例（网上例子）：

$>>blastp -query CEA21_HUMAN.FASTA -query_loc 147-231 -subject CEAM5_HUMAN.FASTA -subject_loc 145-675 -outfmt 7 | less
双序列比对.png

9）PHI-BLAST : 模式发现迭代BLAST

PHI-BLAST, 模式发现迭代BLAST, 用蛋白查询来搜索蛋白数据库的一个程序。仅仅找出那些查询序列中含有的特殊模式的对齐。
PHI的语法详细介绍看这里：http://www.ncbi.nlm.nih.gov/blast/html/PHIsyntax.html （不过，不是太明白）

10）NCBI-blast（在NCBI上blast）

在百度查询；
https://wenku.baidu.com/view/7a716acca1c7aa00b52acb14?aggId=138471d707a1b0717fd5360cba1aa81144318fc3&fr=catalogMain&wkts=1673851517291&bdQuery=blast%E6%AF%94%E5%AF%B9%E5%AD%97%E9%95%BF
https://blast.ncbi.nlm.nih.gov/Blast.cgi

参考资料：

（重要网址）
https://blog.csdn.net/withbeginner/article/details/124126747?spm=1001.2101.3001.6650.15&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-15-124126747-blog-112264605.pc_relevant_3mothn_strategy_and_data_recovery&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-15-124126747-blog-112264605.pc_relevant_3mothn_strategy_and_data_recovery&utm_relevant_index=18

（blast中文手册）
https://blog.csdn.net/yin1331102028yin/article/details/126285678?spm=1001.2101.3001.6650.11&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-11-126285678-blog-88992486.pc_relevant_multi_platform_whitelistv3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-11-126285678-blog-88992486.pc_relevant_multi_platform_whitelistv3&utm_relevant_index=12

https://xuzhougeng.blog.csdn.net/article/details/115511912?spm=1001.2101.3001.6650.4&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-4-115511912-blog-112264605.pc_relevant_3mothn_strategy_and_data_recovery&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-4-115511912-blog-112264605.pc_relevant_3mothn_strategy_and_data_recovery&utm_relevant_index=7

http://www.360doc.com/content/21/0714/12/76149697_986499307.shtml
https://www.jianshu.com/p/e5527735f163?utm_campaign=maleskine&ivk_sa=1024320u
https://www.metagenomics.wiki/tools/blast/evalue
https://blog.csdn.net/weixin_43364556/article/details/107099313
https://blog.csdn.net/withbeginner/article/details/124125746
http://blog.shenwei.me/local-blast-installation/
https://www.ncbi.nlm.nih.gov/books/NBK62051/

上一篇下一篇

猜你喜欢

热点阅读