文本比對演算法
① 比對演算法總結(二)——基於BWT索引結構的比對演算法-Bowite1
這是美國馬里蘭大學計算機研究所、生物信息學和計算生物學中心於2009年發表在《Genome Biology》雜志的一篇經典文章，至此以後依賴於BWT索引的比對演算法成為主流。 Bowite 是一款超快速、內存佔用低的短序列比對軟體，適用於將短reads比對至大型參考基因組。採用Burrows-Wheeler 演算法建立索引的Bowite軟體可以在1 CPU時內，將2000萬條reads 比對至人參考基因組，且內存只佔有1.3Gb。於此同時Bowite 採用了新的quality-aware backtracking（質量回溯）演算法，比對過程允許錯配。
在此之前都是採用對reads (SHRiMP, Maq, RMAP,ZOOM) 或者參考基因組 (SOAP)構建哈希表的演算法進行序列比對，該演算法已在上篇文章中進行了介紹 https://www.jianshu.com/p/f5ccff73b181 。
Bowite 採用了一種完全新的索引構建策略，適用於哺乳動物重測序。根據千人基因組計劃數據，Bowite 在35bp PE 序列上的比對速度要比Maq 軟體快35 倍，比SOAP軟體快300倍。Bowite 採用 Burrows-Wheeler 演算法對 full-text minute-space (FM) 構建索引，人參考基因組佔用的內存為1.3 GB。
為了追求速度，Bowite 針對哺乳動物重測序項目進行了很多合理的折中。例如，如果一條reads有多條最優匹配，Bowite 只會輸出一條最優匹配。當輸出的最優匹配也不是完全匹配時，Bowite並不能保證在所有情況下都能輸出最高質量的匹配。在設定了較高的匹配閾值時，一小部分含有多個錯配的reads可能會比對失敗。在默認參數條件下，Bowite 的靈敏度與SOAP 相當，略低於Maq。可以在命令行手動改變參數，在犧牲更多時間的情況下，增加靈敏度，給出reads所有可能的比對結果。目前Bowite 比對的reads長度范圍為4bp - 1024bp。
Bowite 對參考基因組建立索引的方法是 Burrows-Wheeler transform (BWT) 和 FM index。Bowite 建立的人類基因組索引在硬碟上的大小為2.2GB，在比對時的內存為1.3GB。FM index 常用的精確查找方法為 Ferragina 和 Manzini 演算法。Bowite 沒有完全使用該演算法，因為該演算法不允許錯配，不能比對含有測序錯誤和變異的reads。針對這種情況，Bowite引入了新的擴展演算法：quality-aware backtracking 演算法，允許錯配並支持高質量比對；double indexing 策略，避免過度回溯；Bowite比對策略與Maq軟體相似，允許小部分的高質量reads 含有錯配，並且對所有的錯配位點的質量值設置了上限閾值。
BWT 轉換是字元串的可逆性排列，它最早應用於文本數據的壓縮，依賴BWT建立的索引，可以在較低內存下，實現大型文本的有效搜索。它被在生物信息學中有廣泛的應用，包括重復區域計數、全基因組比對、微陣列探針設計、Smith-Waterman 比對到人參考基因組。Burrows-Wheeler transform (BWT) 的轉換步驟如圖1所示：
1、輪轉排序。如圖1a 所示，（1）將字元$ 添加到文本 T （acaacg）的末尾，但需注意其中字元$ 並未實際添加到文本 T 中，且其在字母表中邏輯順序小於 T 中所有出現過的字元。（2） 然後將當前字元串的第一個字元移到最後一位，形成一個新的字元串，再將新的字元串的第一位移到最後一位形成另一個新的字元串，就這樣不斷循環這個過程，直到字元串循環完畢（即$處於第一位）,這樣就形成了一個基於原字元串的字元矩陣M（這一步原圖1a 進行了省略，見下方小圖）。（3） 然後對矩陣M的各行字元按照字典先後順序排序，獲得排序後的字元矩陣 BWM（T），矩陣的最後一列定義為 BWT（T）。 前期經過一個小復雜的過程獲得了BWT（T）列，那這一列到底有什麼用呢？其實BWT（T）列通過簡單的演算法就可以推算出原始文本T的所有信息。而經過轉換之後的BWT（T）列大量重復字元是靠近的，只儲存該列信息，可以大大提高字元壓縮比例。
2、LF-Mapping。圖1a 轉換矩陣 BWM（T）含有一種 'last first (LF) mapping' 的特性，即最後一列L中出現某字元出現的順序與第一列F某字元出現的次序時一致的。根據Supplementary1 圖中演算法1 STEPLEFT 和 演算法2 UNPERMUTE 就可以推算出BWT（T）到 T 的過程， 圖1 b記錄了整個推算過程。 詳細推算過程可參考這個博客介紹： https://blog.csdn.net/stormlovetao/article/details/7048481 。
3、reads精確匹配。使用BWT演算法的最終目的是要將短reads比對到參考基因組上，確定短reads在參考基因組上的具體位置。轉換後的BWT（T）序列，可以利用Supplementary1 圖中演算法3 EXACTMATCH 實現reads的精確匹配。圖1c 列出了 字元串 aac 比對至acaacg 的過程 。 詳細推算過程可參考這篇介紹： https://zhuanlan.hu.com/p/158901556 。
上述的BWT轉換只能用於精確的匹配，但是測序reads是含有測序錯誤和突變的，精確匹配並不適用。這里應用了 backtracking 搜索的演算法，用於允許錯配快速比對 。含有錯配的reads只是一小部分。測序reads的每個鹼基都含有唯一的測序量值，測序質量值越該位點是測序錯誤的可能越大，只有當一條read 的所有錯配的測序質量值總和小於一定閾值時可以允許錯誤匹配。
圖2顯示了精確匹配和非精確匹配的過程，backtracking 搜索過程類似於 EXACTMATCH ，首先計算連續較長的後綴矩陣。如果矩陣中沒有搜索到相應的reads，則演算法會選擇一個已經匹配的查詢位置，替換一個不同鹼基，再次進行匹配。EXACTMATCH搜索從被替換位置之後開始，這樣就可以比對就可以允許一定的錯配。backtracking 過程發生在堆棧結構的上下文中，當有替換產生時，堆棧的結構會增長；當所有結果都不匹配時，堆棧結構會收縮。
Bowite 軟體的搜索演算法是比較貪婪的，Bowite軟體會報出遇到的第一個有效比對，並不一定是在錯配數目和變異質量上的「最佳比對」。沒有查詢最優比對的原因是尋找「最佳比對」會比現有的模型慢2-3倍。而在重測序項目上，速度是更重要的因素。Bowite 也設置了可以輸出多個比對位置（-k）和所有比對位置（-a）的參數，添加這些參數後，比對速度會顯著變慢。
目前的比對軟體會有過度回溯的情況，在reads的3『端花費大量無用時間去回溯。Bowite利用『double indexing』技術減少了過度回溯的發生。簡單來說就是對正向參考基因組進行BWT轉換，稱為 『Forward index』，同時對反向（注意不是互補配對序列,是反向序列）參考基因組也進行BWT轉換，稱為『Mirror index』。 當只允許一個錯配時，比對根據reads是前半段出現錯配，還是後半段出現錯配會有兩種情況：（1）Phase1 將Forward index 載入入內存，不允許查詢reads右半段出現錯配；（2）Phase2 將Mirror index 載入如內存，不允許查詢序列的反向reads右半段（原查詢序列的左半端） 出現錯配。這樣可以避免過度回溯，提高比比對的靈敏度。 但是，如果比對軟體允許一個reads有多個錯配時，仍然會有過度回溯的現象發生，為了減少過度回溯現象的發生，這里將回溯的上限進行了限定（默認值為：125次）。
Bowite 允許使用者在高質量reads的末端（默認是28bp）設置錯配數目（默認的錯配數目是2）。高質量reads末端的28bp序列被稱為 '種子' 序列。這個『種子』序列又可分為兩等份：14bp的高質量末端稱為 『hi-half』（通常位於5『端）,14bp的低質量末端稱為『lo-half』。 如果種子序列只允許2bp 的錯配，比對會出現4 種情況：（1）種子序列中沒有錯配(case1)；（2）hi-half區域沒有錯配，lo-half區域有一個或兩個錯配（case2）;（3）lo-half區域沒有錯配，hi-half區域有一個或兩個錯配（case3）;（4）lo-half區域有一個錯配，hi-half區域有一個錯配（case4）；
在所有情況下，reads的非種子部分允許任意數目的錯配。如圖3所示，Bowite 演算法會根據上面4 種情況交替變化『Forward index』和『Mirror index』比對策略，主要會有三種比對策略。
Bowite 建立一次參考基因組索引後，後續的比對可反復使用該索引。表1和表2列出了在默認參數條件下，Bowite、SOAP、Maq軟體性能的比較。在reads比對率相近的條件下，Bowite軟體的比對速度速度相對於SOAP、Maq軟體有較大的提升。
1、將reads 比對至人參考基因組上，Bowite相對於SOAP和Maq軟體有較大的優勢。它運行的內存非常小（1.2GB），在相同靈敏度下，速度有了較大的提升。
2、Bowite 軟體建立一次參考基因組索引後，後續的比對可反復使用該索引。
3、Bowite 速度快、內存佔用小、靈敏度高主要是因為使用了BWT演算法構建索引、利用回溯演算法允許錯配、採用Double index策略避免過度回溯。
4、Bowite 軟體目前並不支持插入、缺失比對，這個是今後需要努力的方向。
[1] Langmead B . Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome biology, 2009, 10(3):R25.
[2] BWT 推算過程參考博客 https://blog.csdn.net/stormlovetao/article/details/7048481
[3] FM index 精確查匹配過程參考文章 https://zhuanlan.hu.com/p/158901556
② 文本、語音相似度演算法
前段時間公司項目用到了語音識別,圖像識別,視頻識別等,其實不能說是識別,應該說是相似度對比吧,畢竟相似度對比還上升不了到識別哈,等以後有了更深的理解再來討論修改下!這次就當做一個總結吧!
其實它的原理和視頻圖像相似度演算法類似，將一系列的向量,特徵,權重,進行合並,然後降維降到一維,其實這個演算法也就是採用降維技術,將所有的特徵都用一個唯一標識來表示.然後這個標識是經過這個演算法內部的計算,再利用海明距離計算相似度，視頻和圖片是經過漢明距離計算的
文本我們是採用simhash演算法：
1.我們給文本裡面的詞進行分詞,我們是用ik演算法,這個演算法就是while循環,讀取一行,然後調用ik智能分詞的類,智能去切割裡面的分詞;
2.根據裡面的詞頻,simhash演算法會加一個權重,當然,得詞頻達到多少個的時候才會有有權重,這也是它的缺點,一般文本數據較少的時候,他是不準確的,一般數據量在500+;演算法內部的話會將一系列的向量,特徵,權重,進行合並,然後降維降到一維,其實這個演算法也就是採用降維技術,將所有的特徵都用一個唯一標識來表示.然後這個標識是經過這個演算法內部的計算,然後得到的一個指紋簽名；
3.然後對比兩個文本的相似度就是將兩個指紋簽名進行海明距離計算,如果海明距離<8(根據業務和場景去判斷這個值，8是建議，參考)的話,表示兩個相似,小於3的話.表示兩個文本重復.
simhash演算法我們還可以做語音相似度,它的基本原理就是根據傅里葉變換處理得到聲波的形狀。
語音的坡度如果向上我們就用1表示,向下我們就用0表示,這樣的話,我們也可以用二進制碼去描述一首歌曲.得到一個唯一的指紋簽名,對比兩個音頻的相似度就是將兩個指紋簽名進行海明距離計算<8的話,我們就默認兩個音頻相似.
總結：都是把特徵降到一維，然後採用海明距離計算。計算的值小於多少時，就當做是相似。我這邊講的太淺了，實在領悟有限，時間有限，觸摸不深，等下次有新的領悟再來補充！
③ 文本比較有哪些演算法
用java比較兩個文本文件的不同
RangeDifferencer
public class RangeDifferencer {
private static final RangeDifference[] EMPTY_RESULT= new RangeDifference[0];
/* (non Javadoc)
* Cannot be instantiated!
*/
private RangeDifferencer() {
// nothing to do
}
/**
* Finds the differences between two <code>IRangeComparator</code>s.
* The differences are returned as an array of <code>RangeDifference</code>s.
* If no differences are detected an empty array is returned.
*
* @param left the left range comparator
* @param right the right range comparator
* @return an array of range differences, or an empty array if no differences were found
*/
public static RangeDifference[] findDifferences(IRangeComparator left, IRangeComparator right) {
int rightSize= right.getRangeCount();
int leftSize= left.getRangeCount();
//
// Differences matrix:
// only the last d of each diagonal is stored, i.e., lastDiagonal[k] = row of d
//
int diagLen= 2 * Math.max(rightSize, leftSize); // bound on the size of edit script
int maxDiagonal= diagLen;
int lastDiagonal[]= new int[diagLen + 1]; // the row containing the last d
// on diagonal k (lastDiagonal[k] = row)
int origin= diagLen / 2; // origin of diagonal 0
// script corresponding to d[k]
LinkedRangeDifference script[]= new LinkedRangeDifference[diagLen + 1];
int row, col;
// find common prefix
for (row= 0; row < rightSize && row < leftSize && rangesEqual(right, row, left, row) == true;)
row++;
lastDiagonal[origin]= row;
script[origin]= null;
int lower= (row == rightSize) ? origin + 1 : origin - 1;
int upper= (row == leftSize) ? origin - 1 : origin + 1;
if (lower > upper)
return EMPTY_RESULT;
//System.out.println("findDifferences: " + maxDiagonal + " " + lower + " " + upper);
// for each value of the edit distance
for (int d= 1; d <= maxDiagonal; ++d) { // d is the current edit distance
if (right.skipRangeComparison(d, maxDiagonal, left))
return EMPTY_RESULT; // should be something we already found
// for each relevant diagonal (-d, -d+2 ..., d-2, d)
for (int k= lower; k <= upper; k += 2) { // k is the current diagonal
LinkedRangeDifference edit;
if (k == origin - d || k != origin + d && lastDiagonal[k + 1] >= lastDiagonal[k - 1]) {
//
// move down
//
row= lastDiagonal[k + 1] + 1;
edit= new LinkedRangeDifference(script[k + 1], LinkedRangeDifference.DELETE);
} else {
//
// move right
//
row= lastDiagonal[k - 1];
edit= new LinkedRangeDifference(script[k - 1], LinkedRangeDifference.INSERT);
}
col= row + k - origin;
edit.fRightStart= row;
edit.fLeftStart= col;
//Assert.isTrue(k >= 0 && k <= maxDiagonal);
script[k]= edit;
// slide down the diagonal as far as possible
while (row < rightSize && col < leftSize && rangesEqual(right, row, left, col) == true) {
++row;
++col;
}
//Assert.isTrue(k >= 0 && k <= maxDiagonal); // Unreasonable value for diagonal index
lastDiagonal[k]= row;
if (row == rightSize && col == leftSize) {
//showScript(script[k], right, left);
return createDifferencesRanges(script[k]);
}
if (row == rightSize)
lower= k + 2;
if (col == leftSize)
upper= k - 2;
}
--lower;
++upper;
}
// too many differences
//Assert.isTrue(false);
return null;
}
/**
* Finds the differences among two <code>IRangeComparator</code>s.
* In contrast to <code>findDifferences</code>, the result
* contains <code>RangeDifference</code> elements for non-differing ranges too.
*
* @param left the left range comparator
* @param right the right range comparator
* @return an array of range differences
*/
public static RangeDifference[] findRanges(IRangeComparator left, IRangeComparator right) {
RangeDifference[] in= findDifferences(left, right);
List out= new ArrayList();
RangeDifference rd;
int mstart= 0;
int ystart= 0;
for (int i= 0; i < in.length; i++) {
RangeDifference es= in[i];
rd= new RangeDifference(RangeDifference.NOCHANGE, mstart, es.rightStart() - mstart, ystart, es.leftStart() - ystart);
if (rd.maxLength() != 0)
out.add(rd);
out.add(es);
mstart= es.rightEnd();
ystart= es.leftEnd();
}
rd= new RangeDifference(RangeDifference.NOCHANGE, mstart, right.getRangeCount() - mstart, ystart, left.getRangeCount() - ystart);
if (rd.maxLength() > 0)
out.add(rd);
return (RangeDifference[]) out.toArray(EMPTY_RESULT);
}
//---- private methods
/*
* Creates a Vector of DifferencesRanges out of the LinkedRangeDifference.
* It coalesces adjacent changes.
* In addition, indices are changed such that the ranges are 1) open, i.e,
* the end of the range is not included, and 2) are zero based.
*/
private static RangeDifference[] createDifferencesRanges(LinkedRangeDifference start) {
LinkedRangeDifference ep= reverseDifferences(start);
ArrayList result= new ArrayList();
RangeDifference es= null;
while (ep != null) {
es= new RangeDifference(RangeDifference.CHANGE);
if (ep.isInsert()) {
es.fRightStart= ep.fRightStart + 1;
es.fLeftStart= ep.fLeftStart;
RangeDifference b= ep;
do {
ep= ep.getNext();
es.fLeftLength++;
} while (ep != null && ep.isInsert() && ep.fRightStart == b.fRightStart);
} else {
es.fRightStart= ep.fRightStart;
es.fLeftStart= ep.fLeftStart;
RangeDifference a= ep;
//
// deleted lines
//
do {
a= ep;
ep= ep.getNext();
es.fRightLength++;
} while (ep != null && ep.isDelete() && ep.fRightStart == a.fRightStart + 1);
boolean change= (ep != null && ep.isInsert() && ep.fRightStart == a.fRightStart);
if (change) {
RangeDifference b= ep;
//
// replacement lines
//
do {
ep= ep.getNext();
es.fLeftLength++;
} while (ep != null && ep.isInsert() && ep.fRightStart == b.fRightStart);
} else {
es.fLeftLength= 0;
}
es.fLeftStart++; // meaning of range changes from "insert after", to "replace with"
}
//
// the script commands are 1 based, subtract one to make them zero based
//
es.fRightStart--;
es.fLeftStart--;
result.add(es);
}
return (RangeDifference[]) result.toArray(EMPTY_RESULT);
}
/*
* Tests if two ranges are equal
*/
private static boolean rangesEqual(IRangeComparator a, int ai, IRangeComparator b, int bi) {
return a.rangesEqual(ai, b, bi);
}
/*
* Tests whether <code>right</code> and <code>left</code> changed in the same way
*/
private static boolean rangeSpansEqual(IRangeComparator right, int rightStart, int rightLen, IRangeComparator left, int leftStart, int leftLen) {
if (rightLen == leftLen) {
int i= 0;
for (i= 0; i < rightLen; i++) {
if (!rangesEqual(right, rightStart + i, left, leftStart + i))
break;
}
if (i == rightLen)
return true;
}
return false;
}
/*
* Reverses the range differences
*/
private static LinkedRangeDifference reverseDifferences(LinkedRangeDifference start) {
LinkedRangeDifference ep, behind, ahead;
ahead= start;
ep= null;
while (ahead != null) {
behind= ep;
ep= ahead;
ahead= ahead.getNext();
ep.setNext(behind);
}
return ep;
}
}
下面是一段關於如何使用這些類的簡單的測試代碼
public class RangeDifferencerTest extends TestCase {
InputStream left = null;
InputStream right = null;
/**
* @see junit.framework.TestCase#setUp()
*/
protected void setUp() throws Exception {
String file1 = "d:/temp/1.txt";
String file2 = "d:/temp/2.txt";
left = new FileInputStream(new File(file1));
right = new FileInputStream(new File(file2));
super.setUp();
}
/**
* @see junit.framework.TestCase#tearDown()
*/
protected void tearDown() throws Exception {
left.close();
right.close();
super.tearDown();
}
public static void main(String[] args) {
}
/*
* Test method for 'com.greatroad.smbnm.compare.RangeDifferencer.findDifferences(IRangeComparator, IRangeComparator)'
*/
public void testFindDifferences() {
try {
RangeDifference[] rds = RangeDifferencer.findRanges(new LineComparator(left,"GBK"),new LineComparator(right,"GBK"));
if(rds != null ){
for(int i=0; i<rds.length; i++){
RangeDifference rd = rds[i];
int length = rd.leftLength();
System.out.println(
"kind = "+rd.kind()
+",left["+rd.leftStart()+"-"+rd.leftEnd()
+"],right["+rd.rightStart()+"-"+rd.rightEnd()+"]");
}
}
} catch (Exception e) {
e.printStackTrace();
}
}
}