这或许是我写的最全的BLAST教程

简介: Basic local alignment search tool (BLAST)包括:blastn, blastp, blastx, tblastn, tblastx等. 使用conda安装即可。

Basic local alignment search tool (BLAST)

包括:blastn, blastp, blastx, tblastn, tblastx等. 使用conda安装即可。

conda install -c bioconda blast
# blast安装perl模块的方法
conda isntall perl-digest-md5

BLAST的主要理念

  • Search may take place in nucleotide and/or protein space or translated spaces where nucleotides are translated into proteins.
  • Searches may implement search “strategies”: optimizations to a certain task. Different search strategies will return different alignments.
  • Searches use alignments that rely on scoring matrices
  • Searches may be customized with many additional parameters. BLAST has many subtle functions that most users never need.

本地BLAST的基本步骤

  1. 用makeblastdb为BLAST提供数据库
  2. 选择blast工具,blastn,blastp
  3. 运行工具,有必要的还可以对输出结果进行修饰

第一步:建立检索所需数据库

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]
    [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
    [-taxid_map TaxIDMapFile] [-version]
-dbtype <String, `nucl', `prot'>

具体以拟南芥基因组作为案例,介绍使用方法:
: 拟南芥的基因组可以在TAIR上下在,也可在ensemblgenomes下载。后者还可以下载其他植物的基因组

# 下载基因组
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# 构建核酸BLAST数据库
makeblastdb -in Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -dbtype nucl -out TAIR10 -parse_seqids

# 下载拟南芥protein数据
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
# 构建蛋白BLAST数据库
gzip -dArabidopsis_thaliana.TAIR10.pep.all.fa.gz
makeblastdb -in  Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR10 -parse_seqids

如果你从NCBI或者其他渠道下载了格式化过的数据库,那么可以用blastdbcmd去检索blast数据库,参数很多,常用就如下几个

  • db string : string表示数据库所在路径
  • dbtype string,: string在(guess, nucl, prot)中选择一个
  • 检索相关参数
    • -entry all 或 555, AC147927 或 gnl|dbname|tag'
    • -entry_batch 提供一个包含多个检索关键字的文件
    • -info 数据库基本信息
  • 输出格式 -outfmt %f %s %a %g ...默认是%f
  • out 输出文件
  • show_blastdb_search_path: blast检索数据库路径

使用案例

# 查看信息
blastdbcmd -db TAIR10 -dbtype nucl -info
# 所有数据
blastdbcmd -db TAIR10 -dbtype nucl -entry all | head
# 具体关键字,如GI号
blastdbcmd -db TAIR10 -dbtype nucl -entry 3 | head

还有其他有意思的参数,可以看帮助文件了解

可选:BLAST安装和更新nr和nt库

安装nt/nr库需要先进行环境变量配置,在家目录下新建一个.ncbirc配置文件,然后添加如下内容

; 开始配置BLAST
[BLAST]
; 声明BLAST数据库安装位置
BLASTDB=/home/xzg/Database/blast
; Specifies the data sources to use for automatic resolution
; for sequence identifiers
DATA_LOADERS=blastdb
; 蛋白序列数据库存本地位置
BLASTDB_PROT_DATA_LOADER=/home/xzg/Database/blast/nr
; 核酸数据库本地存放位置
BLASTDB_NUCL_DATA_LOADER=/home/xzg/Database/blast/nt
[WINDOW_MASKER]
WINDOW_MASKER_PATH=/home/xzg/Database/blast/windowmasker

配置好之后,使用BLAST+自带的update_blastdb.pl脚本下载nr和nt等库文件(不建议下载序列文件,一是因为后者文件更大,二是因为可以从库文件中提取序列blastdbcmd -db nr -dbtype prot -entry all -outfmt "%f" -out nr.fa ,最主要是建库需要花费很长时间),直接运行下列命令即可自动下载。

nohup time update_blastdb.pl nt nr > log &

如果你不像通过update_blastdb.pl下载nr和nt等库文件,也可以是从ncbi上直接下载一系列nt/nr.xx.tar.gz文件,然后解压缩即可,后续还可以用update_blastdb.pl进行数据更新。

报错: 使用update_blastdb.pl更新和下载数据库时候出现模块未安装的问题。解决方法,首先用conda安装对应的模块,然后修改update_blastdb.pl的第一行,即shebang部分,以conda的perl替换,或者按照如下方法执行。

perl `which update_blastdb.pl`

下载过程中请确保网络状态良好,否则会出现Downloading nt.00.tar.gz...Unable to close datastream报错。

第二步:选择blast工具

根据不同的需求,比如说你用的序列是氨基酸还是核苷酸,你要查找的数据是核甘酸还是氨基酸,选择合适的blast工具。不同需求的对应关系可以见下图(来自biostars handbook)

img_dfd8b295c8ba386848f74fa444a5fab9.png
BLAST工具

不同工具的应用范围虽然不同,但是基本参数都是一致的,学会blastn,基本上其他诸如blastp,blastx也都会了。

blastn的使用参数很多 blastn [-h] ,但是比较常用是如下几个

  • -db : 数据库在本地的位置,或者是NCBI上数据库的类型,如 -db nr
img_fc3a6df8a0146dd344c1d4303d516876.png
BLAST库
  • -query: 检索文件
  • -query_loc : 指定检索的位置
  • -strand: 搜索正义链还是反义链,还是都要
  • out : 输出文件
  • -remote: 可以用NCBI的远程数据库, 一般与 -db nr
  • -evalue 科学计数法,比如说1e3,定义期望值阈值。E值表明在随机的情况下,其它序列与目标序列相似度要大于这条显示的序列的可能性。 与S值有关,S值表示两序列的同源性,分值越高表明它们之间相似的程度越大

E值总结: 1.E值适合于有一定长度,而且复杂度不能太低的序列。2. 当E值小于10-5时,表明两序列有较高的同源性,而不是因为计算错误。3. 当E值小于10-6时,表时两序列的同源性非常高,几乎没有必要再做确认。

  • 联配计分相关参数: -gapopen,打开gap的代价;-gapextend, gap延伸的代价;-penalty:核算错配的惩罚; -reward, 核酸正确匹配的奖励;

  • 结果过滤:-perc_identity, 根据相似度

BLAST还提供一个task参数,感觉很有用的样子,好像会针对任务进行优化速度。

第三步:运行blast,调整输出格式。

我随机找了一段序列进行检索

echo '>test' > query.fa
echo TGAAAGCAAGAAGAGCGTTTGGTGGTTTCTTAACAAATCATTGCAACTCCACAAGGCGCCTGTAATAGACAGCTTGTGCATGGAACTTGGTCCACAGTGCCCTACCACTGATGATGTTGATATCGGAAAGTGGGTTGCAAAAGCTGTTGATTGTTTGGTGATGACGCTAACAATCAAGCTCCTCTGGT >> query.fa

用的是blastn 检索核酸数据库。最简单的用法就是提供数据库所在位置和需要检索的序列文件

blastn -db BLAST/TAIR10 -query query.fa
# 还可以指定检索序列的位置
blastn -db BLAST/TAIR10 -query query.fa  -query_loc 20-100
# 或者使用远程数据库
blastn -db nr -remote -query query.fa
blastn -db nt -remote -query query.fa

以上是默认输出,blast的-outfmt选项提供个性化的选择。一共有18个选择,默认是0。
0 = Pairwise,
1 = Query-anchored showing identities,
2 = Query-anchored no identities,
3 = Flat query-anchored showing identities,
4 = Flat query-anchored no identities,
5 = BLAST XML,
6 = Tabular,
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1),
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
其中6,7,10,17可以自定输出格式。默认是

qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore

简写 含义
qaccver 查询的AC版本(与此类似的还有qseqid,qgi,qacc,与序列命名有关)
saccver 目标的AC版本(于此类似的还有sseqid,sallseqid,sgi,sacc,sallacc,也是序列命名相关)
pident 完全匹配百分比 (响应的nident则是匹配数)
length 联配长度(另外slen表示查询序列总长度,qlen表示目标序列总长度)
mismatch 错配数目
gapopen gap的数目
qstart 查询序列起始
qstart 查询序列结束
sstart 目标序列起始
send 目标序列结束
evalue 期望值
bitscore Bit得分
score 原始得分
AC: accession

以格式7为实例进行输出,并且对在线数据库进行检索

blastn -task blastn -remote -db nr -query query.fa  -outfmt 7 -out query.txt
head -n 15 query.txt
img_843748ec1cb40c6d551b8439aa716a2f.png
BLAST BEST HIT

如果想输入序列,增加对应的格式qseq, sseq

blastn -query query.fa -db nr -outfmt ' 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore'

对已有序列进行注释时常见的best hit only模式命令行

blastn -query gene.fa -out gene.blast.txt -task megablast -db nt -num_threads 12 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt "7 std stitle" -perc_identity 50 -max_target_seqs 1
# 参数详解
-task megablast : 任务执行模式,可选有'blastn' 'blastn-short' 'dc-megablast' 'megablast' 'rmblastn'
-best_hit_score_edge 0.05 :  Best Hit 算法的边界值,取值范围为0到0.5,系统推荐0.1
-best_hit_overhang 0.25 : Best Hit 算法的阈值,取值范围为0到0.5,系统推荐0.1
-perc_identity 50 : 相似度大于50
-max_target_seqs 1 : 最多保留多少个联配

仅仅看参数依旧无用,还需要知道BLAST的Best-Hits的过滤算法。假设一个序列存在两个match结果,A和B,无论A还是B,他们的HSP(High-scoring Segment Pair, 没有gap时的最高联配得分)一定要高于best_hit_overhang,否则被过滤。如果满足下列条件则保留A

  • evalue(A) >= evalue (B)
  • score(A)/length(A) < (1.0-score_edge)*score(B)/length(B)

搭建网页BLAST

曾经的BLAST安装后提供wwwblast用于构建本地的BLAST网页工具,但是BLAST+没有提供这个工具,好在BLAST足够出名,也就有人给它开发网页版工具。如viroBLASTSequenceserver, 目前来看似乎后者更受人欢迎。

有root安装起来非常的容易

sudo apt install ruby ncbi-blast+ ruby-dev rubygems-integration npm
sudo gem install sequenceserver

数据库准备,前面步骤已经下载了拟南芥基因组的FASTA格式数据

sequenceserver -d /到/之前/建立/BLAST/文件路径

然后就可以打开浏览器输入IP:端口号使用了。

Sequenceserver高级用法

开机启动

新建一个/etc/systemd/system/sequenceserver.service文件,添加如下内容。注意修改ExecStart.

[Unit]
Description=SequenceServer server daemon
Documentation="file://sequenceserver --help" "http://sequenceserver.com/doc"
After=network.target

[Service]
Type=simple
User=seqservuser
ExecStart=/path/to/bin/sequenceserver -c /path/to/sequenceserver.conf
KillMode=process
Restart=on-failure
RestartSec=42s
RestartPreventExitStatus=255

[Install]
WantedBy=multi-user.target

然后重新加载systemctl

## let systemd know about changed files
sudo systemctl daemon-reload
## enable service for automatic start on boot
systemctl enable sequenceserver.service
## start service immediately
systemctl start sequenceserver.service

nginx反向代理:我承认没有基本的nginx的知识根本搞不定这一步,所以我建议组内使用就不要折腾这个。简单的说就是在nginx的配置文件的server部分添加如下内容即可。

location /sequenceserver/ {
    root /home/priyam/sequenceserver/public/dist;
    proxy_pass http://localhost:4567/;
    proxy_intercept_errors on;
    proxy_connect_timeout 8;
    proxy_read_timeout 180;
}

参考资料:http://www.sequenceserver.com/doc/

术语列表

引自BLAST Glossary

参考资料

目录
相关文章
|
算法 数据挖掘
群体遗传学研究荐读丨应知应会(下)
群体遗传学研究荐读丨应知应会(下)
|
算法 调度
详解操作系统四大常用的作业调度算法(FCFS丨SJF丨HRRN丨RR)
详解操作系统四大常用的作业调度算法(FCFS丨SJF丨HRRN丨RR)
8437 0
|
4月前
|
JSON NoSQL API
如何通过API获取贝壳找房二手房详情数据
本文详解贝壳找房二手房详情API调用:含RESTful接口结构、JSON响应解析、Python代码示例及关键实践要点(认证、限频、重试、Redis缓存),强调合规使用与数据安全,助力开发者高效构建房产应用。(239字)
595 1
|
运维 监控 测试技术
自动化运维实践:CI/CD流程详解
【6月更文挑战第30天】CI/CD实践推动软件开发自动化,通过持续集成确保代码质量,自动部署提升交付速度。核心流程包括:代码管理(Git等)、自动化构建与测试、代码审查、部署。关键点涉及选择工具、测试覆盖率、监控及团队协作。采用CI/CD能减少错误,但需应对挑战,如工具选型、全面测试和团队沟通。
5179 2
|
Linux Shell
Linux 进程前台后台切换与作业控制
进程前台/后台切换及作业控制简介: 在 Shell 中,启动的程序默认为前台进程,会占用终端直到执行完毕。例如,执行 `./shella.sh` 时,终端会被占用。为避免不便,可将命令放到后台运行,如 `./shella.sh &`,此时终端命令行立即返回,可继续输入其他命令。 常用作业控制命令: - `fg %1`:将后台作业切换到前台。 - `Ctrl + Z`:暂停前台作业并放到后台。 - `bg %1`:让暂停的后台作业继续执行。 - `kill %1`:终止后台作业。 优先级调整:
1356 5
|
Ubuntu Python
全网最简约的Vscode配置Anaconda环境(百分百成功)
全网最简约的Vscode配置Anaconda环境(百分百成功)
37503 0
全网最简约的Vscode配置Anaconda环境(百分百成功)
|
机器学习/深度学习 开发者 Python
Python中进行特征重要性分析的9个常用方法
在Python机器学习中,特征重要性分析是理解模型预测关键因素的重要步骤。本文介绍了九种常用方法:排列重要性、内置特征重要性(如`coef_`)、逐项删除法、相关性分析、递归特征消除(RFE)、LASSO回归、SHAP值、部分依赖图和互信息。这些方法适用于不同类型模型和场景,帮助识别关键特征,指导特征选择与模型解释。通过综合应用这些技术,可以提高模型的透明度和预测性能。
1859 0
|
机器学习/深度学习 人工智能 弹性计算
阿里云GPU服务器租用费用_GPU服务器详解_A10、V100、T4、P4、P100
阿里云GPU云服务器提供NVIDIA A10、V100、T4、P4、P100等多种GPU卡,适合深度学习、科学计算等场景。实例性能强劲,单实例可达5PFLOPS混合精度计算,VPC网络支持2400万PPS和160Gbps内网带宽。GPU实例包括A10卡GN7i(3213.99元/月起)、V100-16G卡GN6v(3830.00元/月起)等,价格因配置而异。阿里云还提供GPU加速软件如AIACC-Training和AIACC-Inference,以及弹性计算实例EAIS。客户案例包括深势科技、流利说和小牛翻译等。
|
安全 Shell Linux
探索Linux命令chsh:更改用户的默认shell
`chsh`是Linux命令,用于更改用户的默认登录shell。它涉及用户环境配置和系统安全,允许用户选择更适合自己的shell以提升效率。命令有交互式选项和参数如`-s`来指定新shell。在使用时要注意新shell的可执行性、权限问题及选择合适的shell。例如,要更改为bash,用户可运行`chsh`后按提示操作,而root用户能用`sudo chsh -s /bin/zsh john`为用户`john`设定zsh。在更改前,确认shell路径、权限,并了解不同shell的特点。
|
Linux 数据库 索引
如何在Linux中进行本地Blast序列比对?
如何在Linux中进行本地Blast序列比对?