基本概念
相似性（similarity）
- 一種很直接的數(shù)量關(guān)系糠馆，比如部分相同或相似的百分比或其他一些合適的度量
- 如：A序列和B序列的相似性是80%
同源性（homology）
- 從一些數(shù)據(jù)中推斷出的兩個(gè)基因或者蛋白序列具有共同祖先的結(jié)論，屬于質(zhì)的判斷
- 可以說A序列和B序列是同源序列怎憋，但不能說同源性80%
常用工具
- BLAST
- BLAT
序列比對(duì)的常用工具：BLAST榨惠，但是其運(yùn)行速度慢的令人捉急。
一、BLAST（Basic Local Alignment Search Tool赠橙，局部相似性基本查詢工具）
BLAST（Basic Local Alignment Search Tool）是一套在蛋白質(zhì)數(shù)據(jù)庫或DNA數(shù)據(jù)庫中進(jìn)行相似性比較的分析工具耽装。BLAST程序能迅速與公開數(shù)據(jù)庫進(jìn)行相似性序列比較。BLAST結(jié)果中的得分是對(duì)一種對(duì)相似性的統(tǒng)計(jì)說明期揪。
https://www.biomart.cn/experiment/599/608/19912_0.htm
-F Filter query sequence (DUST with blastn, SEG with others) [String] default = T
查詢序列過濾掉奄，將那些 給出影響比對(duì)結(jié)果的低復(fù)雜度區(qū)域過濾掉。用blastn進(jìn)行查詢的序列用DUST程序過濾凤薛，其他的用SEG過濾 姓建。對(duì)DUST和SEG的詳細(xì)情況，用戶可以自己查詢資料缤苫。
使用此參數(shù) -F F 即可獲得完整的比對(duì)
在對(duì)核苷酸或蛋白質(zhì)序列數(shù)據(jù)庫進(jìn)行Blast搜索之前速兔，必須要對(duì)所使用的序列數(shù)據(jù)庫進(jìn)行formatdb， 即對(duì)序列數(shù) 據(jù)庫進(jìn)行格式化活玲，這是所有使用BLAST所必須的一步涣狗。
1.格式化序列數(shù)據(jù)庫— —formatdb
formatdb 簡(jiǎn)單介紹：formatdb處理的都是格式為 ASN.1和 FASTA，而且不論是核苷酸序列數(shù)據(jù)庫舒憾，還是蛋白質(zhì)序列數(shù)據(jù)庫镀钓；不論是使用Blastall ，還是Blastpgp镀迂，Mega Blast應(yīng)用程序丁溅，這一步都是不可少的。
主要參數(shù)的說明
-i 輸入需要格式化的源數(shù)據(jù)庫名稱 Optional
-p 文件類型探遵，是核苷酸序列數(shù)據(jù)庫窟赏，還是蛋白質(zhì)序列數(shù)據(jù)庫 T：protein , F：nucleotide
-n 重命名數(shù)據(jù)庫文件的名稱 ；
-t 數(shù)據(jù)庫的標(biāo)題【可選】箱季；
-l 日志文件名涯穷，formatdb.log
-o 解析選項(xiàng). T - True: 解析序列標(biāo)識(shí)并且建立目錄,F - False: 與上相反,[T/F] Optional default = F
示例：
formatdb -i uniref100.fasta -n uniref100 -t uniref100 -l uniref100.log -p T
formatdb -i uniref90.fasta -n uniref90 -t uniref90 -l uniref90.log -p T
formatdb -i uniref50.fasta -n uniref50 -t uniref50 -l uniref50.log -p T
2.運(yùn)行blastall
blastall -i query.fasta -d uniref50 -o blast.out -p blastn -F F -e 1e-5 -m 8
1．-e參數(shù)
用來過濾比對(duì)較差的結(jié)果的，用"-e"參數(shù)指定一個(gè)實(shí)數(shù)规哪，blast 會(huì)過濾掉期望值大于這個(gè)數(shù)的比對(duì)結(jié)果求豫。這樣不但簡(jiǎn)化了結(jié)果，還縮短了運(yùn)行時(shí)間和結(jié)果占用的空間诉稍。
2. -p參數(shù)
-p blastp：蛋白序列與蛋白庫做比對(duì)蝠嘉。
-p blastx：核酸序列對(duì)蛋白庫的比對(duì)。
-p blastn：核酸序列對(duì)核酸庫的比對(duì)杯巨。
-p tblastn：蛋白序列對(duì)核酸庫的比對(duì)蚤告。
3.-F 參數(shù)
-F(T/F)參數(shù)是用來屏蔽簡(jiǎn)單重復(fù)和低復(fù)雜度序列的。如果選“T”服爷，程序在比對(duì)過程中會(huì)屏蔽query中的簡(jiǎn)單重復(fù)和低復(fù)雜度序列杜恰；選"F"則不會(huì)屏蔽获诈。缺省值為"T"。
4. -m
m8格式如下圖心褐，12列：
1舔涎、Query id：查詢序列ID標(biāo)識(shí)
2、Subject id：比對(duì)上的目標(biāo)序列ID標(biāo)識(shí)
3逗爹、% identity：序列比對(duì)的一致性百分比
4亡嫌、alignment length：符合比對(duì)的比對(duì)區(qū)域的長度
5、mismatches：比對(duì)區(qū)域的錯(cuò)配數(shù)
6掘而、gap openings：比對(duì)區(qū)域的gap數(shù)目
7挟冠、q. start：比對(duì)區(qū)域在查詢序列(Query id)上的起始位點(diǎn)
8、q. end：比對(duì)區(qū)域在查詢序列(Query id)上的終止位點(diǎn)
9袍睡、s. start：比對(duì)區(qū)域在目標(biāo)序列(Subject id)上的起始位點(diǎn)
10知染、s. end：比對(duì)區(qū)域在目標(biāo)序列(Subject id)上的終止位點(diǎn)
11、e-value：比對(duì)結(jié)果的期望值斑胜，解釋是大概多少次隨機(jī)比對(duì)才能出現(xiàn)一次這個(gè)score,Evalue越小控淡，表明這種情況從概率上越不可能發(fā)生，那么這個(gè)比對(duì)的可靠性越高伪窖。
12逸寓、bit score：比對(duì)結(jié)果的bit score值
一般情況我們看第3居兆、11覆山、12兩列，e值越小越可靠泥栖。
blastall 老版本對(duì)應(yīng)的參數(shù)是 -m 8
blast+對(duì)應(yīng)的參數(shù)是-outfmt 6
m6 格式
查找了一下簇宽，列名分別為：
-
qseqidquery (e.g., unknown gene) sequence id -
sseqidsubject (e.g., reference genome) sequence id -
pidentpercentage of identical matches -
lengthalignment length (sequence overlap) -
mismatchnumber of mismatches -
gapopennumber of gap openings -
qstartstart of alignment in query -
qendend of alignment in query -
sstartstart of alignment in subject -
sendend of alignment in subject -
evalueexpect value -
bitscorebit score
二、BLAT
https://blog.csdn.net/alnx37271/article/details/101358955?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2aggregatepagefirst_rank_ecpm_v1~rank_v31_ecpm-2-101358955.pc_agg_new_rank&utm_term=blastall+%E4%BD%BF%E7%94%A8&spm=1000.2123.3001.4430
1. 簡(jiǎn)介
Blat吧享，全稱 The BLAST-Like Alignment Tool魏割，可以稱為"類BLAST 比對(duì)工具"。
- 對(duì)于DNA序列钢颂，BLAT是用來設(shè)計(jì)尋找95%及以上相似至少40個(gè)堿基的序列钞它。
- 對(duì)于蛋白序列，BLAT是用來設(shè)計(jì)尋找80%及以上相似至少20個(gè)氨基酸的序列殊鞭。
特點(diǎn)：
- 速度快（直接把數(shù)據(jù)庫索引讀入內(nèi)存遭垛，無需訪問硬盤）；
- 共線性輸出結(jié)果簡(jiǎn)單易讀操灿；
- 對(duì)于比較小的序列和大基因組的比對(duì)锯仪，BLAT是首選；
對(duì)于比較小的序列（如 cDNA 等）對(duì)大基因組的blat與blast比較比對(duì)趾盐，blat 無疑是首選庶喜。Blat 把相關(guān)的呈共線性的比對(duì)結(jié)果連接成為更大的 比對(duì)結(jié)果小腊，從中也可以很容易的找到 exons 和 introns。因此久窟，在相近物種的基因同源性分析和EST 分析中秩冈，blat 得到了廣泛的應(yīng)用。
Blat 的比對(duì)速度之所以能比Blast快幾百倍,是因?yàn)榇藘烧咧g的比對(duì)機(jī)制有著本質(zhì)的差別斥扛。
- Blast是將查詢序列索引化漩仙，然后線性搜索龐大的目標(biāo)數(shù)據(jù)庫，期間頻繁地訪問硬盤數(shù)據(jù)犹赖，時(shí)間和空間上的數(shù)據(jù)相關(guān)性較卸铀；
- Blast 則將龐大的目標(biāo)數(shù)據(jù)庫索引化峻村，然后線性搜索查詢序列麸折，這種搜索方式在時(shí)間和空間上的數(shù)據(jù)相關(guān)性比較大。
Blat將數(shù)據(jù)庫索引一次性讀入內(nèi)存粘昨，可以反復(fù)地高速調(diào)用垢啼，無需訪問硬盤，占用的系統(tǒng)資源很少张肾。只要索引建立芭析，查詢序列的量越大，Blat的優(yōu)勢(shì)就越明顯吞瞪。
2. 原理
首先 blat 將參考序列拆分成tiles/kmers馁启，其拆分的方式取決于兩個(gè)參數(shù)：
-
-tileSize決定tiles/kmers的大小，一般設(shè)定范圍是：8-12芍秆，預(yù)設(shè)DNA為11惯疙，蛋白質(zhì)為5； -
-stepSize決定tiles/kmers移動(dòng)的步長妖啥；
3. 軟件下載與安裝
簡(jiǎn)單方便霉颠，只需無腦輸入以下命令：
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat
chmod u+x blat
4. 軟件運(yùn)行（比對(duì)）
命令如下：
用法：
blat database query [-ooc=11.ooc] output.psl
blat nt test.fa test.out -out=blast8
說明：
- blat有很多參數(shù)可以選擇，但大部分時(shí)候我們按照默認(rèn)的即可荆虱。
- blat的輸入文件和待查詢數(shù)據(jù)庫的格式蒿偎，可以是fa/nib/2bit三種格式中的任意一種。運(yùn)行時(shí)十分簡(jiǎn)單怀读，不需要進(jìn)行建庫就可以直接比對(duì)诉位。
- 軟件運(yùn)行時(shí)，依次輸入軟件名（軟件執(zhí)行路徑）愿吹，待比對(duì)的數(shù)據(jù)庫不从，待搜索的序列，輸出結(jié)果犁跪。順序前后不能顛倒椿息。這樣就可以開始比對(duì)了
- 程序開始運(yùn)行后歹袁，會(huì)在讀完database 中的所有subject 序列時(shí)在屏幕輸出database的統(tǒng)計(jì)結(jié)果。
以下是一些常用的參數(shù)：
-
-noHead: 不輸出表頭信息寝优，有助于結(jié)果文件的后續(xù)處理 -
-out： 指定輸出的文件格式条舔，此處使用的是blast的m8格式 -
-t： 指定數(shù)據(jù)庫的類型，dna/prot/dnax -
-q: 指定序列類型乏矾，dna/rna/prot/dnax/rna
blat詳細(xì)參數(shù)
用法：blat database query [-ooc=11.ooc] output.psl
database 輸入文件必須是其中一種類型：a .fa , .nib or .2bit file
query 輸入文件必須是其中一種類型：a .fa , .nib or .2bit file
output.psl 輸出文件
-t=type 數(shù)據(jù)庫類型孟抗，可選項(xiàng): dna/prot/dnax
-q=type 查詢序列的類型，可選項(xiàng)：dna/prot/dnax/rnax
-prot 等同于 -t=prot -q=prot
-ooc=N.ooc Use overused tile file N.ooc. N should correspond to the tileSize
-tileSize=N 設(shè)定tiles/kmers的大小
-stepSize=N 設(shè)定tiles/kmers在比對(duì)時(shí)移動(dòng)的步長钻心，即兩個(gè)相鄰tiles/kmers之間的距離凄硼，預(yù)設(shè)值是tileSize
-oneOff=N 如果設(shè)定為 1 ，則表示在比對(duì)到tile上允許有一個(gè)錯(cuò)配堿基(mismatch)捷沸，預(yù)設(shè)值是0
-minMatch=N 設(shè)定至少匹配的tile的個(gè)數(shù)摊沉，一般設(shè)置值的范圍是2-4，通常核苷酸的預(yù)設(shè)值為2痒给，蛋白質(zhì)的預(yù)設(shè)值為1
-minScore=N 設(shè)定最小分值说墨。 由于indel通常會(huì)對(duì)序列的功能產(chǎn)生影響，所以空位在比對(duì)過程中總是對(duì)應(yīng)于一個(gè)負(fù)分苍柏，也就是所謂的空位罰分(Gap penalty)尼斧。根據(jù)打分機(jī)制，這個(gè)分值等于匹配堿基分值減去替換分值(mismatch)和空位罰分试吁。預(yù)設(shè)值為30
-minIdentity=N 設(shè)置序列相似度(sequence identity)最小百分比棺棵。通常核苷酸(nucleotide searches)預(yù)設(shè)值為90，蛋白質(zhì)和翻譯蛋白(protein or translated protein searches)預(yù)設(shè)值為25
-maxGap=N 在一定長度序列中潘悼，設(shè)定兩個(gè)tiles/kmers之間的允許最大的空位(gap)大小律秃。通常設(shè)定范圍是0-3爬橡，預(yù)設(shè)值為2治唤，且僅在minMatch > 1時(shí)搭配使用
-noHead 抑制.psl頭文件的輸出，內(nèi)容全部均是以制表符為分隔符的文件
-makeOoc=N.ooc Make overused tile file. Target needs to be complete genome.
-repMatch=N 在一段序列被標(biāo)記為overused之前糙申，設(shè)定允許tiles/kmers重復(fù)次數(shù)宾添。如果超過設(shè)定值，該tiles/kmers將會(huì)被標(biāo)記為overused柜裸。通常當(dāng)tileSize設(shè)定為12時(shí)缕陕，repMatch則設(shè)定為256；當(dāng)tileSize設(shè)定為11時(shí)疙挺，repMatch則設(shè)定為1024扛邑；當(dāng)tileSize設(shè)定為10時(shí)，repMatch則設(shè)定為4096铐然。
-mask=type Mask out repeats. Alignments won't be started in masked region but may extend through it in nucleotide searches. Masked areas are ignored entirely in protein or translated searches. Types are
lower - mask out lower cased sequence
upper - mask out upper cased sequence
out - mask according to database.out RepeatMasker .out file
file.out - mask database according to RepeatMasker file.out
-qMask=type Mask out repeats in query sequence. 類型選擇與參數(shù)-mask相同
-repeats=type 類型選擇與參數(shù)-mask相同蔬崩。無論如何重復(fù)堿基不會(huì)被掩蓋(masked),但是在匹配重復(fù)區(qū)域時(shí)將會(huì)在psl輸出文件中會(huì)單獨(dú)展示其匹配結(jié)果恶座，即與其他區(qū)域的匹配結(jié)果是分開的。
-minRepDivergence=NN - minimum percent divergence of repeats to allow them to be unmasked. Default is 15. Only relevant for masking using RepeatMasker .out files.
-dots=N 每N個(gè)序列就輸出一個(gè)點(diǎn)沥阳，用于展示程序運(yùn)行的進(jìn)度
-trimT 剪切首部的poly-T
-noTrimA 不剪切尾部的poly-A
-trimHardA 從psl輸出文件中的qSize和alignments中移除poly-A尾巴
-fastMap 快速的DNA/DNA remapping跨琳，要求查詢序列長度不超過5000、高相似度和不進(jìn)行內(nèi)含子的比對(duì)
-out=type 輸出文件格式桐罕，格式如下：
psl - Default. Tab separated format, no sequence
pslx - Tab separated format with sequence
axt - blastz-associated axt format
maf - multiz-associated maf format
sim4 - similar to sim4 format
wublast - similar to wublast format
blast - similar to NCBI blast format
blast8- NCBI blast tabular format
blast9 - NCBI blast tabular format with comments
-fine 對(duì)于高質(zhì)量的mRNAs搜索small initial和terminal exons更為嚴(yán)苛脉让。此選項(xiàng)不推薦應(yīng)用于ESTs
For high quality mRNAs look harder for small initial and terminal exons.
-maxIntron=N 設(shè)定內(nèi)含子最大的序列長度. Default is 750000
-extendThroughN 允許序列的比對(duì)可以從大段N區(qū)域延伸
使用案例
# 處理單個(gè)job
blat chr11.fa human/test.fa test.psl # 輸出不含序列
blat chr11.fa human/test.fa -out=pslx test.pslx # 輸出含序列
blat chr11.fa human/test.fa -out=blast test.blast # 輸出格式同NABI的blast格式
# 并行處理多個(gè)jobs
time parallel blat chr{}.fa human/human.fa test_{}.psl ::: {1..22} X Y M
5. 問題
1. 內(nèi)存溢出
值得關(guān)注的是，因?yàn)閎lat的算法原理功炮，它是將整個(gè)帶查詢數(shù)據(jù)庫全部讀入內(nèi)存進(jìn)行比對(duì)溅潜。因此如果你的服務(wù)器內(nèi)存不大，不建議使用blat進(jìn)行nt/nr庫的比對(duì)薪伏。
下面給出一個(gè)簡(jiǎn)單的計(jì)算方法伟恶，用于評(píng)估自己的服務(wù)器是否可以順溜的run blat。對(duì)于帶查詢基因組毅该，平均每個(gè)堿基博秫，需要2bytes的內(nèi)存。wiki給出的人類基因組大小為3100Mb眶掌，因此我們大概需要6.2G的內(nèi)存才可以順利的對(duì)人類基因組進(jìn)行blat查詢挡育。
年少無知的我，用128G內(nèi)存的服務(wù)器跑nt數(shù)據(jù)庫朴爬，理所當(dāng)然的把我們課題組的服務(wù)器跑崩了~
https://www.biostars.org/p/9479310/#9479314
參考鏈接：Blat-The BLAST-Like Alignment Tool (詳細(xì)的使用教程)
6. GNU Parallel
安裝編譯
wget ftp://ftp.gnu.org/gnu/parallel/parallel-20170822.tar.bz2
tar -jxvf parallel-20170822.tar.bz2
cd parallel-20170822/
cat README
./configure && make && sudo make install
使用
- parallel教程: http://www.gnu.org/software/parallel/parallel_tutorial.html
- parallel中文版教程: http://my.oschina.net/enyo/blog/271612
- parallel與其他Linux命令的搭配使用: http://www.vaikan.com/use-multiple-cpu-cores-with-your-linux-commands/
7. 網(wǎng)絡(luò)版
鏈接 ：http://genome.ucsc.edu/cgi-bin/hgBlat
操作方法
- Genome：選擇物種即寒，比如人
- Assembly：版本號(hào)
- Query type：用于查詢的序列類型（DNA/蛋白質(zhì)）
- Sort output：結(jié)果排序方式
- Output type：輸出格式
- hyperlink：指向結(jié)果的超鏈接，便于可視化
- psl：制表符分隔的表格召噩，便于數(shù)據(jù)處理
查詢結(jié)果（hyperlink）
| ACTIONS | QUERY | SCORE | START | END | QSIZE | IDENTITY | CHROM | STRAND | START | END | SPAN | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| browser | details | CRP_HUMAN | 671 | 1 | 224 | 224 | 100.0% | chr1 | +- | 159713528 | 159714485 | 958 |
| browser | details | CRP_HUMAN | 105 | 119 | 183 | 224 | 77.0% | chr1 | +- | 159705131 | 159705325 | 195 |
| browser | details | CRP_HUMAN | 54 | 117 | 188 | 224 | 62.5% | chr1 | ++ | 159276797 | 159277012 | 216 |
詳情
點(diǎn)擊Browser可以進(jìn)入詳情界面
查詢結(jié)果（psl）
| match | mismatch | rep. match | N's | Q gap count | Q gap bases | T gap count | T gap bases | strand | Q name | Q size | Q start | Q end | T name | T size | T start | T end | block count | blockSizes | qStarts | tStarts |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 224 | 0 | 0 | 0 | 0 | 0 | 1 | 286 | +- | CRP_HUMAN | 224 | 0 | 224 | chr1 | 248956422 | 159713527 | 159714485 | 2 | 19,205, | 0,19, | 89241937,89242280, |
| 50 | 15 | 0 | 0 | 0 | 0 | 0 | 0 | +- | CRP_HUMAN | 224 | 118 | 183 | chr1 | 248956422 | 159705130 | 159705325 | 1 | 65, | 118, | 89251097, |
| 45 | 27 | 0 | 0 | 0 | 0 | 0 | 0 | ++ | CRP_HUMAN | 224 | 116 | 188 | chr1 | 248956422 | 159276796 | 159277012 | 1 | 72, | 116, | 159276796, |
三母赵、diamond
diamond作為一個(gè)和BLAST具有相似功能的軟件，具有以下特點(diǎn)：
- 比BLAST快500到20,000倍
- 長序列的移框聯(lián)配分析(frameshift alignment)
- 資源消耗小具滴，普通臺(tái)式機(jī)和筆記本都能運(yùn)行
- 輸出格式多樣
軟件的安裝很簡(jiǎn)單凹嘲，以下是具體命令：
wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar -zxzf diamond-linux64.tar.gz
diamond 的功能是將蛋白序列或者其翻譯后的核苷酸和蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對(duì)，與blast相比功能單一构韵，但也讓它的使用格外的簡(jiǎn)單周蹭。
不推薦將核苷酸序列與蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對(duì)，因?yàn)榭赡苡性S多序列是比對(duì)到非編碼序列上的疲恢，利用蛋白質(zhì)數(shù)據(jù)庫進(jìn)行比對(duì)凶朗，將導(dǎo)致假陰性過高。
下載nr數(shù)據(jù)庫并建庫
具體命令如下：
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
gunzip nr.gz
diamond makedb --in nr --db nr
-
--in: 后面跟蛋白質(zhì)數(shù)據(jù)庫 -
--db： 指定生成的diamond數(shù)據(jù)庫名稱
比對(duì)搜索
相當(dāng)簡(jiǎn)單显拳，只有兩個(gè)子命令棚愤，blastx和blastp，前者比對(duì)DNA序列杂数，后者比對(duì)蛋白
diamond blastx --db nr -q reads.fna -o dna_matches_fmt6.txt
diamond blastp --db nr -q reads.faa -o protein_matches_fmt6.txt
參數(shù)簡(jiǎn)單介紹：
-
-q: 輸入的待檢索序列 -
-o：指定輸出文件宛畦，默認(rèn)以 --outfmt 6 格式進(jìn)行輸出 -
--db: 后面跟著我們建立好的diamond 蛋白數(shù)據(jù)庫
參考鏈接
http://www.reibang.com/p/f6f868010949
http://www.reibang.com/p/7d536d8da3fb
