全局比对与局部比对

Needleman-Wunsch算法和Smith-Waterman算法分别对应全局比对和局部比对，这两个生物序列比对算法比较简单，这里就不介绍了。

使用MUSCLE进行全局比对

wget https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
$ tar zxvf muscle3.8.31_i86linux64.tar.gz

Basic usage

muscle -in <inputfile> -out <outputfile>

Common options

    -in <inputfile>    Input file in FASTA format (default stdin)
    -out <outputfile>  Output alignment in FASTA format (default stdout)
    -diags             Find diagonals (faster for similar sequences)
    -maxiters <n>      Maximum number of iterations (integer, default 16)
    -maxhours <h>      Maximum time to iterate in hours (default no limit)
    -html              Write output in HTML format (default FASTA)
    -msf               Write output in GCG MSF format (default FASTA)
    -clw               Write output in CLUSTALW format (default FASTA)
    -clwstrict         As -clw, with 'CLUSTAL W (1.81)' header
    -log[a] <logfile>  Log to file (append if -loga, overwrite if -log)
    -quiet             Do not write progress messages to stderr
    -version           Display version information and exit

使用举例

查看待查询文件
$ lsx msa.fasta #in ~/learn_alignment/MUSCLE/data1
>sp|Q9BXR5|HUMAN_TLR10_TIR
CLHFDLPWYLRMLGQCTQTWHRVRKTTQEQLKRNVRFHAFISYSEHDSLWVKNELIPNLE
KEDGSILICLYESYFDPGKSISENIVSFIEKSYKSIFVLSPNFVQNEWCHYEFYFAHHNL
FHENSDHIILILLEPIPFYCIPTRYHKLKALLEKKAYLEWPKDRRKCGLFWANLRAAINV
NVLATREMYELQTFTELNEESRGSTISLMRTDCL
>sp|Q9NR96|HUMAN_TLR9_TIR
GWDLWYCFHLCLAWLPWRGRQSGRDEDALPYDAFVVFDKTQSAVADWVYNELRGQLEECR
GRWALRLCLEERDWLPGKTLFENLWASVYGSRKTLFVLAHTDRVSGLLRASFLLAQQRLL
EDRKDVVVLVILSPDGRRSRYVRLRQRLCRQSVLLWPHQPSGQRSFWAQLGMALTRDNHH
FYNRNFCQGPTAE
...

查询
$ muscle3.8.31_i86linux64 -in msa.fasta -out msa_out.fasta -clw -quiet
查看结果文件
$ lsx msa_out.fasta
MUSCLE (3.8) multiple sequence alignment
sp|O60602|HUMAN_TLR5_TIR       --TKFRGFCFICYKTA-QRLVFKDHPQGTE------PDMYKYDAYLCFSSKDF---TWVQ
sp|O15455|HUMAN_TLR3_TIR       ---EGWRISFY-WNVSVHRV--LGFKEIDRQ-----TEQFEYAAYIIHAYKDK---DWVW
sp|Q9NR96|HUMAN_TLR9_TIR       ----GWDL-WYCFHLCLAWLPWRGRQSGRD------EDALPYDAFVVFDKTQSAVADWVY
sp|Q9NR97|HUMAN_TLR8_TIR       HHLFYWDV-WFIYNVCLAKV--KGYRSLS-------TSQTFYDAYISYDTKDASVTDWVI
sp|Q9NYK1|HUMAN_TLR7_TIR       -HLYFWDV-WYIYHFCKAKI--KGYQRLI-------SPDCCYDAFIVYDTKDPAVTEWVL
sp|O00206|HUMAN_TLR4_TIR       --KFY----FHLMLLA-GCI--KYGR----------GENI-YDAFVIYSSQDE---DWVR
sp|O60603|HUMAN_TLR2_TIR       --HRFHGL-WYMKMMW-AWL--QAKRKPRKAP----SRNICYDAFVSYSERDA---YWVE
sp|Q9BXR5|HUMAN_TLR10_TIR      --CLHFDLPWYLRMLG-QCT--QTWHRVRKTTQEQLKRNVRFHAFISYSEHDS---LWVK
sp|Q9Y2C9|HUMAN_TLR6_TIR       ----YLDLPWYLRMVC-QWT--QTRRRARNIPLEELQRNLQFHAFISYSEHDS---AWVK
sp|Q15399|HUMAN_TLR1_TIR       ---SYLDLPWYLRMVC-QWT--QTRRRARNIPLEELQRNLQFHAFISYSGHDS---FWVK
                                        :                               : *:: .   :     **

不加-clw就生成FASTA，这个格式可以用于构建进化树。
当然这个软件也有在线版的。

使用BLAST系列进行局部比对

$ wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.7.1+-x64-linux.tar.gz
$ tar zxvf ncbi-blast-2.7.1+-x64-linux.tar.gz
$ cd ncbi-blast-2.7.1+/bin/
$ ls
blastdb_aliastool*  blastn*             deltablast*       makembindex*    rpstblastn*  update_blastdb.pl*
blastdbcheck*       blastp*             dustmasker*       makeprofiledb*  segmasker*   windowmasker*
blastdbcmd*         blastx*             legacy_blast.pl*  psiblast*       tblastn*
blast_formatter*    convert2blastmask*  makeblastdb*      rpsblast*       tblastx*

BLAST，图片来自山东大学《生物信息学》

添加到环境变量中

$ vi ~/.bashrc
...
export PATH="/ifs1/Grp3/huangsiyuan/learn_alignment/BLAST/ncbi-blast-2.7.1+/bin:$PATH"
...
$ source ~/.bashrc

创建一个目录来存储将使用的序列数据库

这是BLAST设置的最后一步

$ mkdir database
$ vi ~/.bashrc
export BLASTDB="/ifs1/Grp3/huangsiyuan/learn_alignment/BLAST/database/"
$ source ~/.bashrc
同样需要添加到~/.bashrc中,要用完整的绝对路径（最后的那个正斜线也不能少），此时BLASTDB就是一个变量了。

从NCBI下载序列数据库

在刚才下载BLAST的FTP站点上，可以找到很多预先准备好的数据库，见precomputed databases。关于每一个子数据库的介绍见同站点的文档。

$ cd $BLASTDB
$ nohup wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/refseq_protein.45.tar.gz &
因为文件基本上都很大，下载又很慢，所以就挑了一个较小的文件来下载
$ tar zxvf refseq_protein.45.tar.gz
$ rm -f refseq_protein.45.tar.gz
$ ls
refseq_protein.45.phr  refseq_protein.45.pni  refseq_protein.45.ppi  refseq_protein.45.psq  taxdb.bti
refseq_protein.45.pin  refseq_protein.45.pog  refseq_protein.45.psd  refseq_protein.pal
refseq_protein.45.pnd  refseq_protein.45.ppd  refseq_protein.45.psi  taxdb.btd

这些文件是BLAST执行查询所需的索引和元数据。

BLAST查询

#blastp --help查看帮助文档
$ blastp -query test.fa -db refseq_protein.45 -out results.txt
$ lsx results.txt
Query= sp|Q9BXR5|HUMAN_TLR10_TIR

Length=214
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

NP_001272534.1  toll-like receptor 1 [Capra hircus]                   259     1e-80
NP_001272469.1  toll-like receptor 6 [Capra hircus]                   233     1e-70
XP_003221759.1  PREDICTED: toll-like receptor 2 [Anolis carolinen...  157     7e-43
XP_003221760.1  PREDICTED: toll-like receptor 2 type-2 isoform X2...  154     7e-42
XP_003230533.1  PREDICTED: toll-like receptor 2 [Anolis carolinen...  143     1e-37
......

有一个问题，上述这种方法准备序列数据库文件不仅要一个一个下载还很占空间，有没有什么方法可以在Linux上联网搜索呢？

---

reference

https://www.exoscale.com/syslog/blast/