6

我有以下任务：从 30 个字符长的模式序列开始（它实际上是 DNA 序列，以免将其称为 P30）我需要在文本文件中找到所有以 (^agacatacag... ) 开头的行，精确的 P30，然后是 29 30、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 行和第 2 行的前 N ​​个字符4，其中 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

但这要复杂得多，到目前为止，“仅读取与正则表达式完全匹配的读取”的方法已经足够好，并且使后续步骤更易于分析。

4

2 回答 2

2

您可以以编程方式生成正则表达式，如下所示。
它只是行开始或下一个字符的渐进交替。

这将为您提供进行单遍搜索的能力。
你所要做的就是得到匹配的字符串长度来告诉你你
在哪里。

注意 - 使用多行模式。 

 #  (?:^|a)(?:^|g)(?:^|a)(?:^|c)(?:^|a)(?:^|t)(?:^|a)(?:^|c)(?:^|a)(?:^|g)(?!^)0123456789

 (?: ^ | a )      # P30
 (?: ^ | g )      # P29
 (?: ^ | a )      # P28
 (?: ^ | c )      # P27
 (?: ^ | a )      # P26
 (?: ^ | t )      # P25
 (?: ^ | a )      # P24
 (?: ^ | c )      # P23
 (?: ^ | a )      # P22
 (?: ^ | g )      # P21
                  # ..
                  # P11
 (?! ^)           # Not beginning of line
 0123456789       # P10 - P1

例如，匹配这些： 

agacatacag0123456789
cag0123456789
gacatacag0123456789
acatacag0123456789
acag0123456789
catacag0123456789
catacag0123456789
0123456789
atacag0123456789
tacag0123456789
ag0123456789
g0123456789

但不是这些： 

agaatacag0123456789
ca0123456789
gacataca0123456789
acaacag0123456789
acg0123456789
cataca0123456789
caacag0123456789
123456789
atacg0123456789
tcag0123456789
ag012456789
g012356789

更新

这是一个图形说明，单个正则表达式可以替换所有 30 个。
实际上不需要 30 个单独的正则表达式，您只需要 1 个常量正则表达式。
在此示例中，集群组被替换为捕获组，因此您可以看到它在做什么。 

 # (^|G)(^|A)(^|T)(^|A)(^|C)(^|C)(?!^)A

 ( ^ | G )        # (1)
 ( ^ | A )        # (2)
 ( ^ | T )        # (3)
 ( ^ | A )        # (4)
 ( ^ | C )        # (5)
 ( ^ | C )        # (6)
 (?! ^ )          # Not beginning of line
 A

输入，6 行： 

GATACCA
ATACCA
TACCA
ACCA
CCA
CA

输出： 

 **  Grp 0 -  ( pos 0 , len 7 ) 
GATACCA  
 **  Grp 1 -  ( pos 0 , len 1 ) 
G  
 **  Grp 2 -  ( pos 1 , len 1 ) 
A  
 **  Grp 3 -  ( pos 2 , len 1 ) 
T  
 **  Grp 4 -  ( pos 3 , len 1 ) 
A  
 **  Grp 5 -  ( pos 4 , len 1 ) 
C  
 **  Grp 6 -  ( pos 5 , len 1 ) 
C  

------------------------------

 **  Grp 0 -  ( pos 9 , len 6 ) 
ATACCA  
 **  Grp 1 -  ( pos 9 , len 0 )  EMPTY 
 **  Grp 2 -  ( pos 9 , len 1 ) 
A  
 **  Grp 3 -  ( pos 10 , len 1 ) 
T  
 **  Grp 4 -  ( pos 11 , len 1 ) 
A  
 **  Grp 5 -  ( pos 12 , len 1 ) 
C  
 **  Grp 6 -  ( pos 13 , len 1 ) 
C  

------------------------------

 **  Grp 0 -  ( pos 17 , len 5 ) 
TACCA  
 **  Grp 1 -  ( pos 17 , len 0 )  EMPTY 
 **  Grp 2 -  ( pos 17 , len 0 )  EMPTY 
 **  Grp 3 -  ( pos 17 , len 1 ) 
T  
 **  Grp 4 -  ( pos 18 , len 1 ) 
A  
 **  Grp 5 -  ( pos 19 , len 1 ) 
C  
 **  Grp 6 -  ( pos 20 , len 1 ) 
C  

------------------------------

 **  Grp 0 -  ( pos 24 , len 4 ) 
ACCA  
 **  Grp 1 -  ( pos 24 , len 0 )  EMPTY 
 **  Grp 2 -  ( pos 24 , len 0 )  EMPTY 
 **  Grp 3 -  ( pos 24 , len 0 )  EMPTY 
 **  Grp 4 -  ( pos 24 , len 1 ) 
A  
 **  Grp 5 -  ( pos 25 , len 1 ) 
C  
 **  Grp 6 -  ( pos 26 , len 1 ) 
C  

------------------------------

 **  Grp 0 -  ( pos 30 , len 3 ) 
CCA  
 **  Grp 1 -  ( pos 30 , len 0 )  EMPTY 
 **  Grp 2 -  ( pos 30 , len 0 )  EMPTY 
 **  Grp 3 -  ( pos 30 , len 0 )  EMPTY 
 **  Grp 4 -  ( pos 30 , len 0 )  EMPTY 
 **  Grp 5 -  ( pos 30 , len 1 ) 
C  
 **  Grp 6 -  ( pos 31 , len 1 ) 
C  

------------------------------

 **  Grp 0 -  ( pos 35 , len 2 ) 
CA  
 **  Grp 1 -  ( pos 35 , len 0 )  EMPTY 
 **  Grp 2 -  ( pos 35 , len 0 )  EMPTY 
 **  Grp 3 -  ( pos 35 , len 0 )  EMPTY 
 **  Grp 4 -  ( pos 35 , len 0 )  EMPTY 
 **  Grp 5 -  ( pos 35 , len 0 )  EMPTY 
 **  Grp 6 -  ( pos 35 , len 1 ) 
C
于 2014-12-12T22:09:48.777 回答
2

如果我对您的理解正确，则您有要匹配的预设：

agacatacagagacatacagagacatacag

然后匹配匹配的行：

re: agacatacagagacatacagagacatacag
30: agacatacagagacatacagagacatacag
29: agacatacagagacatacagagacatacac
28: agacatacagagacatacagagacataccc

您实际上并不需要正则表达式，您只需要找到两行之间的差异，因为它是 DNA，我假设该字符串abcde并且aebcd差异为 4，因为所有序列都需要位于正确的位置。

如果顺序无关紧要，并且您只想搜索至少匹配 28 个字符的行，您可以只计算字符串中的差异。

reg = 'agacatacagagacatacagagacatacag'
for row in file:
    letters = diff_letters(reg, row.strip())
    if letters == 30:   # complete match
    elif letters == 29: # one different character
                        # so on

如果您需要实际上以正确顺序开始的匹配，您可以获取字符串之间的差异点，如果第一个差异在点>=28

reg = 'agacatacagagacatacagagacatacag'
for row in file:
    diffs = list(i for i,(a1,a2) in enumerate(zip(s1,s2)) if a1!=a2)
    if not len(diffs):
        difference = len(reg)
    else:
        difference = diffs[0]

    if difference == 30: # First difference is at last offset
于 2014-12-12T22:11:45.950 回答