Ⅰ 基因組序列比對演算法介紹(一)
基因組重測序中序列比對介紹
重測序基因組數據比對,是指將測序儀下機fastq數據(NGS read序列,通常100-150bp),與人類參考基因組(reference)進行匹配,允許錯配(mismatch),插入缺失(indel),目的是在參考基因組找到序列最相似的位置,通常是基因組分析(包括 variation calling,ChIP-seq,RNA-seq,BS-seq)流程的第一步。
常用演算法
圖一
漢明距離(Hamming distance)表示兩個(相同長度)字對應位置不同的數量,我們以d(x,y)表示兩個字x,y之間的漢明距離。對兩個字元串進行異或運算,並統計結果為1的個數,那麼這個數就是漢明距離。圖中read1最佳位置的方法,就是通過查找最小漢明距離的實現的。
編輯距離(Edit distance)是針對二個字元串(例如英文字)的差異程度的量化量測,量測方式是看至少需要多少次的處理才能將一個字元串變成另一個字元串。圖中read3最佳位置,通過查找最我輯距離的方法實現。
圖二
全局比對(Global alignment):全局比對是指將參與比對的兩條序列裡面的所有字元進行比對。全局比對在全局范圍內對兩條序列進行比對打分,找出最佳比對,主要被用來尋找關系密切的序列。其可以用來鑒別或證明新序列與已知序列家族的同源性,是進行分子進化分析的重要前提。其代春梁脊表是Needleman-Wunsch演算法。圖一中,read3使用全部比對。
局部比對(Local alignment):與全局比對不同,局部比對不必對兩個完整的序列進行比對,而是在每個序列中使用某些扒滲局部區域片段進行比對。其產生的需求在於、人們發現有的蛋白序列雖然在序列整體上表現出較大的差異性,但是在某些局部區域能獨立的發揮相同的功能,序列相當保守。這時候依靠全局比對明顯不能得到這些局部相似序列的。其次,在真核生物的基因中,內含子片段表現出了極大變異性,外顯子區域卻較為保守,這時候全局比對表現出了其局限性,無法找出這些局部相似性序列。其代表是Smith-Waterman局部比對演算法。圖一中,read2使用局部比對。
圖三
Smith-Waterman演算法介紹
Smith-Waterman是由Temple F. Smith和Michael S. Waterman於1981年提出的一種進行局部序列比對(相對於全局比對)的演算法,用於找出兩個核苷酸序列或蛋白質序列之間的相似區域。該演算法的目的不是進行全序列的比對,而是找出兩個序列中具有高相似度的片段。S-W演算法基於動態規劃,它接受任意長度、任意位渣戚置、任意序列的對齊,並確定是否能找到最優的比對。
簡單地說就是,動態規劃找到問題中較小部分的解,然後把它們放在一起,形成整個問題的一個完整的最優最終解。
它優於BLAST和FASTA演算法,因為它搜索了更大的可能性,具有更高的敏感性。
S-W演算法不是一次查看整個序列,而是對多個長度的片段進行比較,尋找能夠最大化得分的片段。演算法本身本質上是遞歸的:
圖四
演算法步驟如下:
基因組分析***** 微信 公眾號推出 《50篇文章深入理解NGS》系列文章, 第三篇文章 《基因組序列比對演算法介紹(一)》,爭取每周更新一篇高質量生信干貨帖子。
關注 "基因組分析" 微信公眾號,了解最新最全生信分析知識。
Ⅱ 編輯距離的演算法
比如要計算cafe和coffee的編輯距離。cafe→caffe→coffe→coffee
先創建一個6×8的表(cafe長度為4,coffee長度為6,各加2)
(1): coffeecafe表1接著,在如下位置填入數字(表2): coffee0123456c1a2f3e4表2從3,3格開始,開始計算。取以下三個值的最小值: 如果最上方的字元等於最左方的字元,則為左上方的數字。否則為左上方的數字+1。(對於3,3來說為0) 左方數字+1(對於3,3格來說為2) 上方數字+1(對於3,3格來說為2) 因此為格3,3為0(表3) coffee0123456c10a2f3e4表3循環操作,推出下表 取右下角,得編輯距離為3 動態規劃經常被用來作為這個問題的解決手段之一。
整數 Levenshtein距離(字元串 str1[1..m], 字元串 str2[1..n])
//聲明變數, d[i , j]用於記錄str1[1...i]與str2[1..j]的Levenshtein距離
int d[0..m, 0..n]
//初始化
for i from 0 to m
d[i, 0] := i
for j from 0 to n
d[0, j] := j
//用動態規劃方法計算Levenshtein距離
for i from 1 to m
for j from 1 to n
{
//計算替換操作的代價,如果兩個字元相同,則替換操作代價為0,否則為1
if str1[i]== str2[j] then cost := 0
else cost := 1
//d[i,j]的Levenshtein距離,可以有
d[i, j] := minimum(
d[i-1, j] + 1, //在str1上i位置刪除字元(或者在str2上j-1位置插入字元)
d[i, j-1] + 1, //在str1上i-1位置插入字元(或者在str2上j位置刪除字元)
d[i-1, j-1] + cost // 替換操作
)
}
//返回d[m, n]
return d[m, n]
wikisource上有不同的編程語言的版本。
Ⅲ 字元串距離
又稱Levenshtein距離,是編輯距離(edit distance)的一種。指兩個字串之間,由一個轉成另一個所需的最少編輯操作次數。許可的編輯操作包括將一個字元替換成另一個字元,插入一個字元,刪除一個字元。
一種局部敏感hash,它也是Google公司進行海量網頁去重使用的主要演算法。
傳統的Hash演算法只負責將原始內容盡量均勻隨機地映射為一個簽名值,原理上僅相當於偽隨機數產生演算法。傳統的hash演算法產生的兩個簽名,如果原始內容在一定概率下是相等的;如果不相等,除了說明原始內容不相等外,不再提供任何信息,因為即使原始內容只相差一個位元組,所產生的簽名也很可能差別很大。所以傳統的Hash是無法在簽名的維度上來衡量原內容的相似度,而SimHash本身屬於一種局部敏感哈希演算法,它產生的hash簽名在一定程度上可以表徵原內容的相似度。
我們主要解決的是文本相似度計算,要比較的是兩個文章是否相似,當然我們降維生成了hash簽名也是用於這個目的。看到這里估計大家就明白了,我們使用的simhash就算把文章中的字元串變成 01 串也還是可以用於計算相似度的,而傳統的hash卻不行。
流程
在資訊理論中,兩個等長字元串之間的漢明距離(英語:Hamming distance)是兩個字元串對應位置的不同字元的個數。換句話說,它就是將一個字元串變換成另外一個字元串所需要替換的字元個數。
漢明重量是字元串相對於同樣長度的零字元串的漢明距離,也就是說,它是字元串中非零的元素個數:對於二進制字元串來說,就是1的個數,所以11101的漢明重量是4。
例如:
1011101與1001001之間的漢明距離是2
xlturing/simhashJava
Ⅳ 怎麼計算一個字元串公式的值
編輯距離
關於兩個字元串s1,s2的差別,可以通過計算他們的最小編輯距離來決定。
所謂的編輯距離:
讓s1和s2變成相同字元串需要下面操作的最小次數。
1.
把某個字元ch1變成ch2
2.
刪除某個字元
3.
插入某個字元
例如
s1
=
「12433」
和s2=」1233」;
則可以通過在s2中間插入4得到12433與s1一致。
即
d(s1,s2)
=
1
(進行了一次插入操作)
編輯距離的性質
計算兩個字元串s1+ch1,
s2+ch2的編輯距離有這樣的性質:
1.
d(s1,」」)
=
d(「」,s1)
=
|s1|
d(「ch1」,」ch2」)
=
ch1
==
ch2
?
0
:
1;
2.
d(s1+ch1,s2+ch2)
=
min(
d(s1,s2)+
ch1==ch2
?
0
:
1
,
d(s1+ch1,s2),
d(s1,s2+ch2)
);
第一個性質是顯然的。
第二個性質:
由於我們定義的三個操作來作為編輯距離的一種衡量方法。
於是對ch1,ch2可能的操作只有
1.
把ch1變成ch2
2.
s1+ch1後刪除ch1
d
=
(1+d(s1,s2+ch2))
3.
s1+ch1後插入ch2
d
=
(1
+
d(s1+ch1,s2))
對於2和3的操作可以等價於:
_2.
s2+ch2後添加ch1
d=(1+d(s1,s2+ch2))
_3.
s2+ch2後刪除ch2
d=(1+d(s1+ch1,s2))
因此可以得到計算編輯距離的性質2。
復雜度分析
從上面性質2可以看出計算過程呈現這樣的一種結構(假設各個層用當前計算的串長度標記,並假設兩個串長度都為
n
)
可以看到,該問題的復雜度為指數級別
3
的
n
次方,對於較長的串,時間上是無法讓人忍受的。
分析:
在上面的結構中,我們發現多次出現了
(n-1,n-1),
(n-1,n-2)……。換句話說該結構具有重疊子問題。再加上前面性質2所具有的最優子結構。符合動態規劃演算法基本要素。因此可以使用動態規劃演算法把復雜度降低到多項式級別。
動態規劃求解
首先為了避免重復計運算元問題,添加兩個輔助數組。
一.
保存子問題結果。
m[
|s1|
,|s2|
]
,
其中m[
i
,
j
]
表示子串
s1(0->i)
與
s2(0->j)
的編輯距離
二.
保存字元之間的編輯距離.
e[
|s1|,
|s2|
]
,
其中
e[
i,
j
]
=
s[i]
=
s[j]
?
0
:
1
三.
新的計算表達式
根據性質1得到
m[
0,0]
=
0;
m[
s1i,
0
]
=
|s1i|;
m[
0,
s2j
]
=
|s2j|;
根據性質2得到
m[
i,
j
]
=
min(
m[i-1,j-1]
+
e[
i,
j
]
,
m[i,
j-1]
,
m[i-1,
j]
);
Ⅳ 如何度量兩個詞之間的語義相似度
如何度量句子的語義相似度,很容易想到的是向量空間模型(VSM)和編輯距離的方法,比如A:「我爸是李剛」,B:「我兒子是李剛」,利用VSM方法A(我,爸,是,李剛)B(我,兒子,是,李剛),計算兩個向量的夾角餘弦值,不贅述;編輯距離就更好說了將「爸」,「兒子」分別替換掉,D(A,B)= replace_cost;
這是兩種相當呆的方法,屬於baseline中的baseline,換兩個例子看一下就知道A:「樓房如何建造?」,B:「高爾夫球怎麼打?」,C:「房子怎麼蓋?」,如果用VSM算很明顯由於B,C中有共同的詞「怎麼」,所以BC相似度高於AC;編輯距離同理;
解決這種問題方法也不難,只要通過同義詞詞典對所有句子進行擴展,「如何」、「怎麼」,「樓房」、「房子」都是同義詞或者近義詞,擴展後再算vsm或者edit distance對這一問題即可正解。這種方法一定程度上解決了召回率低的問題,但是擴展後引入雜訊在所難免,尤其若原句中含有多義詞時。例如:「打醬油」、「打毛衣」。在漢字中有些單字詞表達了相當多的意義,在董振東先生的知網(hownet)中對這種類型漢字有很好的語義關系解釋,通過hownet中詞語到義元的樹狀結構可以對對詞語粒度的形似度進行度量。
問題到這里似乎得到了不錯的解答,但實際中遠遠不夠。VSM的方法把句子中的詞語看做相互獨立的特徵,忽略了句子序列關系、位置關系對句子語義的影響;Edit Distance考慮了句子中詞語順序關系,但是這種關系是機械的置換、移動、刪除、添加,實際中每個詞語表達了不同的信息量,同樣的詞語在不同詞語組合中包含的信息量或者說表達的語義信息大不相同。What about 句法分析,計算句法樹的相似度?這個比前兩種方法更靠譜些,因為句法樹很好的描述了詞語在句子中的地位。實際效果要待實驗證實。
對了,還有一種方法translation model,IBM在機器翻譯領域的一大創舉,需要有大量的語料庫進行訓練才能得到理想的翻譯結果。當然包括中間詞語對齊結果,如果能夠利用web資源建立一個高質量的語料庫對兩兩相似句對通過EM迭代詞語對齊,由詞語對齊生成句子相似度,這個。。想想還是不錯的方法!