Python - 检查两个大文本文件之间的一致性

Posted 2023-04-18

【中文标题】Python - 检查两个大文本文件之间的一致性

【英文标题】：Python - Checking concordance between two huge text files

【发布时间】：2015-11-18 07:50:28

【问题描述】：

所以，这个让我很难受！ 我正在处理 HUGE 文本文件，我的意思是 100Gb+。具体来说，它们位于fastq format 中。此格式用于 DNA 测序数据，由四行记录组成，如下所示：

@REC1
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT  
+  
!''*((((***+))%%%++)(%%%%).1***-+*''))*55CCF>>>>>>CCCCCCC65  
@REC2
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT  
+  
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65  
.  
.  
.

为了这个问题，只关注标题行，以“@”开头。

因此，出于 QA 的目的，我需要比较两个这样的文件。这些文件应该有匹配的标题，所以另一个文件中的第一条记录也应该有标题'@REC1'，下一个应该是'@REC2'等等。在进行大量下游分析之前，我想确保情况确实如此。 由于文件很大，一个简单的迭代字符串比较会花费很长时间，但是这个 QA 步骤会运行很多次，我不能等那么久。所以我认为更好的方法是从文件中的几个点采样记录，例如每 10% 的记录。如果记录的顺序搞砸了，我很可能会发现它。 到目前为止，我已经能够通过估计文件大小来处理此类文件，而不是使用 python 的file.seek() 访问文件中间的记录。例如，要访问大约在中间的一行，我会这样做：

file_size = os.stat(fastq_file).st_size
start_point = int(file_size/2)
with open(fastq_file) as f:
    f.seek(start_point)
    # look for the next beginning of record, never mind how

但现在问题更复杂了，因为我不知道如何在两个文件之间进行协调，因为字节位置不是文件中行索引的指示符。换句话说，我如何访问两个文件中的第 10,567,311 行以确保它们相同，而无需遍历整个文件？

将不胜感激任何想法\提示。也许并行迭代？但具体如何？ 谢谢！

【问题讨论】：

我缩进了您的文件示例以防止 SO 将其格式化为粗体/斜体等 - 我希望结果是正确的。请检查我是否搞砸了。

请求澄清：如果相应的@REC123 行出现在两个文件中的相同行号，您将认为两个文件一致。没有其他标准？

@TimPietzcker - 感谢您的编辑，是的，这是唯一的标准。很简单...

顺便说一句，要使质量控制成功，不应该允许任何差异，对吧？换句话说，我们想在第一个错误时中止？

你又是对的。

【参考方案1】：

抽样是一种方法，但您要靠运气。此外，Python 不是这项工作的错误工具。您可以使用标准的 Unix 命令行工具以不同的方式做事并以仍然相当有效的方式计算出准确的答案：

线性化您的 FASTQ 记录：将前三行中的换行符替换为制表符。 在一对线性化文件上运行diff。如有差异，diff 会报告。

要线性化，您可以通过awk 运行您的 FASTQ 文件：

$ awk '\
    BEGIN  \
      n = 0; \
     \
     \
      a[n % 4] = $0; \
      if ((n+1) % 4 == 0)  \
        print a[0]"\t"a[1]"\t"a[2]"\t"a[3]; \
       \
      n++; \
    ' example.fq > example.fq.linear 

比较一对文件：

$ diff example_1.fq.linear example_2.fq.linear

如果有任何差异，diff 会找到并告诉您哪条 FASTQ 记录不同。

您可以直接在这两个文件上运行diff，而无需进行额外的线性化工作，但如果您首先进行线性化，则更容易看出哪个读取有问题。

所以这些都是大文件。写入新文件在时间和磁盘空间上都是昂贵的。有一种方法可以改进这一点，使用 streams。

如果您将awk 脚本放入一个文件（例如、linearize_fq.awk），您可以像这样运行它：

$ awk -f linearize_fq.awk example.fq > example.fq.linear

这可能对您的 100+ Gb 文件有用，因为您现在可以通过 bash process substitutions 设置两个 Unix 文件流，并直接在这些流上运行 diff：

$ diff <(awk -f linearize_fq.awk example_1.fq) <(awk -f linearize_fq.awk example_2.fq)

或者你可以使用named pipes:

$ mkfifo example_1.fq.linear
$ mkfifo example_2.fq.linear
$ awk -f linearize_fq.awk example_1.fq > example_1.fq.linear &
$ awk -f linearize_fq.awk example_2.fq > example_2.fq.linear &
$ diff example_1.fq.linear example_2.fq.linear
$ rm example_1.fq.linear example_2.fq.linear

命名管道和进程替换都避免了创建额外（常规）文件的步骤，这可能是您的输入类型的问题。将 100+ Gb 文件的线性化副本写入磁盘可能需要一段时间，而且这些副本也可能会使用您可能没有太多的磁盘空间。

使用流解决了这两个问题，这使得它们对于以有效的方式处理生物信息学数据集非常有用。

您可以使用 Python 重现这些方法，但几乎可以肯定它的运行速度要慢得多，因为 Python 在此类 I/O 繁重的任务中非常慢。

【讨论】：

谢谢！在我标记为答案之前，我会尝试并与其他方法进行比较。

【参考方案2】：

并行迭代可能是在 Python 中执行此操作的最佳方式。我不知道它的运行速度有多快（快速的 SSD 可能是加快速度的最佳方式），但由于无论如何你都必须在两个文件中计算换行符，我看不出有什么办法：

with open(file1) as f1, open(file2) as f2:
    for l1, l2 in zip(f1,f2):
        if l1.startswith("@REC"):
            if l1 != l2:
                print("Difference at record", l1)
                break
    else:
        print("No differences")

这是为 Python 3 编写的，其中 zip 返回一个迭代器；在 Python 2 中，您需要改用 itertools.izip()。

【讨论】：

谢谢，但这是我试图避免的那种幼稚的解决方案。我会尝试这里建议的其他方法，看看哪种方法效果最好。

当然，我明白这一点。我只是怀疑没有办法解决这个问题。无论如何，您都必须读取每个文件中的每一行（否则无法找出当前行号），并且由于您无法将它们读入内存，因此您别无选择，只能并行迭代......除非我错过了一些东西。

@soungalo 同意 TimPietzcker，所有三个答案都建议您逐行读取文件。

【参考方案3】：

您是否考虑过使用rdiff 命令。 rdiff 的优点是：

对于相同的 4.5GB 文件，rdiff 只占用了大约 66MB 的 RAM，并且扩展性非常好。迄今为止，它从未崩溃过。 它也比 diff 快得多。 rdiff 本身结合了差异和补丁功能，因此您可以创建增量并使用相同的程序应用它们

rdiff 的缺点是：

它不是标准 Linux/UNIX 发行版的一部分——您必须 安装 librsync 包。 rdiff 生成的 delta 文件与 diff 的格式略有不同。 delta 文件稍大（但不够重要）。 使用 rdiff 生成 delta 时使用了一种略有不同的方法，这有好有坏 - 需要 2 个步骤。这 第一个生成一个特殊的签名文件。第二步，一个 delta 是使用另一个 rdiff 调用创建的（如下所示）。尽管 两步过程可能看起来很烦人，它的好处是 提供比使用 diff 时更快的增量。

见：http://beerpla.net/2008/05/12/a-better-diff-or-what-to-do-when-gnu-diff-runs-out-of-memory-diff-memory-exhausted/

【讨论】：

这看起来很有趣，感谢您发布此答案。

【参考方案4】：

import sys
import re

""" To find of the difference record in two HUGE files. This is expected to  
use of minimal memory. """

def get_rec_num(fd):
    """ Look for the record number. If not found return -1"""
    while True:
        line = fd.readline()
        if len(line) == 0: break
        match =  re.search('^@REC(\d+)', line)
        if match:
            num = int(match.group(1))
            return(num)
    return(-1)

f1 = open('hugefile1', 'r')
f2 = open('hugefile2', 'r')

hf1 = dict()
hf2 = dict()
while f1 or f2:
    if f1:
        r = get_rec_num(f1)
        if r < 0:
            f1.close()
            f1 = None
        else:
            # if r is found in f2 hash, no need to store in f1 hash  
            if not r in hf2:
                hf1[r] = 1
            else:
                del(hf2[r])
        pass
    pass
    if f2:
        r = get_rec_num(f2)
        if r < 0:
            f2.close()
            f2 = None
        else:
            # if r is found in f1 hash, no need to store in f2 hash  
            if not r in hf1:
                hf2[r] = 1
            else:
                del(hf1[r])
        pass
    pass

print('Records found only in f1:')
for r in hf1:
    print(', '.format(r));
print('Records found only in f2:')
for r in hf2:
    print(', '.format(r));

【讨论】：

谢谢！在我标记为答案之前，我会尝试并与其他方法进行比较。

【参考方案5】：

从我的角度来看，@AlexReynolds 和 @TimPietzcker 的两个答案都非常出色，但我想投入两分钱。您可能还想加快硬件速度：

带 SSD 的 Raplace 硬盘 使用n SSD 并创建一个 RAID 0。在理想情况下，您的磁盘 IO 将获得n 倍的速度。 调整从 SSD/HDD 读取的块大小。例如，我希望执行一次 16 MB 的读取比执行十六次 1 MB 的读取要快。 （这适用于单个 SSD，对于 RAID 0 优化，必须查看 RAID 控制器选项和功能）。

最后一个选项与 NOR SSD 尤其相关。不要追求最小的 RAM 使用率，而是尽可能多地读取以保持磁盘快速读取。例如，从两个文件中并行读取单行可能会降低读取速度 - 想象一个 HDD，其中两个文件的两行始终位于同一个磁盘的同一侧。