biostar handbook(七)| BLAST

2024-06-23 21:18
文章标签 blast biostar handbook

本文主要是介绍biostar handbook(七)| BLAST,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

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
USAGEmakeblastdb [-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)

2013053-a10280e8ef6d389b.jpg
blast类型

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

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

  • -db : 数据库在本地的位置,或者是NCBI上数据库的类型,如 -db nr
2013053-d6e43f4c62a60bc4.jpg
image
  • -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错配数目
gapopengap的数目
qstart查询序列起始
qstart查询序列结束
sstart目标序列起始
send目标序列结束
evalue期望值
bitscoreBit得分
score原始得分

AC: accession

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

blastn -task blastn -remote -db nr -query query.fa  -outfmt 7 -out query.txt
head -n 15 query.txt
2013053-eb6ae1ea9f4e846e.jpg
image

术语列表

引自BLAST Glossary

参考资料

  • Linux系统中NCBI BLAST+本地化教程
  • BLAST Command Line Applications User Manual

这篇关于biostar handbook(七)| BLAST的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1088304

相关文章

biostar handbook(六)|正则表达式和K-mers在模式匹配中的使用

模式匹配中的正则表达式和k-mer 模式匹配指的是在看似杂乱无章的系统中找到符合要求的部分。比如说你想从基因组中寻找motif,转录因子结合位点,CDS, 或者检测测序结果里是否有接头等,这些行为都可以解读为根据已有的模式去寻找目标序列中符合要求的片段。后面介绍的序列模式匹配的两种方法,正则表达式和k-mers。 正则表达式的基础概念 正则表达式(regular expression)的概念

biostar handbook: 第五周笔记汇总+第六周任务布置

简单总结下第五周的笔记: biostarhandbook(五)|序列从何而来和质量控制YXF-测序仪和质量控制Biostar学习笔记(5)测序平台、测序原理及质量控制Biostar第六课 测序仪和QC 顺便发布第六周的任务: 第六周的任务是第十章和第十一章。其中第十一章介绍的是序列模式,是对正则表达式的温故知新。第十一章则是开始序列分析的重要一步:序列联配。任务如下: 什么是正则表达式什么

biostar handbook(五)|序列从何而来和质量控制

测序仪 2017年一篇发表在Nature的综述"DNA sequencing at 40: past, present and future"介绍了DNA测序这40年的发展历程。1976年,Sanger和Coulson同时发表了2种方法用于对上百个DNA碱基进行解码,这就是第一代测序技术。到了2005年,罗氏的454平台揭开了高通量测序的序幕,后面则是SOLiD,454和Illumina三方对抗

biostar handbook: 第四周笔记汇总+第五周任务布置

不知不觉已经过去了四周,这个系列的开篇语写于2017年10月14日,距离今天差不多是一个月的时间了。这个月的时间学的内容并不算多,大致也就是如下几个内容: *nux基础: 这个尤为重要,学会了*nix(Linux或unix)之后,如果能够在日常科研生活中进行使用,那么你的效率将会大大提高。数据格式: 数据的保存具有一定的格式,处理数据的前提在于知己知彼,这样才能选择合适的工具。这个部分内容包括知

biostar handbook(四)|生物数据及其下载和基本操作

2017/11/9 第一版: 生物数据库,基本数据类型(genbank, fasta/fastq),数据上传站点 2017/11/12 第二版:如何利用esearch, efecth快速获取SRR序列号 生物数据库 目前绝大部分数据由NCBI, EMBL-EBI, DDBJ三大机构托管,可划分为五类: (表格数据来源于INSDC) Data typeDDBJEMBL-EBINCBIN

biostar handbook: 第三周笔记汇总+第四周任务布置

就目前来看,这一周主要以补交上次作业为主,所以上一周的总结增加了很多新内容。但是关于本体论和富集分析,大家还是有点小困难。当然我自己交东西也慢了很多,因为时间也总是不太够,完全靠周末挤出来。 第三周笔记汇总如下 Biostar学习笔记(3)Gene set analysis related topics.本体论biostarhandboo(四)|本体论和功能分析基因本体论 第四周作业发布

biostar handbook: 第二周笔记汇总+第三周任务布置

第二周已经结束了,我不确定大家对Linux到底了解到了什么程度,但是我觉得再给一周时间让初学者去熟悉Linux肯定是必要的。于是这一周的任务不会太难,只需要让大家去理解本体论(ontology)。 笔记汇总 这周有一些小伙伴开始遭遇人生抉择,有一些则是要出差赶路,所以上交作业不算太多。可能大家对自己的要求有点高了,其实我一直强调的是笔记的不断迭代,只要你写了一点内容就可以发出来,后来不断修改,

biostar handbook(十)|如何进行变异检测

变异检测流程 什么是基因组变异 基因组变异是一个定义比较模糊的概念. 所谓的变异是相对于一个完美的“参考基因组”而言。但是其实完美的“参考基因组”并不存在,因为我们只是选择某一个物种里的其中似乎比较正常的个体进行测序组装,然后基于它进行后续的研究。简单的说,参考和变异是相对而言,变异也可能完全正常。 常见的基因组变异一般可以归为如下几类: SNP, 单核苷酸多态性, 一

biostar handbook|如何模拟NGS测序结果

如何用软件模拟NGS数据 为了评价一个工具的性能,通常我们都需要先模拟一批数据。这样相当于有了参考答案,才能检查工具的实际表现情况。因此对于我们而言,面对一个新的功能,可以先用模拟的数据测试下不同工具的优缺点。有如下几个工具值得推荐一下: 'wgsim/dwgsim': 从全基因组中获取测序reads'msbar': EMBOSS其中一个工具,能够从单个序列中模拟随机突变'biosed': E

「JCVI教程」如何基于物种的CDS的blast结果绘制点图(dotplot)

这是唐海宝老师GitHub上的JCVI工具的非官方说明书。 该工具集的功能非常多,但是教程资料目前看起来并不多,因此为了能让更多人用上那么好用的工具,我就一边探索,一边写教程 这一篇文章教大家如何利用JCVI里面的工具绘制点图,展现两个物种之间的共线性关系。 在分析之前,你需要从PhytozomeV11 下载A.thaliana和Alyrata的CDS序列,保证文件夹里有如下内容 Al