搜索
查看: 3023|回复: 0

[Other] Biostar: 课程11、12-下载序列,建blast数据库以及本地blast

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-7-20 10:40:09 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:http://mp.weixin.qq.com/s/q7yT0Z-8d2nD8bjleJHpsw


翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

#为作者的注释
#*为译者的注释





Lecture 11 - 安装和使用Blast

# 获取并安装BLAST+
# 下边的URLs 是MAC系统的,Linux的需自行修改。
# 转到src文件夹
cd ~/src

# 获取Mac的blast可执行版本。
#*译者已将原版本下载链接更改为最新版本的blast
tar zxvf ncbi-blast-2.6.0+-x64-macosx.tar.gz

# Linux 则下载另一个包。
# tar zxvf ncbi-blast-2.6.0+-x64-linux.tar.gz

# 我可以把可执行的软件或者整个文件夹加到我的路径下。通常我都是添加整个文件夹。
#*Mac下
echo 'export PATH=$PATH:~/src/ncbi-blast-2.6.0+/bin' >> ~/.profile

# Linux下
# echo 'export PATH=$PATH:~/src/ncbi-blast-2.6.0+/bin' >> ~/.bashrc

# 加载新设定到当前终端(新开的终端会自动加载)
#*这是mac下的。
source ~/.profile

#*对应Linux下
source ~/.bashrc

# 检测是否可用了。简单版本的帮助文档。
blastn -h

# 完整的帮助文档。
blastn -help

# 检测清单:
#   1. 我们想找什么?-> 查询序列
#   2. 到哪里去找? -> 数据库 -> 目标序列
#   3. 如何找到? -> 搜索类型

# 我们来建一个blast的数据库
# 获取埃博拉病毒的基因组序列
# 当我们索引数据库是,可能会产生多个文件。最好还是把这些文件单独放到一个文件夹里。我们给这文件夹命名为refs。
mkdir -p refs
efetch -db nucleotide -id KM233118 --format fasta > refs/KM233118.fa

# 参考序列可能有一个或多个fastq记录。要把它变成一个blast的数据库,我们需要把它转成程序可以操作的格式。
# blast数据库建立索引的完整或者简单的帮助文档
makeblastdb -h
makeblastdb -help

# 从基因组创建一个核苷酸数据库。
makeblastdb -in refs/KM233118.fa -dbtype nucl

# 来看一下发生了什么。
ls refs

# 创建查询序列
head -2 refs/KM233118.fa > query.fa

# 搜索这条序列:核苷酸序列 vs 核苷酸数据库
# 注意一下格式。
blastn -db refs/KM233118.fa -query query.fa | more

# 我们可以修改输出格式。 当碱基对碱基的比对信息不需要时,我们可以采用表格格式。
blastn -db refs/KM233118.fa -query query.fa -outfmt 6

# 这个同时会显示表头信息。
blastn -db refs/KM233118.fa -query query.fa -outfmt 7

# 你也可以自己定义输出格式。
blastn -db refs/KM233118.fa -query query.fa -outfmt '6 qseqid sseqid qlen length mismatch'

#*在blastn -help里输出的帮助文档内有详细的格式说明。分别告诉你不同的数字表示的输出格式,以及自定义格式时,每个参数的含义。
# 默认情况下将采用MEGABLAST工具搜索.
# 将查询序列缩短为:
# >short
# AATCATACCTGG

#*你可以这么生成一个新的query.fa:
echo ">short" > query.fa
echo "AATCATACCTGG" >> query.fa

# blast有不同的搜索策略或task。
#*task里的参数有'blastn','blastn-short' ,'dc-megablast','megablast', 'rmblastn',默认的是'megablast'
#*使用不同的task,出来的结果内会有文献引用,有兴趣的朋友可以去看一下。
#*blastn:需要11bp以上的精确匹配。
#*blast-short:对于小于50bp的序列最优。
#*megablast:一般用于非常相似的序列。
#*dc-megablast:dc(discontiguous,不连续的),一般用于找距离远(相似性低的,比如种间同源)的序列。

# 这些是一些预设参数集的典型。
blastn -db refs/KM233118.fa -query query.fa
blastn -db refs/KM233118.fa -query query.fa -task blastn

# 把序列弄得更短些。
# >short
# AATCATA
# 现在 -task blastn 也不能用了, 有另一个task叫blastn-short
# Blast的参数调整本身就是个小世界。
blastn -db refs/KM233118.fa -query query.fa -task blastn-short

# 你在目标数据库中,可以有多条序列。
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta > refs/all-genomes.fa

# 对所有的基因组创建blast数据库。
makeblastdb -in refs/all-genomes.fa -dbtype nucl

# 随意挑取基因组的任意区域。
efetch -db nucleotide -id KM233118 -format fasta -seq_start 1 -seq_stop 1000 > 1000.fa

#*大家可以试一下
#*blastn -db refs/all-genomes.fa -query 1000.fa -outfmt 6



Lecture 12 - blastdbcmd, blastp, blastx, tblastn 的用法

# blastdbcmd, blastp, blastx, tblastn 的用法
# 比对1999年爆发的新病毒。
# RefSeq 访问号 NC_002549.1, 核蛋白例子: NP_066243.1
# 这个Bioproject 数据来自1999年。
# 获取特征表格
esearch -db protein -query PRJNA14703 | efetch -format ft

# 搜索基因组的特征, ft -> feature table(特征表格)的简写。
efetch -db nucleotide -id NC_002549 -format ft | more

# 从蛋白质数据库里获取这个项目的数据
esearch -db protein -query PRJNA14703 | efetch -format fasta > NC_002549-prot.fa

# 添加parse eqids选项来makeblastdb.
# 新建一个核苷酸和蛋白质的Blast数据库。

#*NC_002549的核苷酸序列,前边并没有下载下来。
#*请先下载:esearch -db nucleotide -query NC_002549 | efetch -format fasta > NC_002549.fa
#*这里面由于前面下载的时候路径没有改到refs/下,所以下面运行会出现文件的路径问题的。
#*我们现在先把NC_002549.fa和 NC_002549-prot.fa移到refs/下
#* mv NC_002549* refs/
#*好了,可以开始了。
makeblastdb -in refs/NC_002549.fa -dbtype nucl -parse_seqids
makeblastdb -in refs/NC_002549-prot.fa -dbtype prot -parse_seqids

# blastdbcmd 有很多参数,可以对Blast数据库内容进行查询和格式设定。
# 列出数据库的所有内容。
blastdbcmd -db refs/NC_002549.fa -entry 'all' | head -3

# 获取数据库的特定条目。
blastdbcmd -db refs/NC_002549.fa -entry '10313991' | head -3
#*报错,暂未解决。

# 获取数据库中的一定范围的核苷酸条目。
blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %s' -range 1-10 -strand minus
#*%a表示accession(访问号),%s 表示 sequence data(序列数据)

# 设定数据库内容的格式。
blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %l'
#*%l 表示 sequence length(序列长度)

# 列出每个蛋白及其长度。
blastdbcmd -db refs/NC_002549-prot.fa -entry 'all' -outfmt '%a %l'

# 从蛋白数据库中提取一个条目并存到一个文件中。
blastdbcmd -db refs/NC_002549-prot.fa -entry 'NP_066243.1' > query-p.fa

# 运行blastp。
#*blastp是蛋白序列比对蛋白数据库。
blastp -db refs/NC_002549-prot.fa -query query-p.fa

# 运行 tblastn。
# 提取蛋白质编码区域。我是自己从GenBank页面去找对应蛋白区域的。
#*自行到NCBI搜索NC_002549。NCBI会给出cds区域信息的
blastdbcmd -db refs/NC_002549.fa -entry 'NC_002549' -range 470-2689  > nucleotide.fa

# 将核苷酸序列比对到蛋白数据库中。
#*blastx是核苷酸序列比对到蛋白数据库
blastx -db refs/NC_002549-prot.fa -query nucleotide.fa | more

# 抓取核蛋白。
#*tblastn是蛋白序列比对到核苷酸数据库。
efetch -db protein -id NP_066243.1 -format fasta > NP_066243.fa
tblastn -db refs/NC_002549.fa -query NP_066243.fa  | more

# 把2014的埃博拉项目的所有蛋白获取下来。
esearch -db protein -query PRJNA257197 | efetch -format fasta > refs/ebola-2014.fa
makeblastdb -in refs/ebola-2014.fa -dbtye prot -parse_seqids

# 抓取核蛋白。
efetch -db protein -id NP_066243.1 -format fasta > NP_066243.fa

#*大家可以试试对这条蛋白进行序列比对。


参考资料(pdf文件):chrome-extension://padcapdkhelngdelppbbjmkmkfceoikg/content/pdf/viewer/viewer.html?file=https%3A%2F%2Fwww.ncbi.nlm.nih.gov%2Fbooks%2FNBK279690%2Fpdf%2FBookshelf_NBK279690.pdf


欢迎到微信公众号订阅我们
生信媛
bio_sxy

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:我为什么要写博客---生信笔记
下一篇:Biostar:课程13、14-常用的mapping软件bwa & sam格式简介
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:32 , Processed in 0.041034 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.