2

所以，这个让我很难受！
我正在处理巨大的文本文件，我的意思是 100Gb+。具体来说，它们采用fastq 格式。这种格式用于 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 行以确保它们相同，而无需遍历整个文件？

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

4

5 回答 5

3

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

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

要进行线性化，您可以通过以下方式运行 FASTQ 文件awk：

$ 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在这两个文件上运行，而无需进行额外的线性化工作，但是如果您首先进行线性化，则更容易看出哪个读取有问题。

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

如果将awk脚本放入文件（例如, linearize_fq.awk）中，则可以像这样运行它：

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

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

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

或者您可以使用命名管道：

$ 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 繁重的任务上非常慢。

于 2015-11-18T08:07:35.043 回答
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()。

于 2015-11-18T08:52:26.033 回答
1

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

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

rdiff 的缺点是：

  • 它不是标准 Linux/UNIX 发行版的一部分——您必须安装 librsync 包。
  • rdiff 生成的 delta 文件的格式与 diff 的格式略有不同。
  • 增量文件稍大（但不足以关心）。
  • 使用 rdiff 生成 delta 时使用了一种稍微不同的方法，这有好有坏——需要 2 个步骤。第一个生成一个特殊的签名文件。在第二步中，使用另一个 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/

于 2015-11-18T10:10:30.953 回答
0

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

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

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

于 2015-11-18T09:41:05.273 回答
0 
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));
于 2015-11-18T09:38:14.723 回答