文本比对算法
① 比对算法总结(二)——基于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();
}
}
}