grep .fastq 中的每四行

Question

我正在 linux 机器上使用 bash。

我的问题是，如何使用 grep 跳过查询文件中的行？

我正在处理一个名为 example.fastq 的 ~16Gb .fastq 大文件，它具有以下格式。

example.fastq

@SRR6750041.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR6750041.2 2/1
CTATANTATTCTATATTTATTCTAGATAAAAGCATTCTATATTTAGCATATGTCTAGCAAAAAAAA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
@SRR6750041.3 3/1
ATCCANAATGATGTGTTGCTCTGGAGGTACAGAGATAACGTCAGCTGGAATAGTTTCCCCTCACAG
+
AAAAA#EE6E6EEEEEE6EEEEAEEEEEEEEEEE//EAEEEEEAAEAEEEAE/EAEEA6/EEA<E/
@SRR6750041.4 4/1
ACACCNAATGCTCTGGCCTCTCAAGCACGTGGATTATGCCAGAGAGGCCAGAGCATTCTTCGTACA
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAE/E/<//AEA/EA//E//
@SRR6750041.5 5/1
CAGCANTTCTCGCTCACCAACTCCAAAGCAAAAGAAGAAGAAAAAGAAGAAAGATAGAGTACGCAG
+
AAAAA#EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEAEEEAEEE<EE/E

我需要提取包含感兴趣的字符串的行 @SRR6750041.2 @SRR6750041.5 存储在名为 IDarray 的 bash 数组中，以及每个匹配项之后的 3 行。下面的 grep 命令允许我这样做

for ID in "${IDarray[@]}";
    do
    grep -F -A 3 "$ID " example.fastq 
    done

这正确输出了以下内容。

@SRR6750041.2 2/1
CTATANTATTCTATATTTATTCTAGATAAAAGCATTCTATATTTAGCATATGTCTAGCAAAAAAAA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
@SRR6750041.5 5/1
CAGCANTTCTCGCTCACCAACTCCAAAGCAAAAGAAGAAGAAAAAGAAGAAAGATAGAGTACGCAG
+
AAAAA#EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEAEEEAEEE<EE/E

我正在寻找加快此过程的方法...一种方法是通过将搜索限制为以 @ 开头的行或跳过不可能包含匹配项的行来减少 grep 搜索的行数@SRR6750041.1 例如第 2、3、4 和 6、7、8 行等。有没有办法使用 grep 来做到这一点？也欢迎使用其他方法！

Answer 1

这里有一些想法和例子。出于测试目的，我创建了测试用例，因为你的 example_mini.fastq 的迷你版本有 145 MB 大，IDarray 有 999 个元素（兴趣）。

您的版本具有此性能（用户 space 超过 2 分钟）：

$ time for i in "${arr[@]}"; do grep -A 3 "${i}" example_mini.fastq; done 1> out.txt
real    3m16.310s
user    2m9.645s
sys     0m53.092s

$ md5sum out.txt
8f199a78465f561fff3cbe98ab792262  out.txt

首次匹配后将 grep 升级为结束 grep -m 1，我假设兴趣 ID 是唯一的。这缩小了 50% 的复杂性，并在用户 space:

中花费了大约 1 分钟

$ time for i in "${arr[@]}"; do grep -m 1 -A 3 "${i}" example_mini.fastq; done 1> out.txt
real    1m19.325s
user    0m55.844s
sys     0m21.260s

$ md5sum out.txt
8f199a78465f561fff3cbe98ab792262  out.txt

这些解决方案与元素数量线性相关。对大文件调用 n 次 grep。

现在让我们只在 AWK 中实现一个运行，我正在将 IDarray 导出到输入文件中，这样我就可以在一个运行中进行处理。我正在将大文件加载到每个 ID 的关联数组中，然后通过你的 ID 数组循环 1x 进行搜索。这是通用场景，您可以在其中定义正则表达式和打印后的行数。这具有复杂性，只有一个运行通过文件 + N 比较。这是 2000% 的加速：

$ for i in "${arr[@]}"; do echo $i; done > IDarray.txt
$ time awk '
(FNR==NR) && (linesafter-- > 0) { arr[interest]=arr[interest] RS [=12=]; next; }
(FNR==NR) && /^@/ { interest=; arr[interest]=[=12=]; linesafter=3; next; }
(FNR!=NR) && arr[] { print(arr[]); }
' example_mini.fastq IDarray.txt 1> out.txt
real    0m7.044s
user    0m6.628s
sys     0m0.307s

$ md5sum out.txt
8f199a78465f561fff3cbe98ab792262  out.txt

如您的标题所示，如果您真的可以确认每四行是感兴趣的 id，之后的三行将被打印。你可以简化成这个并再加速 20%:

$ for i in "${arr[@]}"; do echo $i; done > IDarray.txt
$ time awk '
(FNR==NR) && (FNR%4==1) { interest=; arr[interest]=[=13=]; next; }
(FNR==NR) { arr[interest]=arr[interest] RS [=13=]; next; }
(FNR!=NR) && arr[] { print(arr[]); }
' example_mini.fastq IDarray.txt 1> out.txt
real    0m5.944s
user    0m5.593s
sys     0m0.242s

$ md5sum out.txt
8f199a78465f561fff3cbe98ab792262  out.txt

对于包含 999 个元素的 1.5 GB 文件，搜索时间为：

real    1m4.333s
user    0m59.491s
sys     0m3.460s

所以根据我在我的机器上的预测，你的 15 GB 示例包含 10k 个元素，用户 space 需要大约 16 分钟来处理。

grep .fastq 中的每四行

grep every fourth line in .fastq

grep

fastq