Biostar: 课程11、12

翻译：生物女博士

注：本系列课程翻译已获得原作者授权。

#为作者的注释

#*为译者的注释

#* 本文需要用到NCBI的工具，详情参见课程3、4

Lecture 11 - 安装和使用Blast

# 获取并安装BLAST+

# 下边的URLs 是MAC系统的，Linux的需自行修改。

# 转到src文件夹

cd ~/src

# 获取Mac的blast可执行版本。

#*译者已将原版本下载链接更改为最新版本的blast

curl -O  ftp://ftp.ncbi.nlm./blast/executables/blast+/LATEST/ncbi-blast-2.6.0+-x64-macosx.tar.gz

tar zxvf ncbi-blast-2.6.0+-x64-macosx.tar.gz

# Linux 则下载另一个包。

# curl -O ftp://ftp.ncbi.nlm./blast/executables/blast+/LATEST/ncbi-blast-2.6.0+-x64-linux.tar.gz

# 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

# http://www.ncbi.nlm./genome/4887

# 这个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.%2Fbooks%2FNBK279690%2Fpdf%2FBookshelf_NBK279690.pdf