【发布时间】:2015-10-14 00:58:38
【问题描述】:
我有以下任务: 从 30 个字符长的模式序列开始(它实际上是 DNA 序列,以免称它为 P30)我需要在文本文件中找到所有以 (^agacatacag... ) 开头的行都具有精确的 P30,然后是 30 个中的最后 29 个字符, 28 和最多 10 个字符。我需要做的只是删除模式的第一个字符并继续搜索。 为简单起见,我目前需要精确匹配,但允许 1 个不匹配的时间更长(20-30 个字符长的模式)会更好。
我目前相当慢的解决方案是创建一个每行有一个截断模式的 shell 文件,然后 grep[1] 它。这意味着我正在阅读 20 倍的巨大、几 GB 的文本文件,这可能需要一天以上的时间。
我可以切换到 python,创建一个包含所有必需模式的列表/元组,然后只读取一次文件,每个序列循环 20 次,使用 pypy 加快速度。
- 问题 1:有没有比这种循环更快的正则表达式?
- 问题 2:通过切换到更快的编译语言来加快速度是否有意义? (我正在尝试理解 Dlang)
[1] 因为它是 DNA 序列并且要搜索的输入是 FASTQ 格式,所以我使用的是 fqgrep:https://github.com/indraniel/fqgrep 带 tre 库:https://github.com/laurikari/tre/
edit_1 变化(缩短模式)的示例。仅显示前几个步骤/较短的模式:
^abcde
^bcde
^cde
或者,如果您更喜欢它作为 DNA:
^GATACCA
^ATACCA
^TACCA
edit_2 简单的 grep 并没有真正削减它。我需要对 FASTQ 格式的每 4 行进行后处理,其中只有第 2 行匹配。如果我不使用 fqgrep,那么我必须:
读取 4 行输入
- 检查第 2 行(序列)是否以 20 种模式中的任何一种开始(P30-P10)
- 如果我得到了匹配,我需要剪掉#2和#4行的前N个字符,其中N代表匹配模式的长度
- 在不匹配的情况下打印输出/写入文件行#1-$4 什么都不做
对于内部解决方案,我可以尝试使用 GNU 并行将输入文件拆分为 4M 的谎言块,并以这种方式加快处理速度。但是,如果我想让每个新软件都能被其他人使用,我会要求最终用户安装广告,这会更加复杂。
** 编辑 3 ** 来自 Vyctor 的正则表达式和匹配行的简单示例:
starting P30 regex
^agacatacagagacatacagagacatacag
matching sequence:
^agacatacagagacatacagagacatacagGAGGACCA
P29:
^gacatacagagacatacagagacatacag
matching sequence:
^gacatacagagacatacagagacatacagGACCACCA
P28:
^acatacagagacatacagagacatacag
matching sequence:
^acatacagagacatacagagacatacagGATTACCA
我从左侧删除字符/DNA 碱基(或 DNA 中的 5 个素数末端),因为这是这些序列被真正的酶降解的方式。一旦被发现,正则表达式序列本身就没有意义了。所需的输出是正则表达式之后的读取序列。在上面的例子中,它是大写的,然后可以在下一步中映射到基因组。 应该强调的是,除了这个玩具示例之外,我正在变得更长,在正则表达式模式之后的先验未知和变化的序列。在现实世界中,我不必处理 DNA 的大写/小写字符(一切都是大写的),但我可能会在我正在搜索模式的序列中遇到 Ns(= 未知 DNA 碱基)。这些可以在第一次近似中忽略,但对于更敏感的算法版本可能应该作为简单的不匹配处理。在理想情况下,我们不会计算给定位置的简单错配,而是考虑存储在以 FASTQ 格式存储的每 4 行长序列记录的第 4 行中的 DNA 序列质量值来计算更复杂的惩罚:http://en.wikipedia.org/wiki/FASTQ_format#Quality
但这要复杂得多,到目前为止,“只读取与正则表达式完全匹配的读取”的方法已经足够好,并且使后续步骤更容易分析。
【问题讨论】:
-
我认为你的问题是非常好的结束语,因为你不清楚你在问什么。请添加一些示例,说明您要如何处理“这是源,这是预期的输出”。
-
如果您添加当前的
sh文件也会有所帮助。并考虑这些很好的答案:stackoverflow.com/questions/9066609/fastest-possible-grep -
只匹配前 20 个字符选项:
^a?g?a?c?a?t?a?c?a?g?... -
我为您发布了一些更新信息。真的,当您只需要 1 个常量时,没有必要经历检查 30 个单独的正则表达式的所有麻烦。例如,在您的 edit_2 中,匹配的长度减去 P30 字符的总长度告诉您要修剪多少。它是一个简单的减法。真的很容易。