加快迭代长字符串中的字符

Posted 2023-03-12

【中文标题】加快迭代长字符串中的字符

【英文标题】：Speed up iterating over characters in long strings

【发布时间】：2022-01-02 16:24:06

【问题描述】：

我有一个采用 DNA 字符串的代码，其中仅找到 4 个字符：A、C、T 和 G，例如“ATACAAG”，并且对于每个字符，如果找到了 3 个其他可能的字符。该代码包括一个用于字符串的循环和另一个用于可能字符列表的循环。问题是字符串很长：多达数十万个字符，因此速度不快并且计算机会发热（风扇开始快速工作）。

我正在寻找一种更快的方法来做到这一点。我尝试了列表理解，但它仍然很慢。我还尝试将代码作为函数从 pandas lambda 调用，每个字符串仍然需要大约一分钟。这是我能得到的最好的吗？

对于每个字符，代码在文件的不同行中记录了 3 个替代项。

代码：

bases = set(list('ACGT'))
alts = base: list(bases.difference(base)) for base in bases

def get_variants(data, output_path):    # pb: position base, b: base
    [open(output_path + f'/data.symbol_variants.txt', 'a').writelines(
        [f'data.chromosome\tdata.end + index\tdata.end + index\tpb/b\tdata.strand\n' for b in alts[pb]])
        for index, pb in enumerate(data.sequence)]

为“ATACAAG”调用函数：

get_variants(pandas.Series('symbol': 'XYZ', 'sequence': 'ATACAAG', 'chromosome': 12, 'start': 9067664, 'end': 9067671, 'strand': '-'),
             'write_an_existing_output_directory_path_here')

输出排列在以下列中的文件中：

chromosome number, start position, end position, original character/alternative character, strand (can + or -)

它在文件 XYZ_variants.txt 中产生以下行：

12  9067664 9067664 A/T -
12  9067664 9067664 A/G -
12  9067664 9067664 A/C -
12  9067665 9067665 T/A -
12  9067665 9067665 T/G -
12  9067665 9067665 T/C -
12  9067666 9067666 A/T -
12  9067666 9067666 A/G -
12  9067666 9067666 A/C -
12  9067667 9067667 C/T -
12  9067667 9067667 C/A -
12  9067667 9067667 C/G -
12  9067668 9067668 A/T -
12  9067668 9067668 A/G -
12  9067668 9067668 A/C -
12  9067669 9067669 A/T -
12  9067669 9067669 A/G -
12  9067669 9067669 A/C -
12  9067670 9067670 G/T -
12  9067670 9067670 G/A -
12  9067670 9067670 G/C -

谢谢。

【问题讨论】：

请提供输入文件/数据框的示例

我有。该问题包括一个pandas.Series 对象，它是数据框中一行的示例。

如果你想遍历字符串中的每个字符，那么不，没有更快的方法。列表推导只会使代码更小，不会加快进程。

【参考方案1】：

我会这样做。

从数据框开始：

  symbol sequence chromosome    start      end strand
0    XYZ  ATACAAG         12  9067664  9067671      -

我会explode 序列，reindex 拥有 A/C/G/T 的所有组合，并且只保留初始碱基不同的组合

import numpy as np

df2 = df.assign(base=df['sequence'].apply(list)).explode('base').reset_index()
df2 = (df2.reindex(df2.index.repeat(4))
          .assign(variant=np.tile(list('ACGT'), len(df2)))
          .loc[lambda d: d['base'].ne(d['variant'])]
          .assign(var=lambda d:d['base']+'/'+d['variant'])
       )

中间输出：

>>> df2.head()
   index symbol sequence chromosome    start      end strand base variant  var
0      0    XYZ  ATACAAG         12  9067664  9067671      -    A       C  A/C
0      0    XYZ  ATACAAG         12  9067664  9067671      -    A       G  A/G
0      0    XYZ  ATACAAG         12  9067664  9067671      -    A       T  A/T
1      0    XYZ  ATACAAG         12  9067664  9067671      -    T       A  T/A
1      0    XYZ  ATACAAG         12  9067664  9067671      -    T       C  T/C

然后导出：

df2[['start', 'end', 'var', 'strand']].to_csv('variants.txt', sep='\t', index=False, header=None)

示例输出（第一行）：

9067664 9067671 A/C -
9067664 9067671 A/G -
9067664 9067671 A/T -
9067664 9067671 T/A -
9067664 9067671 T/C -
9067664 9067671 T/G -
9067664 9067671 A/C -
9067664 9067671 A/G -
9067664 9067671 A/T -
9067664 9067671 C/A -

优化

现在我们删除所有不需要的东西以保持最小尺寸：

df2 = (df.drop(columns=['symbol', 'chromosome'])
         .assign(sequence=df['sequence'].apply(list))
         .explode('sequence').reset_index(drop=True)
      )
df2 = (df2.reindex(df2.index.repeat(4))
          .assign(var=np.tile(list('ACGT'), len(df2)))
          .loc[lambda d: d['sequence'].ne(d['var'])]
          .assign(var=lambda d:d['sequence']+'/'+d['var'])
       )
df2[['start', 'end', 'var', 'strand']].to_csv('variants.txt', sep='\t', index=False, header=None)

【讨论】：

谢谢。考虑到序列有数十万个字符，这是一个好方法吗？我尝试使用您的代码并获得了大约 330 万行的数据框。将来它将在具有大量内存和内核的服务器上运行，因此可能没问题。

好吧，如果你有数百万行，无论如何都需要时间。这里的好处是您的行可能是独立的，因此将输入拆分为较小的数据帧，应用于不同的线程/核心并加入输出是一个完美的案例；）让我知道它是怎么回事！

另外，输出包括第一列的染色体，所以第一行代码为df2 = (df.drop(columns=['symbol'])，第五行代码为df2 = (df2.reindex(df2.index.repeat(4))，最后一行代码为df2[['chromosome', 'start', 'end', 'var', 'strand']]?

似乎正确，在测试样本上试一试；）

我希望将每个符号/基因名称保存在一个单独的文件中，因此我将数据帧拆分为多个数据帧，并分别在每个数据帧上使用您建议的代码。你会采取不同的方式吗？