字符串匹配

喜欢喝冰茶

字符串匹配，string match，这个是计算机里面常见的问题，例如：
  k5 ]) y9 l1 _  I% W& A( e
1 e6 c3 B, Z6 E3 t0 |% Ystring1: TACGGCATGGCTATCGTAGCTAG

string2: GCTAT
7 O6 c- F& }) i6 m9 J
要求在string1里找到string2的位置，如果存在多个的话，都要找出来。

/ y  d) r7 U& u3 j  j) l8 D8 ~, b$ G可以自己估计一下时间复杂度，真实的例子是，String1长达3billion，或者6个billion。string2长约一、二百，但是数目可以是以billion计的。

先扛着。

喜欢喝冰茶

最野蛮也是最简单的办法：一个一个比。

5 ]2 k2 N" I" z/ w" p, q" n: w

string1: TACGGCATGGCTATCGTAGCTAG

string2: GCTAT

% x7 l/ p  N7 O# R要求在string1里找到string2的位置，如果存在多个的话，都要找出来。

拿string2和string1比，至少需要比string1的长度减去string2的长度再加1次。在实际应用中，如果string1的长度是10^9，而string2只有一二百，那么做一次基本上就是比10^9次。当然如果很幸运，string2在string1开始的地方，那一次就够了。所以平均来说，要比10^9／2次，也就是O(10^9)。
9 ~* m6 Y: }' z+ `( C5 V- a, n
+ O; C! L& w- v0 r$ Q但是如果实际情况中，有10^6到10^9个string2s，那总共要比多少次？10^15到10^18次。这什么概念？不考虑所有的overhead，比一次只需一个时钟，那3G的CPU，意味着一秒可以比10^9次，要完成这样一个工作，需要10^6到10^9秒，1年＝365天 x 24小时 x 3600秒＝31Millon秒。也就是说，最短大约需要12天，最长需要30年。如果这样的操作做十次，一台CPU要算至少120天到300年！！！人都死几次还没比完，太郁闷了，所以不可接受。
8 o3 y! J; y5 r* \3 k# Z. F" X
那如果是这个样子
" u- `6 n" j. c7 ?8 j, e8 G
string1: AAAAAATTTTCCCCCGGGTTTTAAAACCCCCCGG
string2: TTAAA

是不是会快很多？

继续扛。

喜欢喝冰茶

1 _* y: f' F1 p/ y3 x0 R1 }

如果两个字符串是这个样子

3 _- W/ n" y- Lstring1: AAAAAATTTTCCCCCGGGTTTTAAAACCCCCCGG
string2: TTAAA

当然要省很多时间，因为不需要对string1一个一个比了！！！

" U1 p9 d( N8 S; N5 ustring1可以写成：
长度 字符 起始位置
6      A     1
2 Z4 b! o* l3 E( a4      T     7
+ ]% I3 e2 C; a: t( U5      C     11
7 k# h! D, [  H, Z/ s6 R3      G     16
7 E+ w0 B  S, O9 r: T4      T     19
......

所以当用string2去比的时候，一开始根本就不用考虑字符为A，G，C的行，因为string2开始是T。因此在这个例子中，不需要去检查string1的每个位置，而是非常有限的几个位置，所以可以省很多时间。
* L. F" f1 _" u( r
- m5 q4 N# ~6 t1 z2 _% m$ ~- q那么如果存在一种这样的转换方法能够将主贴里的字符串转化成这种，势必会省很多时间。有这样一种方法吗？哪里去找？
0 i8 u2 z( E& J( r7 j7 i
如果你是有心人，你觉得这个东西最常用在哪里？
! ]: s! I! P7 f) O8 B2 U
对了，文件压缩。

# K+ a% p& T6 o+ n6 R6 L2 Y" V事实上，真正的解决方法就是借鉴了最早用于文件压缩的一种算法，称为Burrows－Wheeler Transform，又称block-sorting compression。这是当年在DEC工作的Michael Burrows和David Wheeler发明的，所以以他们的名字命名，bzip2的压缩文件就是基于该算法的。它的转换其实很简单，如果感兴趣大家可以google／wiki（wikipedia上很详细的操作细节）上去看细节，但简单的来说，就是把一个字符串头围相接，不停的移动一位，然后排序，最后取出最后一列就行了。Burrows－Wheeler Transform的特性就是转换后的字符串相对于原始字符串含有大量的重复字符片段，所以就可以使得我们的问题变的相对快捷。
; {& o* V# z' y; C
$ Q% C, c! o( @9 L5 S2 T那么是否就十全十美，万事大吉了呢？这个需要从实际的具体需求来看。

扛吧，没什么好说的。

喜欢喝冰茶

MacArthur 发表于 2013-10-6 16:17
( F8 @5 c9 ?: t3 i4 |这就是基因BLAST算法的标准定义嘛。。。
4 z& H- J1 k0 D- _2 T+ ~# v
Google ＂Blast算法＂，一堆开源程序。

1 z4 [+ y' I; Q/ t$ E- E) v3 c' Wblast并非只是解决基因的问题，蛋白同样适用，只不过相对于DNA／RNA而言，蛋白质要复杂得多。

喜欢喝冰茶

MacArthur 发表于 2013-10-6 21:35
& ~: M6 L- S( X% s) _( X不明觉厉。。。 上Billion字节的操作，想想就头大。。。 这个规模是不是得上并行计算了？
2 @, M7 s( {! y, B! f. u- Q ...

blast好像支持吧，主要不是分布式的问题，而是blast使用的算法对蛋白的分数的定义是很有讲究的，这些分数的定义大概要涉及到另外一个和进化相关的计算领域。

喜欢喝冰茶

已经有同校指出了blast，一个用于对生物大分子的测序程序。这东西正式诞生于1990年，同年人类基因组计划启动。它不仅可以用于DNA/RNA的测序，还可以对蛋白测序。六年前在河里曾经写过一个搭积木的小玩意儿，其中涉及到利用计算方法来预测实验中很难测量的蛋白质空间结构，最有效的计算方法就是同源模型。在这个模型里面，我们允许寻找的字符串,string2可以不用严格的和string1，也就是模板序列，匹配的。换句话说，就是字符串匹配是容错的。具体到主贴里提到的真实问题的挑战，就是模板字符串是人类的基因组（genome），单链长达3billion结构单位，也就是3billion的字符长度，而DNA是双链的，也就是6B的长度。想一想，每个人都是unique的，所以即使都是“健康人”，每个人的基因组都不会完全一样，因此，这种字符串的匹配一定是允许错误的，否则的话，如果大家都一样，即使算法速度再慢都不是问题，因为只做一次，从时间复杂度上，就是0。
$ E3 W: z$ b* J- B8 @) R
7 q. s4 I) w! b+ M! t2 S那么如果考虑容错的匹配，基于bwt方法的算法将面临一个巨大的挑战，因为BWT是严格转换的，可是容错的实际需要，就要考虑转换前和转换后的字符串的错误问题，这会使的问题非常复杂。因此，现在的解决方案就是，当我们需要更多的关注容错的问题时，我们仍然使用传统的，也就是blast所使用的Smith－Waterman算法。具体的细节，感兴趣的可以google／wiki，基本来讲，这是一种被称之为动态编程的计算机算法，可以很好的考虑容错问题，诸如substitution、insertion、deletion或者indel（也就是insertion和deletion的混合体），而这些“错误”，则广泛的存在于病人的基因组中，特别是癌症基因组。当然，现在所使用的smith－waterman基于的方法都是修改了的，主要是计算机算法上的加速，诸如hash的SW算法。

$ m4 P# t$ w8 C4 G7 o但是，bwt基于的方法，运算速度非常快。曾经有朋友几年前做过测试，拿上一代mac book pro，对于5万个长达100字符的输入string2s，string1是人类第一号染色体的DNA序列。blast要跑a couple of minutes，但是第一代bwt based的应用程序，在用手捋了一下头发之后，结果就出来了，当时朋友还以为自己错了。所以相对于blast所使用的算法而言，以bwt为基础的应用，诸如bowtie／bowtie2，bwa（short read version），soap（主要运行于HPC上，由位于深圳的华大基因开发，他们大概是曙光的最大用户），这些闻名遐迩的应用程序（从文章引用次数上，bowtie是最早的，发表于2009年，引用数已经接近3千，bwa大约2千五，soap约为700。引用它们的文章，不仅包括Nature，Science，Cell还有医学领域里的一些著名杂志，像新英格兰医学杂志），使得新型测量技术得以广泛的应用。如果不考虑监管、法律和伦理方面的问题，在可以预见的将来，你的mobile device里存储自己的DNA序列将不再是梦想而是现实，而人们同时希望，DNA序列的检测成为新生儿的常规检查。

喜欢喝冰茶

7 E/ |; z: N. g! E  ?
8 L( E/ Z0 d; W这个帖子里的东西，看起来似乎是一个简单的计算问题，但是却很可能是一场改变人类健康革命的基础。正是由于高速有效的匹配算法的发展，才使得新的测序技术得以广泛应用，不仅是基础研究，而更重要的是临床应用，这里简单的提到一些。等有空的时候，写一些有关基于genome和tanscriptome的技术以及实际应用，像分子诊断技术，这将对癌症的诊断和预期提供巨大的帮助。感兴趣的同学，可以看看AML，也就是急性骨髓性白血病的亚型分类标准，特别是WHO的新标准，还有预期，是否能看懂。
2 m3 Z: O, s. Y; {7 ^5 V; B

喜欢喝冰茶

chalet 发表于 2013-10-10 22:06
9 N1 b  M5 t: U5 B2 m欢迎欢迎。
" x/ q% ?% w. V/ T) I# j
如果说硬要分类，这个应该还是属于生物信息学领域吧，在科教学圃应该很对口啊。

; U! M# j: a8 V1 _- ^. g/ }+ x握握手，看起来也是生物计算的啊，现在在做什么？

- O' \$ ?4 Z1 F% f还没想好怎么写，涉及的范围得控制一下，要不太大了，而且需要用些篇幅介绍一下疾病得复杂性，像癌症，这个得加好多分子遗传方面的科普，否则的话会很难理解分子诊断方面所面临的挑战。
3 `1 `) P% \8 X, t+ i" `
9 @+ u  p: O& v4 ?# {5 O. c+ X$ \昨天刚在八卦里贴了一个，http://www.aswetalk.org/bbs/thread-25733-1-1.html，这是一个很好的分子诊断的例子，并且符合现代医学的发展趋势。以后尽可能的在八卦里发一些小东西，最后写的东西应该都能用的上，慢慢来吧。

喜欢喝冰茶

chalet 发表于 2013-10-11 20:29
我可不是作生物计算的，是学临床医学出身的，而数学正是我的致命伤，哭啊~~~
) I+ d  s1 Y6 ?那是上一轮生物技术泡沫的时 ...

/ C0 N7 k2 l5 q9 h$ [$ \  G7 w兄弟原来是医生，幸会幸会，我的很多合作者都是MD。

& f+ h, a; w) H% ~呵呵，03年太早了，HGP刚完成，那会儿还没看出个所以然来。到07年，新一代测序技术引入市场，一场“革命”就真得来了。对于癌症，人们了解的越来越多，这方面TCGA功不可没。
/ L$ I) [( J" E
! |& T1 k& K7 \9 V事实上，大家谈论的是分子诊断，而非基因诊断。虽然大多数生物学家和医生仍然以基因作为基本单元考虑疾病，但在做生物医学信息的人看来，基因是个太大的概念，也是太不准确的概念，例如，临床常用的mRNA的表达，一般是以基因为基本单位的。但是，基因里面有exon，如果一个exon正向表达，另一个负向表达，从基因层次上可能没区别，但是其实已经涉及到了isoform或者splicing的问题，也就是蛋白的编码有可能被改变了，自然功能就会发生改变，那么这种异常就可以成为该种疾病的一种特征，也就是biomarker。

幸运的是，越来越多的人开始意识到这个问题。我给医学院代的一门课的听众就是MD student，MD fellow和一些即是医生也是faculty的同学们，大家也越来越意识到新的技术正在颠覆我们的概念。

喜欢喝冰茶

chalet 发表于 2013-10-12 03:12
1 F; b, ?7 {8 |# D" G$ X非常赞同你说的。确实当年我去那家公司的时候，他们拿手的是cDNA表达谱芯片，后面的事实证明，这个层面的 ...

% f8 a8 A' S7 t% N# l" r. a
) j9 Q: L  m( M, c0 j
cDNA算是经典的array了，该不会在程氏公司吧~_~。现在mRNA的expression，特别是经典的几十个癌症基因的表达，仍然是对一些癌症进行相对早期诊断的手段之一，只不过受益者比较小众而已。
) m. t8 ^5 I* x! }% |8 U
3 y( F; {, q) ?7 w( Z/ g分子诊断在临床上的大规模应用现在确实没开始，但是这方面的工作很早就铺开了，而且很多临床机构已经参加进来。这个月cell上刚发了一篇TCGA网络的作者们有关glioblastoma的文章，还没看完，刚看了个开头，不过又看到一个朋友的名字在上面。除了像Washington Univ at St louis(它有和Johnhopkins，harvard齐名的医学院)，broad institute这样的学术机构(这个机构其实是很应用的)，如果对美国的癌症治疗了解的话，会看到很多著名的医疗机构，诸如sloan-kettering, dana-farber,md anderson, mayo clinic, fred hutchinson等等这些名闻遐迩的癌症中心的临床人员也参与进来。大量的临床信息和分子测量信息被综合考虑加以分析，事实上，TCGA的array数据是公开可以下载的。TCGA因为是NCI直接支持的，资源非常充足。无论是这次的GBM还是五月份在New England Journal of Medicine发表的AML，动辄就是几百个病人一下子上五六个平台同时测量分析，使用了当前最先进的设备和分析手段，知道了很多以前不了解的信息，这些东西势必对以后的应用提供很大的帮助。另外一方面，不仅是对第一个问题的补充，同时也算是回答了第二个问题。那就是现在数得着的药厂，都已经或者正在建立分子诊断部门，而他们的实力绝对不可小觑，只不过人家闷声发财罢了。就我所了解的情况，很多医生还是很感兴趣的，不仅是美国的，还有中国的，因为已经有越来越多的病例是依赖新技术而得到治疗的。

然而，诚如你所说，癌症的heterogeneous特性，使得这一领域所面临的挑战远远超过了我们的想象。在AML的研究中，recurrent 的定义是5%，可是看看有多少variants是可以被5%的病人所共享的？这东西不像别的学科，了解的越多，会越来越明朗，而是知道的越多，会越来越困惑。没办法，人自己设计的东西和自然选择出来的东西，论精巧和控制真的不能比。但是，现在的努力仍然对以后是有很大帮助的。例如现在被部分人批评的GWAS研究，不可否认那东西烧了不少钱，而结果相对有限。但是23&me能卖100块钱的疾病风险评估的基础，不就是GWAS鉴定出来的一两千个和疾病关联的SNP吗？

坦白的说，不要说personalized medicine了，就是diagnosis在我们这代人的有生之年都没戏，因为疾病太复杂。不要试图解决所有的问题，但是现在的工作仍然有很大的意义。例如，根据NCI和ACA的数据，美国刚刚过去的财年，确诊了大约一万五千个AML病人，但是也死掉了一万多点儿的病人，当然不都是这一万五里的。你也许知道这种癌症是现在常规手段根本无法诊断和发现的，因为癌细胞在骨髓里。所以不要指望能一下子解决它，但是如果现在的工作能够使一百个病人受益，虽然统计上没有什么意义，毕竟不到1%，但是对这一百个家庭，那又意味着什么？

4 W/ O) T& L/ I7 o& S9 u* |+ V2 T
1 S2 r; Z: j- j) s0 g* g
! r% Y5 h( v8 Z. D$ f% E, y

喜欢喝冰茶

chalet 发表于 2013-10-13 21:20
' T" ^1 ]3 Z( d2 P我对当前这个领域的研究有2个观点：
1. 当前这种研究思路，有点花大本钱做剥丝抽茧的小事的味道。一方面 ...

. b" t" @2 e# j呵呵，那你觉得什么样的思路才是花小钱办大事的呢？HGP刚开始的时候，也有人觉得是浪费钱。没有大量的片面，何来所谓的全面。这样做的动机，其实很简单，就是大量的证据表明，很多疾病，特别是癌症，大部分获得性癌症病人的DNA，RNA序列上有非常明显的异常，有些时候，连形态学上都能看出来。那么大家自然就想既然生出来是“正常的”，发病了就不正常了，所以DNA／RNA序列上的变化一定是时间的“函数”。既然现在没什么别的技术系统地可以用于癌症的早期诊断，那分析DNA/RNA异常大概是最可能也是最现实的方法。

# K1 z+ Q' H: o! D& ~从来没有人，再没有确认更好的方法前，去准备取代现在的手段。而且一定要明确得是，基本上，现在得研究给出的是风险，而非确认。毕竟DNA/RNA上观测到的异常只是癌症一个因素，癌症的发展还和“环境”有很大关联的。对现在这种diagnosis方面的工作，相当一部分人存在一种误区。有点像一些实验生物的人看计算生物学一样，开始以为这玩意儿什么都能干，结果发现蛮不是那回事，然后就弃之不用。都不想想，真要是啥都能做，做实验的不早失业了？要真是没用，这门学科也早死了。整天吵吵garbage in garbage out的同学们，你们怎么知道人家input的就是garbage呢？远在这帮做实验的同学们去质疑计算是否正确之前，做biocomputing的早想到了这个问题，并且已经从实验中去寻找间接证据去确认了。把新生的东西当成一种手段就是了，不要排斥，和已有的成熟手段一起使用就行，只要能提高癌症病人的生存率，那每一种方法都有意义。事实上新技术在临床上的成功应用都是和其他手段一起使用的，毕竟这些技术不能治疗，像著名的Nic案例，真正的治疗手段仍然是骨髓移植，但是新技术对最后确定手术起了很大的作用。

' g- R$ m) g+ V
# f' l7 S: ^1 `! R" t* r  \! o

喜欢喝冰茶

一叶飞刀 发表于 2014-11-15 07:01
/ P# q6 d; B0 p4 L/ s关于字符串匹配，应当已经解决完毕了，大概不会有更高级的算法了。
$ F4 }/ \$ x. b5 B
$ {/ e" J) {0 C3 l8 H7 X+ H9 r# ^从S中找ss简单匹配算法为用ss的第一个 ...

: W" k$ S5 j' O6 r9 b8 Q嗯，perfect match不是大问题，问题在于容错匹配，甚至有时候只有头上或者尾巴上的部分匹配。这部分在RNA seq或者DNA seq中的translocation中应用很广。