biostar handbook(六)| 序列聯(lián)配

序列聯(lián)配

序列聯(lián)配是生物信息學(xué)最基礎(chǔ)的概念，因?yàn)榇蠖鄶?shù)數(shù)據(jù)分析分析策略都需要使用聯(lián)配得到的信息。

舉個(gè)簡(jiǎn)單的例子，假設(shè)你手頭上有一些片段'THIS','LI','NE','ISALIGNED', 已知他們來(lái)自于一個(gè)詞，那么原來(lái)這個(gè)詞應(yīng)該是什么樣子。

--ISALIGNED
THIS-------
-----LI----
--------NE-

組成原來(lái)單詞的過(guò)程就是序列聯(lián)配的所做的事情，就是讓碎片化的信息根據(jù)相似的部分進(jìn)行排列，從而推測(cè)出原來(lái)的數(shù)據(jù)。這個(gè)工作看起來(lái)簡(jiǎn)單，但是涉及到上萬(wàn)條碎片的時(shí)候，人腦估計(jì)就崩潰了，而且還容易出錯(cuò)，所有就開(kāi)發(fā)相應(yīng)算法進(jìn)行實(shí)現(xiàn)。

聯(lián)配的表示方法

序列聯(lián)配結(jié)果的展示結(jié)果有兩種，一種是比較可視化適合人類閱讀，一種適合機(jī)器處理。默認(rèn)BLAST得到的結(jié)果就是第一種展示方式，如下

ATGCAAATGACAAATAC
||||   |||.||.|
ATGC---TGATAACT--
  • '-': 表示gap，可以認(rèn)為是候選的InDel.
  • '|': 表示匹配，一摸一樣的那種
  • '.': 表示錯(cuò)配，可能有堿基發(fā)生了突變。

一般而言，在閱讀BLAST結(jié)果信息時(shí)，都是把上面的序列當(dāng)作參考序列，下面的序列是自己提供的序列。既然序列比對(duì)都是相對(duì)而言，那么其實(shí)存放形式還可以以更加緊縮的形式存放，這就是CIGAR(compact idiosyncratic gapped alignment report,緊湊的異質(zhì)間隙排列報(bào)告)，處于SAM格式的第六列。

所以第二列的信息，用 CIGAR 可以表示為: 4M3D3M1X2M1X1M2D ，解釋一下：

  • 4 M(matches)
  • 3 D(deletions),
  • 3 M(matches),
  • 1 X(mismatch),
  • 2 M(matches),
  • 1 X(mismatch),
  • 1 M(match),
  • 2 D(deletions)

根據(jù)這個(gè)CIGAR值和參考序列就能重構(gòu)出我們自己提供的序列

如何計(jì)算聯(lián)配得分

假設(shè)你的序列是ATGAA，和4個(gè)目標(biāo)序列進(jìn)行了比較，那么哪個(gè)目標(biāo)序列和自己的序列最相近呢？

  1      2    3     4
ATGAA ATGAA ATGAA AT-GAA
|.|.| |||.| |||.| || |||
ACGCA ATGCA ATGTA ATCGAA

為了挑選出最好的"聯(lián)配"，就需要按照一定規(guī)則對(duì)聯(lián)配結(jié)果進(jìn)行打分，比如說(shuō)匹配得5分，錯(cuò)配扣4分，出現(xiàn)一個(gè)gap后扣10分，并且在已有g(shù)ap上延伸會(huì)多口0.5分。按照這個(gè)規(guī)則，計(jì)算上面聯(lián)配的得分分別是"7,16,16,15".因此你或許認(rèn)為第二個(gè)和第三個(gè)可能才是你的目標(biāo)序列，畢竟看起來(lái)也是這樣。但是這只是生信分析得到的一個(gè)猜測(cè)而已。其實(shí)沒(méi)有最好的聯(lián)配結(jié)果，也沒(méi)有最好聯(lián)配相對(duì)SCORING選擇。改變計(jì)分規(guī)則，最好的聯(lián)配結(jié)果也會(huì)改變。

但是根據(jù)一些基本信息，還是能夠推理出最優(yōu)選擇。因此在聯(lián)配的時(shí)候，請(qǐng)注意挑選合適的得分矩陣：

另外你還可以去ftp://ftp.ncbi.nlm.nih.gov/blast/matrices查看目前的得分矩陣類型.

拓展閱讀：

序列聯(lián)配的算法

在后續(xù)的操作需要預(yù)先安裝一些軟件和工具

conda install -c emboss
mkdir -p ~/bin
curl http://data.biostarhandbook.com/align/global-align.sh > ~/bin/global-align.sh
curl http://data.biostarhandbook.com/align/local-align.sh > ~/bin/local-align.sh
chmod +x ~/bin/*-align.sh

下面的例子中使用兩條神奇的蛋白序列THISLINE和ISALIGNED。之所以說(shuō)是神奇，是因?yàn)檫@完全是認(rèn)為按照語(yǔ)言習(xí)慣對(duì)蛋白縮寫字母進(jìn)行排序，當(dāng)然這個(gè)例子最早來(lái)自于Marketa Zvelebil和Jeremy Baum的《理解生物信息學(xué)》。

這兩個(gè)shell腳本實(shí)際上是封裝了emboss的needle方法，修改了標(biāo)準(zhǔn)輸入和標(biāo)準(zhǔn)輸出。

全局聯(lián)配：Global alignments

所謂全局聯(lián)配，就是盡可能保證兩條序列的每個(gè)堿基都能配對(duì)，因此會(huì)有比較多的錯(cuò)配和gap，而且不會(huì)對(duì)序列兩端的gap進(jìn)行懲罰。

# 默認(rèn)
./global-align.sh THISLINE ISALIGNED
# 自定義gap打開(kāi)的罰分
./global-align.sh THISLINE ISALIGNED -gapopen 7
image

局部聯(lián)配：Local alignments

局部聯(lián)配不是盡可能把兩條序列進(jìn)行配對(duì)，而是盡可能的找到那些子區(qū)域是能最優(yōu)聯(lián)配。并且按照得分矩陣，產(chǎn)生分?jǐn)?shù)在閾值之內(nèi)的比對(duì)結(jié)果

./local-align.sh THISLINE ISALIGNED
a                  7 NE      8
                     ||
b                  7 NE      8

這個(gè)比對(duì)是依賴于EBLOSUM62得分矩陣，默認(rèn)得分低于11都會(huì)被淘汰掉。我們可以下載其他得分矩陣進(jìn)行比較。

wget -4 ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/NUC.4.4 -q &
wget -4 ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM30 -q &
wget -4 ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM62 -q &
wget -4 ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM90 -q &

比較之后發(fā)現(xiàn)，在BLOSUM90下得到的序列長(zhǎng)度最長(zhǎng)

./local-align.sh THISLINE ISALIGNED -data BLOSUM90
4 SLI-NE
 :|| ||
3 ALIGNE

局部聯(lián)配適用于尋找兩個(gè)序列最相似的區(qū)域。

如果兩個(gè)序列完全相似，那么局部聯(lián)配的結(jié)果和全局聯(lián)配的結(jié)果將會(huì)一致。畢竟每一個(gè)地方都是最好的，總體而言也是最好的。如果只是為了實(shí)現(xiàn)總體最佳，而不顧及每一個(gè)堿基的感受，有一些堿基就得為了實(shí)現(xiàn)大局而舍棄自己了。

在全局聯(lián)配和局部聯(lián)配之間還有一種混合算法，叫做半全局聯(lián)配(semi-global alignments), 僅僅知道這么一回事就行。

聯(lián)配的窘境

由于沒(méi)有上帝視角，我們只能以手頭的信息去推測(cè)未知的部分。序列聯(lián)配的時(shí)候得到的最優(yōu)結(jié)果僅僅是根據(jù)經(jīng)驗(yàn)設(shè)計(jì)的算法運(yùn)行后的結(jié)果。給你兩條序列“CCAAACCCCCCCTCCCCCGCTTC“和”CCAAACCCCCCCCTCCCCCCGCTTC”，

到底是./global-align.sh CCAAACCCCCCCTCCCCCGCTTC CCAAACCCCCCCCTCCCCCCGCTTC得到的結(jié)果接近真相

1 CCAAACCCCCCC--TCCCCCGCTTC     23
  ||||||||||||  .||||||||||
1 CCAAACCCCCCCCTCCCCCCGCTTC     25

還是./global-align.sh CCAAACCCCCCCTCCCCCGCTTC CCAAACCCCCCCCTCCCCCCGCTTC -gapopen 9得到的結(jié)果接近真相呢？

1 CCAAA-CCCCCCCT-CCCCCGCTTC     23
  ||||| |||||||| ||||||||||
1 CCAAACCCCCCCCTCCCCCCGCTTC     25

比對(duì)未必能發(fā)現(xiàn)真實(shí)的插入和缺失，因?yàn)樗惴▋A向于用最簡(jiǎn)單的解釋去協(xié)調(diào)不同的序列，以一句話作為總結(jié)吧：

Alignment reliability depends on the information content of the aligned sequence itself. Alignments that include low complexity regions are typically less reliable. Additional analysis is typically needed to confirm these results.

多序列聯(lián)配

如果要不同物種同一個(gè)基因做進(jìn)化分析，那么就需要先進(jìn)行多序列聯(lián)配，然后才能用特定的軟件確定不同序列之前的進(jìn)化關(guān)系。常用的軟件為: mafft, muscle, clusta-omega, t-coffee等。算法核心思想大多數(shù)基于序列遞增，但是具體實(shí)現(xiàn)有所不同。

以植物轉(zhuǎn)錄因子STAT家族的同源序列作為例子，數(shù)據(jù)來(lái)自于北京大學(xué)開(kāi)發(fā)的planttfdb, 自行去http://planttfdb.cbi.pku.edu.cn/family.php?fam=STAT下載序列。我將序列存放在STAT文件夾下。

一共有214條序列

# 序列所在文件夾
$ seqkit stat STAT/seq.fas
file     format  type     num_seqs  sum_len  min_len  avg_len  max_len
seq.fas  FASTA   Protein       214  152,000      330    710.3    2,301

僅僅挑選出前10條序列

seqkit seq -n STAT/seq.fas | cut -f 1 -d ' '  | head -n 10 > STAT/ids.txt
seqkit grep --pattern-file STAT/ids.txt STAT/seq.fas > STAT/small.fa

使用mafft進(jìn)行多序列聯(lián)配

mafft --clustalout STAT/small.fa > STAT/alignment.maf

結(jié)果可以用less,head,tail查看，僅選取部分用作展示

image
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容