BLAST(Basic local alignment search tool)奄侠,即基于局部比對(duì)算法的搜索工具
BLAST包括blastn, blastp, blastx, tblastn, tblastx等. 使用conda安裝即可寂祥,亦可使用源碼(https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+)奄喂。
想安裝某個(gè)特定版本可以使用
conda install -c bioconda blast==版本號(hào)
# blast安裝perl模塊的方法
conda isntall perl-digest-md5
BLAST的基本步驟
- 用makeblastdb為BLAST建立數(shù)據(jù)庫(kù)
- 選擇BLAST工具赞别,blastn, blastp, blastx, tblastn, tblastx等
- 運(yùn)行BLAST工具,在有必要的情況下還可以對(duì)輸出結(jié)果進(jìn)行修飾以得到自己想要的結(jié)果格式
第一步 BLAST數(shù)據(jù)庫(kù)的構(gòu)建
格式化蛋白質(zhì)數(shù)據(jù)庫(kù):
makeblastdb -in spinach.fa -parse_seqids -hash_index -out spinach_DB -dbtype prot
格式化核酸數(shù)據(jù)庫(kù):
makeblastdb -in spinach.fa -parse_seqids -hash_index -out spinach_DB -dbtype nucl
#一般帶上-parse_seqids -hash_index
將基因組和想要構(gòu)建的數(shù)據(jù)庫(kù)名字替換即可峰档,一般我們做生物信息都會(huì)用本地BLAST寿谴,因此不介紹遠(yuǎn)程數(shù)據(jù)庫(kù)。
第二步 選擇BLAST工具
- blastn nucleotide to nucleotide
- blastp protein to protein
- blastx translated nucleotide to protein
- tblastn protein to translated nucleotide
blastp的常用參數(shù)
-query 后接查詢序列的文件名稱疏旨;
-db 后接格式化好的數(shù)據(jù)庫(kù)名稱很魂;
-out 后接要輸出的文件名稱及格式
blastn的常用參數(shù)
-db: 指定blast搜索用的數(shù)據(jù)庫(kù)
-query: 用來(lái)查詢的輸入序列,fasta格式
-out:輸出結(jié)果文件
-evalue: 設(shè)置e值cutoff
-max_target_seqs: 設(shè)置最多的目標(biāo)序列匹配數(shù)
-num_threads:指定多少個(gè)cpu運(yùn)行任務(wù)(依賴于你的系統(tǒng)檐涝,同于以前的-a參數(shù))
-outfmt format "7 qacc sacc evalue length pident" :這個(gè)是新BLAST+中最叼的功能了遏匆,直接控制輸出格式,不用再用parser啦谁榜, 7表示帶注釋行的tab格式的輸出幅聘,可以自定義要輸出哪些內(nèi)容,用空格分格跟在7的后面窃植,并把所有的輸出控制用雙引號(hào)括起來(lái)帝蒿,其中qacc查詢序列的acc,sacc表示目標(biāo)序列的acc,evalue即是e值巷怜,length即是匹配的長(zhǎng)度葛超,pident即是序列相同的百分比暴氏,其他可用的特征(代碼塊中顯示)如下:
Formatting options
-outfmt <String>
alignment view options:
0 = pairwise,
1 = query-anchored showing identities,
2 = query-anchored no identities,
3 = flat query-anchored, show identities,
4 = flat query-anchored, no identities,
5 = XML Blast output,
6 = tabular,
7 = tabular with comment lines,
8 = Text ASN.1,
9 = Binary ASN.1
10 = Comma-separated values
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
Options 6, 7, and 10 can be additionally configured to produce
a custom format specified by space delimited format specifiers.
一般來(lái)說(shuō)用的就是6、7绣张、10答渔、17,因?yàn)橐ㄟ^(guò)覆蓋度和得分來(lái)篩選結(jié)果的話用這三種比較方便侥涵,而7的話是比6多了一些注釋部分沼撕,結(jié)果部分是一樣的。
qseqid means Query Seq-id
qgi means Query GI
qacc means Query accesion
sseqid means Subject Seq-id
sallseqid means All subject Seq-id(s), separated by a ';'
sgi means Subject GI
sallgi means All subject GIs
sacc means Subject accession
sallacc means All subject accessions
qstart means Start of alignment in query #這個(gè)位置在最起始為1芜飘,而不是0务豺,注意這個(gè)位置和別的的計(jì)算未知的差別
qend means End of alignment in query
sstart means Start of alignment in subject #同理,這個(gè)也是最起始為1嗦明,不是0
send means End of alignment in subject
qseq means Aligned part of query sequence
sseq means Aligned part of subject sequence
evalue means Expect value
bitscore means Bit score
score means Raw score
length means Alignment length
pident means Percentage of identical matches
nident means Number of identical matches
mismatch means Number of mismatches
positive means Number of positive-scoring matches
gapopen means Number of gap openings
gaps means Total number of gaps
ppos means Percentage of positive-scoring matches
frames means Query and subject frames separated by a '/'
qframe means Query frame
sframe means Subject frame
第三步 運(yùn)行BLAST
blastp
blastp -query spinach_1.fa -db spinach_DB -out spinach_blast
blastn
自己經(jīng)常用的命令如下所示
blastn -task blastn -db datebase_name -query input_file -outfmt 6 -out output_file
blastn -task blastn -db database_name -query input_file -out output_file -evalue evalue -max_target_seqs num_sequences -num_threads int_value -outfmt format format_string
blastn -task blastn -db database_name -query input_file -out output_file -evalue evalue -max_target_seqs num_sequences -num_threads int_value -outfmt format " 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore"
一般默認(rèn)的使用輸出結(jié)果格式6的話輸出的結(jié)果文件的數(shù)據(jù)如下所示:
query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
注:第二條命令可以增加對(duì)應(yīng)的qseq, sseq