0

我在这方面花了太多时间，正在寻找建议。我的文件太大（对于感兴趣的人，来自 Illumina 测序运行的 FASTQ 文件）。我需要做的是匹配两个文件之间共有的模式，并将该行加上它下面的 3 行打印到两个单独的文件中，没有重复（存在于原始文件中）。Grep 做得很好，但文件大约 18GB，它们之间的匹配速度非常慢。我需要做的例子如下。

档案A：

@DLZ38V1_0262:8:1101:1430:2087#ATAGCG/1
NTTTCAGTTAGGGCGTTTGAAAACAGGCACTCCGGCTAGGCTGGTCAAGG
+DLZ38V1_0262:8:1101:1430:2087#ATAGCG/1
BP\cccc^ea^eghffggfhh`bdebgfbffbfae[_ffd_ea[H\_f_c
@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/1
NAGGATTTAAAGCGGCATCTTCGAGATGAAATCAATTTGATGTGATGAGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/1
BP\ccceeggggfiihihhiiiihiiiiiiiiihighiighhiifhhhic
@DLZ38V1_0262:8:2316:21261:100790#ATAGCG/1
TGTTCAAAGCAGGCGTATTGCTCGAATATATTAGCATGGAATAATAGAAT
+DLZ38V1_0262:8:2316:21261:100790#ATAGCG/1
__\^c^ac]ZeaWdPb_e`KbagdefbZb[cebSZIY^cRaacea^[a`c

您可以看到 3 个独特的标题，@以 3 行开头

文件B：

@DLZ38V1_0262:8:1101:1430:2087#ATAGCG/2
GAAATCAATGGATTCCTTGGCCAGCCTAGCCGGAGTGCCTGTTTTCAAAC
+DLZ38V1_0262:8:1101:1430:2087#ATAGCG/2
_[_ceeeefffgfdYdffed]e`gdghfhiiihdgcghigffgfdceffh
@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
GCCATTCAGTCCGAATTGAGTACAGTGGGACGATGTTTCAAAGGTCTGGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
_aaeeeeegggggiiiiihihiiiihgiigfggiighihhihiighhiii
@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
GCCATTCAGTCCGAATTGAGTACAGTGGGACGATGTTTCAAAGGTCTGGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
_aaeeeeegggggiiiiihihiiiihgiigfggiighihhihiighhiii
@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
GCCATTCAGTCCGAATTGAGTACAGTGGGACGATGTTTCAAAGGTCTGGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
_aaeeeeegggggiiiiihihiiiihgiigfggiighihhihiighhiii

这里有 4 个标题，但只有 2 个是唯一的，因为其中一个重复了 3 次

我需要两个文件之间没有重复的公共标题加上它们下面的 3 行。每个文件中的顺序相同。

这是我到目前为止所拥有的：

grep -E @DLZ38V1.*/ --only-matching FileA | sort -u -o FileA.sorted
grep -E @DLZ38V1.*/ --only-matching FileB | sort -u -o FileB.sorted
comm -12 FileA.sorted FileB.sorted > combined

结合

@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/
@DLZ38V1_0262:8:1101:1430:2087#ATAGCG/

这只是两个文件之间的公共标题，没有重复。这就是我要的。现在我需要将这些标题与原始文件相匹配，并抓住它们下面的 3 行，但只有一次。

如果我使用 grep 我可以得到我想要的每个文件

while read -r line; do
   grep -A3 -m1 -F $line FileA
done < combined > FileA.Final

文件A.Final

@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/1
NAGGATTTAAAGCGGCATCTTCGAGATGAAATCAATTTGATGTGATGAGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/1
BP\ccceeggggfiihihhiiiihiiiiiiiiihighiighhiifhhhic
@DLZ38V1_0262:8:1101:1430:2087#ATAGCG/1
NTTTCAGTTAGGGCGTTTGAAAACAGGCACTCCGGCTAGGCTGGTCAAGG
+DLZ38V1_0262:8:1101:1430:2087#ATAGCG/1
BP\cccc^ea^eghffggfhh`bdebgfbffbfae[_ffd_ea[H\_f_c

重复while循环以生成FileB.Final

@DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
GCCATTCAGTCCGAATTGAGTACAGTGGGACGATGTTTCAAAGGTCTGGC
+DLZ38V1_0262:8:1101:1369:2106#ATAGCG/2
_aaeeeeegggggiiiiihihiiiihgiigfggiighihhihiighhiii
@DLZ38V1_0262:8:1101:1430:2087#ATAGCG/2
GAAATCAATGGATTCCTTGGCCAGCCTAGCCGGAGTGCCTGTTTTCAAAC
+DLZ38V1_0262:8:1101:1430:2087#ATAGCG/2

这可行，但 FileA 和 FileB 约为 18GB，而我的组合文件约为 2GB。有人对我如何显着加快最后一步有任何建议吗？

4

2 回答 2

1

根据您需要多久运行一次：

  1. 您可以将您的数据转储（您可能希望使用之后构建的索引进行批量插入）到 Postgres（sqlite？）数据库中，在其上构建索引，并享受 40 年来对关系数据库的有效实现的研究成果几乎没有你的投资。

  2. 您可以通过使用 unix 实用程序“join”来模拟拥有一个关系数据库，但不会有太多乐趣，因为这不会给您一个索引，但它可能比“grep”更快，您可能会遇到物理限制...我从未尝试加入两个 18G 文件。

  3. 你可以写一点 C 代码（把你最喜欢的编译（到机器代码）语言放在这里），它将你的字符串（只有四个字母，对吗？）转换为二进制并基于它建立一个索引（或更多）。由于您的 50 个字符串仅占用两个 64 位字，因此这可以实现闪电般的快速和较小的内存占用。

于 2014-05-21T12:43:39.073 回答
1

以为我应该发布我为此提出的修复程序。一旦我获得组合文件（上图），我使用 perl 哈希引用将它们读入内存并扫描文件 A。文件 A 中的匹配被哈希并用于扫描文件 B。这仍然需要大量内存，但运行速度非常快。从使用 grep 的 20 多天到约 20 分钟。

于 2014-06-07T19:13:23.117 回答