基本概念

相似性（similarity）

• 一種很直接的數(shù)量關(guān)系，比如部分相同或相似的百分比或其他一些合適的度量
• 如：A序列和B序列的相似性是80%

同源性（homology）

• 從一些數(shù)據(jù)中推斷出的兩個(gè)基因或者蛋白序列具有共同祖先的結(jié)論，屬于質(zhì)的判斷
• 可以說(shuō)A序列和B序列是同源序列，但不能說(shuō)同源性80%

常用工具

• BLAST
• BLAT
  序列比對(duì)的常用工具：BLAST，但是其運(yùn)行速度慢的令人捉急。

一、BLAST（Basic Local Alignment Search Tool，局部相似性基本查詢工具）

BLAST（Basic Local Alignment Search Tool）是一套在蛋白質(zhì)數(shù)據(jù)庫(kù)或DNA數(shù)據(jù)庫(kù)中進(jìn)行相似性比較的分析工具。BLAST程序能迅速與公開數(shù)據(jù)庫(kù)進(jìn)行相似性序列比較。BLAST結(jié)果中的得分是對(duì)一種對(duì)相似性的統(tǒng)計(jì)說(shuō)明。

https://www.biomart.cn/experiment/599/608/19912_0.htm
-F Filter query sequence (DUST with blastn, SEG with others) [String] default = T
查詢序列過(guò)濾，將那些 給出影響比對(duì)結(jié)果的低復(fù)雜度區(qū)域過(guò)濾掉。用blastn進(jìn)行查詢的序列用DUST程序過(guò)濾，其他的用SEG過(guò)濾 。對(duì)DUST和SEG的詳細(xì)情況，用戶可以自己查詢資料。
使用此參數(shù) -F F 即可獲得完整的比對(duì)

在對(duì)核苷酸或蛋白質(zhì)序列數(shù)據(jù)庫(kù)進(jìn)行Blast搜索之前，必須要對(duì)所使用的序列數(shù)據(jù)庫(kù)進(jìn)行formatdb， 即對(duì)序列數(shù) 據(jù)庫(kù)進(jìn)行格式化，這是所有使用BLAST所必須的一步。

1.格式化序列數(shù)據(jù)庫(kù)— —formatdb

formatdb 簡(jiǎn)單介紹：formatdb處理的都是格式為 ASN.1和 FASTA，而且不論是核苷酸序列數(shù)據(jù)庫(kù)，還是蛋白質(zhì)序列數(shù)據(jù)庫(kù)；不論是使用Blastall ，還是Blastpgp，Mega Blast應(yīng)用程序，這一步都是不可少的。

主要參數(shù)的說(shuō)明
-i 輸入需要格式化的源數(shù)據(jù)庫(kù)名稱 Optional
-p 文件類型，是核苷酸序列數(shù)據(jù)庫(kù)，還是蛋白質(zhì)序列數(shù)據(jù)庫(kù) T：protein , F：nucleotide
-n 重命名數(shù)據(jù)庫(kù)文件的名稱 ；
-t 數(shù)據(jù)庫(kù)的標(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ù)
用來(lái)過(guò)濾比對(duì)較差的結(jié)果的，用"-e"參數(shù)指定一個(gè)實(shí)數(shù)，blast 會(huì)過(guò)濾掉期望值大于這個(gè)數(shù)的比對(duì)結(jié)果。這樣不但簡(jiǎn)化了結(jié)果，還縮短了運(yùn)行時(shí)間和結(jié)果占用的空間。

2. -p參數(shù)
-p blastp：蛋白序列與蛋白庫(kù)做比對(duì)。
-p blastx：核酸序列對(duì)蛋白庫(kù)的比對(duì)。
-p blastn：核酸序列對(duì)核酸庫(kù)的比對(duì)。
-p tblastn：蛋白序列對(duì)核酸庫(kù)的比對(duì)。

3.-F 參數(shù)
-F(T/F)參數(shù)是用來(lái)屏蔽簡(jiǎn)單重復(fù)和低復(fù)雜度序列的。如果選“T”，程序在比對(duì)過(guò)程中會(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ū)域的長(zhǎng)度
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 格式

查找了一下，列名分別為：

1. qseqid query (e.g., unknown gene) sequence id
2. sseqid subject (e.g., reference genome) sequence id
3. pident percentage of identical matches
4. length alignment length (sequence overlap)
5. mismatch number of mismatches
6. gapopen number of gap openings
7. qstart start of alignment in query
8. qend end of alignment in query
9. sstart start of alignment in subject
10. send end of alignment in subject
11. evalue expect value
12. bitscore bit score

二、BLAT

https://blog.csdn.net/alnx37271/article/details/101358955?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2_{aggregatepage}first_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是用來(lái)設(shè)計(jì)尋找95%及以上相似至少40個(gè)堿基的序列。
• 對(duì)于蛋白序列，BLAT是用來(lái)設(shè)計(jì)尋找80%及以上相似至少20個(gè)氨基酸的序列。

特點(diǎn)：

• 速度快（直接把數(shù)據(jù)庫(kù)索引讀入內(nèi)存，無(wú)需訪問硬盤）；
• 共線性輸出結(jié)果簡(jiǎn)單易讀；
• 對(duì)于比較小的序列和大基因組的比對(duì)，BLAT是首選；

對(duì)于比較小的序列（如 cDNA 等）對(duì)大基因組的blat與blast比較比對(duì)，blat 無(wú)疑是首選。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ù)庫(kù)，期間頻繁地訪問硬盤數(shù)據(jù)，時(shí)間和空間上的數(shù)據(jù)相關(guān)性較??；
• Blast 則將龐大的目標(biāo)數(shù)據(jù)庫(kù)索引化，然后線性搜索查詢序列，這種搜索方式在時(shí)間和空間上的數(shù)據(jù)相關(guān)性比較大。

Blat將數(shù)據(jù)庫(kù)索引一次性讀入內(nèi)存，可以反復(fù)地高速調(diào)用，無(wú)需訪問硬盤，占用的系統(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)的步長(zhǎng)；

image

3. 軟件下載與安裝

簡(jiǎn)單方便，只需無(wú)腦輸入以下命令：

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

說(shuō)明：

• blat有很多參數(shù)可以選擇，但大部分時(shí)候我們按照默認(rèn)的即可。
• blat的輸入文件和待查詢數(shù)據(jù)庫(kù)的格式，可以是fa/nib/2bit三種格式中的任意一種。運(yùn)行時(shí)十分簡(jiǎn)單，不需要進(jìn)行建庫(kù)就可以直接比對(duì)。
• 軟件運(yùn)行時(shí)，依次輸入軟件名（軟件執(zhí)行路徑），待比對(duì)的數(shù)據(jù)庫(kù)，待搜索的序列，輸出結(jié)果。順序前后不能顛倒。這樣就可以開始比對(duì)了
• 程序開始運(yùn)行后，會(huì)在讀完database 中的所有subject 序列時(shí)在屏幕輸出database的統(tǒng)計(jì)結(jié)果。

以下是一些常用的參數(shù)：

• -noHead: 不輸出表頭信息，有助于結(jié)果文件的后續(xù)處理
• -out： 指定輸出的文件格式，此處使用的是blast的m8格式
• -t： 指定數(shù)據(jù)庫(kù)的類型，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ù)庫(kù)類型，可選項(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)的步長(zhǎ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ì)過(guò)程中總是對(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     在一定長(zhǎng)度序列中，設(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ù)。如果超過(guò)設(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相同。無(wú)論如何重復(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，要求查詢序列長(zhǎng)度不超過(guò)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)含子最大的序列長(zhǎng)度. 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ù)庫(kù)全部讀入內(nèi)存進(jìn)行比對(duì)。因此如果你的服務(wù)器內(nèi)存不大，不建議使用blat進(jìn)行nt/nr庫(kù)的比對(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查詢。

年少無(wú)知的我，用128G內(nèi)存的服務(wù)器跑nt數(shù)據(jù)庫(kù)，理所當(dāng)然的把我們課題組的服務(wù)器跑崩了~

https://www.biostars.org/p/9479310/#9479314

參考鏈接：Using blat

參考鏈接：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）

詳情

點(diǎn)擊Browser可以進(jìn)入詳情界面

image

查詢結(jié)果（psl）

三、diamond

diamond作為一個(gè)和BLAST具有相似功能的軟件，具有以下特點(diǎn)：

• 比BLAST快500到20,000倍
• 長(zhǎng)序列的移框聯(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ù)庫(kù)進(jìn)行比對(duì)，與blast相比功能單一，但也讓它的使用格外的簡(jiǎn)單。

不推薦將核苷酸序列與蛋白質(zhì)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，因?yàn)榭赡苡性S多序列是比對(duì)到非編碼序列上的，利用蛋白質(zhì)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，將導(dǎo)致假陰性過(guò)高。

下載nr數(shù)據(jù)庫(kù)并建庫(kù)

具體命令如下：

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ù)庫(kù)
• --db： 指定生成的diamond數(shù)據(jù)庫(kù)名稱

比對(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ù)庫(kù)

參考鏈接
http://www.itdecent.cn/p/f6f868010949
http://www.itdecent.cn/p/7d536d8da3fb