【问题标题】:Convert String to 64bit integer mapping characters to custom two-bit values mapping将字符串转换为 64 位整数映射字符到自定义两位值映射
【发布时间】:2018-10-18 14:54:09
【问题描述】:

我正在尝试将一串字符(A、T、C、G)映射为一个 64 位整数,其中每个字母使用此映射表示为两位:

mapping = {'A': 0b00, 'C': 0b01, 'G': 0b10, 'T': 0b11} 

“sequence”字符串不会超过 28 个字符,并且我打算将整数在开头补零以使其成为 64 位。目前,我使用以下功能,但速度非常慢。然后我通过调用转换输出:

int(result, 2)

这目前有效,但我想让这个功能非常快。我不太了解 C++,所以我很难移植到它。我现在正在尝试 Cython,但我也不熟悉。任何有助于在 Python(甚至 C++ 或 Cython 等价物)中提高效率的帮助将不胜感激。

下面是我的代码,之后我再次调用 int()。

def seq_to_binary(seq):
    values = [mapping[c] for c in seq]
    BITWIDTH = 2
    return "".join(map(lambda x: bin(x)[2:].zfill(BITWIDTH), values)).encode();

在典型的序列输入中会是这样的:'TGTGAGAAGCACCATAAAAGGCGTTGTG'

【问题讨论】:

    标签: python python-3.x binary cython


    【解决方案1】:

    您正在将一个由 4 个不同“数字”组成的字符串解释为一个数字,因此 base 4 表示法。如果你有一串实际数字,在 0-3 范围内,你可以让 int() 非常快地生成一个整数。

    def seq_to_int(seq, _m=str.maketrans('ACGT', '0123')):
        return int(seq.translate(_m), 4)
    

    上述函数使用str.translate() 将4 个字符中的每一个替换为一个匹配的数字(我使用静态str.maketrans() function 创建翻译表)。然后将生成的数字字符串解释为以 4 为底的整数。

    请注意,这会产生一个整数对象,而不是零和一个字符的二进制字符串:

    >>> seq_to_int('TGTGAGAAGCACCATAAAAGGCGTTGTG')
    67026852874722286
    >>> format(seq_to_int('TGTGAGAAGCACCATAAAAGGCGTTGTG'), '016x')
    '00ee20914c029bee'
    >>> format(seq_to_int('TGTGAGAAGCACCATAAAAGGCGTTGTG'), '064b')
    '0000000011101110001000001001000101001100000000101001101111101110'
    

    这里不需要填充;只要您的输入序列是 32 个字母或更少,生成的整数将适合无符号 8 字节整数表示。在上面的输出示例中,我使用 format() 字符串将该整数值分别格式化为十六进制和二进制字符串,并将这些表示零填充为 64 位数字的正确位数。

    为了衡量这是否更快,让我们抽取 100 万个随机生成的测试字符串(每个 28 个字符长):

    >>> from random import choice
    >>> testvalues = [''.join([choice('ATCG') for _ in range(28)]) for _ in range(10 ** 6)]
    

    在我的 Macbook Pro 上使用 2.9 GHz Intel Core i7,在 Python 3.6.5 上,上述函数可以在不到 3/4 秒的时间内产生 100 万次转换:

    >>> from timeit import timeit
    >>> timeit('seq_to_int(next(tviter))', 'from __main__ import testvalues, seq_to_int; tviter=iter(testvalues)')
    0.7316284350017668
    

    所以每次调用需要 0.73 微秒。

    (之前我提倡预计算版本,但经过实验后,我想到了 base-4 的想法)。

    为了与迄今为止发布的其他方法进行比较,有些方法也需要进行调整以产生整数,并被包装到函数中:

    def seq_to_int_alexhall_a(seq, mapping={'A': b'00', 'C': b'01', 'G': b'10', 'T': b'11'}):
        return int(b''.join(map(mapping.__getitem__, seq)), 2)
    
    def seq_to_int_alexhall_b(seq, mapping={'A': b'00', 'C': b'01', 'G': b'10', 'T': b'11'}):
        return int(b''.join([mapping[c] for c in seq]), 2)
    
    def seq_to_int_jonathan_may(seq, mapping={'A': 0b00, 'C': 0b01, 'G': 0b10, 'T': 0b11}):
        result = 0
        for char in seq:
            result = result << 2
            result = result | mapping[char]
        return result
    

    然后我们可以比较这些:

    >>> testfunctions = {
    ...     'Alex Hall (A)': seq_to_int_alexhall_a,
    ...     'Alex Hall (B)': seq_to_int_alexhall_b,
    ...     'Jonathan May': seq_to_int_jonathan_may,
    ...     # base_decode as defined in https://stackoverflow.com/a/50239330
    ...     'martineau': base_decode,
    ...     'Martijn Pieters': seq_to_int,
    ... }
    >>> setup = """\
    ... from __main__ import testvalues, {} as testfunction
    ... tviter = iter(testvalues)
    ... """
    >>> for name, f in testfunctions.items():
    ...     res = timeit('testfunction(next(tviter))', setup.format(f.__name__))
    ...     print(f'{name:>15}: {res:8.5f}')
    ...
      Alex Hall (A):  2.17879
      Alex Hall (B):  2.40771
       Jonathan May:  3.30303
          martineau: 16.60615
    Martijn Pieters:  0.73452
    

    我建议的 base-4 方法很容易在这个比较中获胜。

    【讨论】:

      【解决方案2】:

      我在 Cython 中笨拙的直接尝试,它的速度是迄今为止最佳解决方案(@MartijnPieters 的)的两倍:

      %%cython
      
      ctypedef unsigned long long ull
      
      cdef ull to_int(unsigned char *data, int n):
          cdef ull res=0
          cdef int i
          cdef unsigned char ch
          for i in range(n):
              res<<=2
              ch=data[i]
              if ch==67: #C
                  res+=1
              if ch==71: #G
                  res+=2
              if ch==84: #T
                  res+=3
          return res
      
      cpdef str_to_int_ead(str as_str):
          s=as_str.encode('ascii')
          return to_int(s, len(s))
      

      与当前@MartijnPieters 的解决方案相比,它在我的机器上快了一倍:

      >>> [str_to_int_ead(x) for x in testvalues] == [seq_to_int(x) for x in testvalues]
      True
      
      >>> tviter=iter(testvalues)
      >>> %timeit -n1000000 -r1 seq_to_int(next(tviter))
      795 ns ± 0 ns per loop (mean ± std. dev. of 1 run, 1000000 loops each)
      
      >>> tviter=iter(testvalues)
      >>> %timeit -n1000000 -r1 str_to_int_ead(next(tviter))
      363 ns ± 0 ns per loop (mean ± std. dev. of 1 run, 1000000 loops each)
      

      这使得整个运行时间为 0.795 秒,而整个运行时间为 0.363 秒(因此可以与@MartijnPieters 测量的时间进行比较)。

      有人可能会问,如果不需要转换 unicode ascii,可以节省多少开销?

      %%cython
      ....
      cpdef bytes_to_int_ead(bytes as_bytes):
          return to_int(as_bytes, len(as_bytes))
      
      
      >>> testbytes=[bytes(x.encode('ascii')) for x in testvalues]
      >>> tviter=iter(testbytes)
      >>> %timeit -n1000000 -r1 bytes_to_int_ead(next(tviter))
      327 ns ± 0 ns per loop (mean ± std. dev. of 1 run, 1000000 loops each)
      

      仅快 10% - 这有点令人惊讶...

      但是,我们不应该忘记,我们还测量了“下一个”迭代器的开销,而没有得到:

      >>> v=testvalues[0]
      >>> %timeit str_to_int_ead(v)
      >>> 139 ns ± 0.628 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
      
      
      >>> v=testbytes[0]
      >>> %timeit bytes_to_int_ead(v)
      97.2 ns ± 1.03 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
      

      现在实际上有大约 40% 的加速!

      另一个有趣的结论:使用迭代器进行测试时也有大约 250ns(或 70%)的开销。如果没有这个开销,cython 会超过 @MartijnPieters 的 140ns 和 550ns,即几乎高出 4 倍。


      与 cython 进行比较的列表函数(@MartijnPieters 回答的当前状态):

      def seq_to_int(seq, _m=str.maketrans('ACGT', '0123')):
          return int(seq.translate(_m), 4)
      

      测试数据:

      from random import choice
      testvalues = [''.join([choice('ATCG') for _ in range(28)]) for _ in range(10 ** 6)]
      

      【讨论】:

        【解决方案3】:
        seq = 'TGTGAGAAGCACCATAAAAGGCGTTGTG'
        
        mapping = {'A': b'00', 'C': b'01', 'G': b'10', 'T': b'11'}
        
        result = b''.join(map(mapping.__getitem__, seq)).zfill(64)
        
        print(result)
        

        这是一些比较选项的计时代码:

        import timeit
        
        setup = """
        seq = 'TGTGAGAAGCACCATAAAAGGCGTTGTG'
        
        mapping = {'A': b'00', 'C': b'01', 'G': b'10', 'T': b'11'}
        """
        
        for stmt in [
            "b''.join(map(mapping.__getitem__, seq)).zfill(64)",
            "b''.join([mapping[c] for c in seq]).zfill(64)",
        ]:
            print(stmt)
            print(timeit.timeit(stmt, setup, number=10000000))
        

        我发现这两个选项大致相同,但你的结果可能会有所不同。

        【讨论】:

          【解决方案4】:

          使用位移运算符和加法。使用字典来保存字符代码是正确的想法:

          mapping = {'A': 0b00, 'C': 0b01, 'G': 0b10, 'T': 0b11}

          为这个例子生成一个 28 个字符的字符串(这样称呼它有点多余,字符串就可以了):

          chars = 'TGTGAGAAGCACCATAAAAGGCGTTGTG'

          定义一个结果并将其设置为零:

          result = 0

          Python 中的字符串实际上只是一个字符数组,您可以像遍历任何数组一样遍历字符串。我们将使用它以及一系列嵌套的位操作来完成您需要的操作:

          for char in chars:
            result = result << 2
            result = result | mapping[char]
          

          这将产生长度为 2*len(chars) 的位,在本例中为 56。获得额外的

          要添加额外的 8 位前导零,其整数表示实际上是一个 QWORD(64 位),并将自动用零填充 8 个最高有效位。

          print(result)
          >> 67026852874722286
          

          如果你想变得更花哨,你可以使用ctypes 来加速你的代码。

          【讨论】:

            【解决方案5】:

            思考这个问题的一种方法是意识到它所做的事情的本质是从以 4 为底的数字转换为以 10 为底的数字。这可以通过多种方式完成,但我喜欢的一种实际上是非常问题Base 62 conversion的通用接受答案。

            以下是它的修改版本,默认进行 base 4 转换:

            def base_decode(astring, alphabet="ACGT"):
                """Decode a Base X encoded astring into the number
            
                Arguments:
                - `astring`: The encoded astring
                - `alphabet`: The alphabet to use for encoding
                """
                base = len(alphabet)
                strlen = len(astring)
                num = 0
                for idx, char in enumerate(astring):
                    power = (strlen - (idx + 1))
                    num += alphabet.index(char) * (base ** power)
            
                return num
            
            seq = 'TGTGAGAAGCACCATAAAAGGCGTTGTG'
            print('seq_to_binary:', seq_to_binary(seq))
            print('base_decode:', format(base_decode(seq), 'b'))
            

            请注意,这实际上返回一个整数,该整数需要任何位长度(整数在 Python 中是可变长度),以将给定的数字存储为打包成二进制整数值的字符串。添加的对format() 的调用将该值转换为二进制字符串,以便可以打印它并与调用返回字符串的seq_to_binary() 函数的结果进行比较,不是提到的64位整数在标题中。

            【讨论】:

            • Base-4 转换可以完成得非常快; int() 可以做到这一点对我们来说是 C,所需要的只是从字母到数字的字符串转换。
            • @Martijn:意识到并指出这只是一个 base-4 转换问题是我回答的重点,而不是它提出的特定实现——我之所以选择它是因为它方便、容易理解并轻松适应进行base-4转换。将字符串的字符映射到“普通”数字并使用内置的int() 是非常聪明且更快的实现。恭喜。
            猜你喜欢
            • 2012-01-23
            • 1970-01-01
            • 1970-01-01
            • 1970-01-01
            • 1970-01-01
            • 2012-08-21
            • 2018-07-09
            • 2014-05-26
            • 2017-01-16
            相关资源
            最近更新 更多