Ⅰ 基因组序列比对算法介绍(一)
基因组重测序中序列比对介绍
重测序基因组数据比对,是指将测序仪下机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迭代词语对齐,由词语对齐生成句子相似度,这个。。想想还是不错的方法!