【问题标题】:Looking for elegant glob-like DNA string expansion寻找优雅的球状 DNA 字符串扩展
【发布时间】:2009-07-08 14:28:03
【问题描述】:

我正在尝试对一组具有多个可能碱基的 DNA 字符串进行球状扩展。

我的 DNA 字符串的碱基包含字母 A、C、G 和 T。但是,我可以有特殊字符,例如 M,可以是 A 或 C。

例如,假设我有字符串:

ATMM

我想把这个字符串作为输入,输出四个可能匹配的字符串:

ATAA ATAC ATCA ATCC

我觉得必须有一些优雅的 Python/Perl/正则表达式技巧才能做到这一点,而不是蛮力解决方案。

感谢您的建议。

编辑,感谢 cortex 的产品运营商。这是我的解决方案:

仍然是 Python 新手,所以我敢打赌,处理每个字典键的方法比另一种 for 循环更好。任何建议都会很棒。

import sys
from itertools import product

baseDict = dict(M=['A','C'],R=['A','G'],W=['A','T'],S=['C','G'],
                  Y=['C','T'],K=['G','T'],V=['A','C','G'],
                  H=['A','C','T'],D=['A','G','T'],B=['C','G','T'])
def glob(str):
    strings = [str]

    ## this loop visits very possible base in the dictionary
    ## probably a cleaner way to do it
    for base in baseDict:
        oldstrings = strings
        strings = []
        for string in oldstrings:
            strings += map("".join,product(*[baseDict[base] if x == base 
                                 else [x] for x in string]))
    return strings

for line in sys.stdin.readlines():
    line = line.rstrip('\n')
    permutations = glob(line)
    for x in permutations:
        print x

【问题讨论】:

    标签: python permutation glob dna-sequence


    【解决方案1】:

    同意其他发帖人的观点,这似乎是一件奇怪的事情。当然,如果你真的想要,在 Python(2.6+)中(一如既往)有一种优雅的方式来做到这一点:

    from itertools import product
    map("".join, product(*[['A', 'C'] if x == "M" else [x] for x in "GMTTMCA"]))
    

    带有输入处理的完整解决方案:

    import sys
    from itertools import product
    
    base_globs = {"M":['A','C'], "R":['A','G'], "W":['A','T'],
                  "S":['C','G'], "Y":['C','T'], "K":['G','T'],
    
                  "V":['A','C','G'], "H":['A','C','T'],
                  "D":['A','G','T'], "B":['C','G','T'],
                  }
    
    def base_glob(glob_sequence):
        production_sequence = [base_globs.get(base, [base]) for base in glob_sequence]
        return map("".join, product(*production_sequence))
    
    for line in sys.stdin.readlines():
        productions = base_glob(line.strip())
        print "\n".join(productions)
    

    【讨论】:

    • 这看起来很有趣,我看看能不能让它做我想做的事。
    • 我想我得到了一些工作。如果您不介意,您能否查看解决方案并让我知道是否有更简洁的方法来处理外部 for 循环?
    • 当然,我添加了完整解决方案的版本,并加入了更多的pythonicisms。顺便说一句:您的 glob 语法中缺少很多排列: from itertools import permutations list(permutations("ACGT", 2)) list(permutations("ACGT", 2))
    • 而不是“base_globs[base] if base in base_globs else [base]”你可以使用“base_globs.get(x, [x])”,它更简洁一点...
    【解决方案2】:

    你可能可以在 python 中使用 yield 操作符来做这样的事情

    def glob(str):
          if str=='':           
              yield ''
              return      
    
          if str[0]!='M':
              for tail in glob(str[1:]): 
                  yield str[0] + tail                  
          else:
             for c in ['A','G','C','T']:
                 for tail in glob(str[1:]):
                     yield c + tail                 
          return
    

    编辑:正如正确指出的那样,我犯了一些错误。这是我试用过的版本。

    【讨论】:

    • 我是 Python 新手,所以对生成器了解不多。我明白了这一点,但我不明白将字符连接到生成器的行。这可能吗?
    • @Rich 它没有将一个字符连接到生成器......生成器每次都会产生连接的结果
    【解决方案3】:

    这并不是一个真正的“扩展”问题,而且几乎可以肯定任何合理的正则表达式都不可行。

    我相信您正在寻找的是“如何生成排列”。

    【讨论】:

    • 我在考虑扩展,你会说 glob 扩展了字符集的匹配。但是,你是对的,它确实感觉更像是一种排列。
    【解决方案4】:

    例如,您可以递归地执行此操作。伪代码:

    printSequences(sequence s)
      switch "first special character in sequence"
        case ...
        case M:
          s1 = s, but first M replaced with A
          printSequences(s1)
          s2 = s, but first M replaced with C
          printSequences(s2)
        case none:
          print s;
    

    【讨论】:

      【解决方案5】:

      正则表达式 match 字符串,它们并不打算变成它们可能匹配的每个字符串。

      此外,您正在查看由此输出的大量字符串 - 例如:

      MMMMMMMMMMMMMMMM (16 M's)
      

      产生 65,536 个 16 个字符串 - 我猜 DNA 序列通常比这更长。

      从计算机科学的角度来看,可以说任何解决方案都是“蛮力”,因为您的算法在原始字符串长度上是 O(2^n)。实际上还有很多工作要做。

      为什么要生成所有组合?你打算怎么处理他们? (如果您想产生每个字符串的可能性,然后在一个大的 DNA 序列中寻找它,那么有很多更好的方法来做到这一点。)

      【讨论】:

      • 幸运的是,我正在查看大约 80 个小的 DNA sn-ps,其中没有一个具有超过 8 个这些生成碱基。扩展整个列表后,它应该很容易适合磁盘。没错,这是将这些字符串与更大的 DNA 序列进行比较的项目的一部分,但我已经通过直接使用压缩格式完成了这项工作。出于政治原因,我现在需要生成整个序列列表。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2013-08-13
      • 2012-08-26
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2013-10-21
      相关资源
      最近更新 更多