Blast

Blast有好几个工具，blastn,blastp,blastx还有tblastn，tblastx

blast 是干啥的

就是你有一条query的reads，他帮你在数据库里找个能配上的，是一个局部对比的工具，寻找两条序列最相似的地方

可以直接用ncbi的网页版，也可以自己搞一个本地的，看个人需求。NCBI的blast说明书

blast并不是最好的序列匹配工具，这个软件的调教偏向于高效率寻找可能的匹配位点，所以有时候不能给出全部可能匹配的位置。其目的是在一个很大的序列已知的数据库中找到匹配

blast的四条要点

• 可以找核苷酸序列，蛋白序列，也可以把核苷酸翻译成蛋白来找蛋白序列
• 可以指定不同的算法，返回的结果也很不同
• 其算法依赖于打分矩阵
• 可以自己进行微调，有很多参数一般用户从来都没有接触过

使用blast的几个基本步骤

• 先建一个数据库，makeblastdb
• 选一个工具帮你在数据库里找，可以是blastn也可以blastp，自己调一下参数
• 下边就让软件跑起来

blast小程序分别是干啥的

输入序列	数据库序列类型	匹配水平	应该叫啥	实际叫啥
核酸	核酸	核酸	blastNN	blastn
肽链	肽链	肽链	blastPP	blastp
核酸	肽链	肽链	blastNP	blastx
肽链	核酸	肽链	blastPN	tblastn
核酸	核酸	肽链	blastNNP	tblastx

blast的黑话

• Query: 你手头的短短小reads
• Target: 我们要去扒拉着找的大数据库
• Subject 我们叫能匹配上的输入的序列
• Score 打分
• E-value 期待值，说的是在一个老大的数据库里面扒拉着找到一个比目前的结果更好的匹配的可能性

咋建数据库

搞一个埃博拉病毒的的基因组序列，建立索引，这会搞出几个我们在别的地方根本用不到的文件，当然在这儿可是必须的。所以最好找个专门的地方放置这些索引文件，一般放到自己的目录的 ~/refs 里面

mkdir -p ~/refs/ebola
efetch -db nucleotide -id KM233118 --format fasta > ~/refs/ebola/KM233118.fa

这个里面有一个或者多个的fasta文件，接下来要把这个fasta搞成数据库，也就是为blast建一个索引，要用makeblastdb

makeblastdb -h
makeblastdb -help

makeblastdb -in ~refs/ebola/KM233118.fa -dbtype nucl -out ~/refs/ebola/KM233118 -parse_seqids

下边就可以来查询建好的数据库了，先建一条test序列

echo ">test" > query.fa
echo "AATCATACCTGGTTTGTTTCAGAGCCATATCACCAAGAT" >> query.fa

跑一下blastn，其实这个过程就是我们在网页粘上一条序列，点一下blast按钮是一样的。

blastn -db ~/refs/ebola/KM233118 -query query.fa

按照自己的需要来设定输出格式

可以改成输出格式6或者7（带评论和header的表格）

  qseqid means Query Seq-id
      qgi means Query GI
        qacc means Query accesion
     qaccver means Query accesion.version
        qlen means Query sequence length
      sseqid means Subject Seq-id
   sallseqid means All subject Seq-id(s), separated by a ';'
          sgi means Subject GI
       sallgi means All subject GIs

blastn -db ~/refs/ebola/KM233118 -query query.fa  -outfmt "6 qseqid sseqid pident"

这样信息就很简约

test    gi|667853353|gb|KM233118.1|    100.000

blast task

task就是对算法的优化和修改，让blast工具使用起来更加麻溜完成特定的工作
但是task本身有时候看起来就跟一个程序差不多
比如
blastn 找出更多的不同的align
mega-blast 找少一些不同algin
blastn-short 短小序列的query

这帮子人把这个task搞得比这个还复杂，比如blastn的默认task是megablast

比如咱搞一段短序列，跑一下blast，啥都找不到，但是咱明明就是从开头部分选的这段一样的序列哈

echo ">short" > short.fa
echo "AATCATACCTGG" >> short.fa

blastn -db ~/refs/ebola/KM233118 -query short.fa

要想找着，就得改参数

blastn -db ~/refs/ebola/KM233118 -query short.fa -task blastn

这样就搞定

咱们要是把序列搞得再短一点，那么就只有blastn-short可以搞定了

echo ">mini" > mini.fa
echo "AATCATA" >> mini.fa

blastn -db ~/refs/ebola/KM233118 -query mini.fa -task blastn-short

就是这么的任性，你能肿么办

两条短序列的对比

这就不用搞啥数据库了，直接来就行

blastn -query query.fa -subject ~/refs/ebola/KM233118.fa 

blast能不能把所有可能的位置都找到

blast是给大众干活的，要求太高对他来说就有点不公平，比如基因组中，blast会自动过滤掉低复杂度的区域

#下载yeast的第一个染色体
efetch  -id NC_001133  -db nucleotide -format fasta > NC_001133.fa  
head -2 NC_001133.fa > start.fa 

#但是你blast回整个基因组的时候，发现找不到match的微店
blastn -query start.fa -subject NC_001133.fa

#这是因为query的序列有很多低复杂度的区域，自动过滤掉了
CCACACCACACCCACACACCCACACACCACACCACACACCACACCACACCCACACACACACATCCTAACA

#绝大多数的时候，这是个好事，但是这儿就不是了，所以要关掉这个自动过滤用 -dust no
blastn -query start.fa -subject NC_001133.fa -dust no

建立blast数据库

其实就是创建一个索引，这个索引是必须的，主要是提高搜索的效率，否则一个一个的扒拉着去找就太费劲了。对于核酸和肽链来说，创建的索引是不一样的。

makeblastdb 创建blast数据库
blastdbcmd 搜索数据库
update_blastdb.pl 更新已有的数据库

我们可不可以直接下载blast数据库？

虽然咱们可以用makeblastdb来创建数据库索引，但是费老大劲，其实ncbi有现成的可以直接下载下来
NCBI的blast的索引

我们可不可以自动化这个下载索引的过程

#先找一个地方放我们的下载下来的数据库
mkdir -p ~/refs/refseq
cd ~/refs/refseq/
##先看看有哪一些数据库可以下载
update_blastdb.pl --showall
##下载16微生物的索引文件
update_blastdb.pl 16SMicrobial --decompress
##下载分类学数据库，这个以后有用
update_blastdb.pl taxdb --decompress

我们也可以把blast的索引路径加入到环境变量中，这样以后就不用一次次的指定了。

export BLASTDB=$BLASTDB:~/refs/refseq/
#上述变量指定之后，我们就可以直接这样打命令了
blastdbcmd -db 16SMicrobial -info

查看已经建立好的数据库的信息

blastdb