用于在输入文件的每一行上执行 sed 表达式的慢 bash 脚本

Question

我有一个简单的bash脚本如下

#!/bin/bash
#This script reads a file of row identifiers separated by new lines 
# and outputs all query FASTA sequences whose headers contain that identifier.
# usage filter_fasta_on_ids.sh fasta_to_filter.fa < seq_ids.txt; > filtered.fa


while read SEQID; do
    sed -n -e "/$SEQID/,/>/ p"  | head -n -1
done

fasta 文件具有以下格式：

> HeadER23217;count=1342
ACTGTGCCCCGTGTAA
CGTTTGTCCACATACC
>ANotherName;count=3221
GGGTACAGACCTACAC
CAACTAGGGGACCAAT

编辑更改了 header 名称以更好地显示其在文件中的实际结构

我上面制作的脚本确实正确地过滤了文件，但是速度很慢。我的输入文件有大约 20,000,000 行包含大约 4,000,000 个序列，并且我有一个包含 80,000 header 的列表，我想对其进行过滤。使用 bash/sed 或其他工具（如 python 或 perl）是否有更快的方法来完成此操作？为什么上面的脚本需要数小时才能完成？

Answer 1

您正在扫描大文件 8 万次。我将建议使用不同工具的不同方法：awk。将选择列表加载到哈希图（awk 数组）中，并在扫描大文件时是否有任何序列匹配打印。

例如

$ awk -F"\n" -v RS=">" 'NR==FNR{for(i=1;i<=NF;i++) a["Sequence ID " $i]; next}
                         in a' headers fasta

-F"\n" 标志将输入文件中的 field separator 设置为新行。 -v RS=">" 将记录分隔符设置为 ">"

Sequence ID 1
ACTGTGCCCCGTGTAA
CGTTTGTCCACATACC

Sequence ID 4
GGGTACAGACCTACAT
CAACTAGGGGACCAAT

headers 文件包含

$ cat headers
1
4

并且 fasta 文件包含更多相同格式的记录。

如果您的 header 已经包含 "Sequence ID" 前缀，请按此调整代码。我没有针对大文件对此进行测试，但只要您没有内存限制来容纳 80K 大小的数组，它应该比您的代码快得多。在那种情况下，将 header 拆分为多个部分并合并结果应该是微不足道的。

要允许 header 的任何格式并使生成的文件成为有效的 FASTA 文件，您可以使用以下命令：

awk -F"\n" -v RS=">" -v ORS=">" -v OFS="\n" 'NR==FNR{for(i=1;i<=NF;i++) a[$i]; next}  in a' headers fasta > out

ORS 和 OFS 标志设置输出字段和记录分隔符，在这种情况下与输入 fasta 文件相同。

Answer 2

你应该利用这个事实（你没有明确说明，但我假设）巨大的 fasta 文件包含顺序（按 ID 排序）的序列。

我还假设 headers 文件是按 ID 排序的。如果不是，那就这样做 - 对 80k 整数进行排序并不昂贵。

当对两个文件进行排序时，它归结为对两个文件进行一次同时线性扫描。由于它在常量内存中运行，因此它可以使用任何大小，这与其他 awk 示例不同。我在 python 中给出了一个示例，因为我对 awk 中的手动迭代不满意。

import sys

fneedles = open(sys.argv[1])
fhaystack = open(sys.argv[2])

def get_next_id():
    while True:
        line = next(fhaystack)
        if line.startswith(">Sequence ID "):
            return int(line[len(">Sequence ID "):])

def get_next_needle():
    return int(next(fneedles))

try:
    i = get_next_id()
    j = get_next_needle()
    while True:
        if i == j:
            print(i)
        while i <= j:
            i = get_next_id()
        while i > j:
            j = get_next_needle()
except StopIteration:
    pass

当然它有点冗长，但它在我的旧机器上大约 10 秒内找到 80k 的 4M 序列（339M 的输入）。（它也可以用 awk 重写，这可能会快得多）。我这样创建了 fasta 文件：

for i in range(4000000):
    print(">Sequence ID {}".format(i))
    print("ACTGTGCCCCGTGTAA")
    print("ACTGTGCCCCGTGTAA")
    print("ACTGTGCCCCGTGTAA")
    print("ACTGTGCCCCGTGTAA")

并且 headers ("needles") 这样：

import random

ids = list(range(4000000))
random.shuffle(ids)
ids = ids[:80000]
ids.sort()
for i in ids:
    print(i)

Answer 3

它很慢，因为您正在多次读取同一个文件，而您本可以 sed 读取一次并处理所有模式。所以你需要生成一个 sed 脚本，其中包含每个 ID 的语句和 />/b 来替换你的 head -n -1.

while read ID; do
    printf '/%s/,/>/ { />/b; p }\n' $ID;
done | sed -n -f - data.fa

用于在输入文件的每一行上执行 sed 表达式的慢 bash 脚本

Slow bash script to execute sed expression on each line of an input file

bash

parsing

sed

fasta