python - 如何更有效地读取 DNA 序列？

Question

我在 python 中编写了一个代码来读取 DNA 序列(稍后对它们进行基序比对)但是，我正在寻找一种更有效的方法来执行此操作。

如果您能提供帮助，请参阅下文:

handle = open("a.fas.txt", "r")
a = handle.readlines()[1:]
a = ''.join([x.strip() for x in a])
with open("Output.txt", "w") as text_file:
    text_file.write(a)

f = 0
z = 100
b = ''
while f < len(a):
    b += a[f:z]+'\n'
    f += 1
    z += 1
with open("2.txt", "w") as runner_mtfs:
   runner_mtfs.write(b)

总而言之，我想对 b 的每一行进行大量分析，但我不知道有什么更有效的方法来执行此操作，而不是将每 100 个碱基对分开。输出文件超过 500 兆字节。有什么建议吗？

代码的第一部分只是一个 DNA 序列，我将所有的线连接在一起，并分离 100 个碱基对。

Answer 1

我在这里看到的主要问题是您将所有内容都写到一个文件中。这样做没有意义。您创建的大型输出文件非常冗余，在您进行分析时将其重新加载没有帮助。

最初加载文件后，您有兴趣查看的每个窗口都是 a[x:x+100] 一些 x。您根本不需要实际显式生成这些窗口:这样做应该没有任何好处。遍历并直接从 a 的每个窗口生成这些矩阵。

如果您真的需要整个东西，请将其生成为一个 numpy 数组。此外，如果我不使用任何退化的基本代码，请使用 0、1、2、3 表示 A、C、G、T 将序列转换为 uint8。这有助于加快速度，尤其是当您需要采取在任何一点补码，这可以通过简单地摆弄位来完成。

Numpy 可以使用 stride_tricks 非常有效地生成数组，如前所述 in this blog post :

def rolling_window(a, window):
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
handle = open("U00096.2.fas.txt", "r")
a = handle.readlines()[1:]
a = ''.join([x.strip() for x in a])
b = numpy.array([x for x in a], dtype=numpy.character)
rolling_window(b,100)

或者，转换为整数:

def rolling_window(a, window):
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
handle = open("U00096.2.fas.txt", "r")
a = handle.readlines()[1:]
a = ''.join([x.strip() for x in a])
conv = {'a': 0, 'c': 1, 'g': 2, 't': 3}
b = numpy.array([conv[x] for x in a], dtype=numpy.uint8)
rolling_window(b,100)

这段代码在我的机器上比你的代码快大约十倍。