全局比对与局部比对
2018-11-01 本文已影响0人
TOP生物信息
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上联网搜索呢？
