解复用文件时如何解决打开文件限制?

How to work around open files limit when demuxing files?

我经常有大型文本文件(10-100GB 解压缩)根据每行中的条形码进行多路分解,实际上生成的单个文件(唯一条形码)的数量在 1K 到 20K 之间。为此,我一直在使用 awk,它完成了任务。但是,我注意到对较大文件(与使用的更多独特条形码相关)进行多路分解的速度要慢得多 (10-20X)。检查 ulimit -n 显示 4096 作为每个进程打开文件的限制,所以我怀疑速度变慢是由于 awk 的开销在解复用文件总数超过时被迫不断关闭和重新打开文件4096.

缺少 root 访问权限(即限制是固定的),可以使用哪些变通方法来规避此瓶颈?

我确实有每个文件中存在的所有条形码的列表,因此我考虑过分叉多个 awk 进程,其中每个进程都分配有一个相互排斥的条形码子集 (< 4096) 以供搜索。但是,我担心必须检查每一行的条形码以获得集合成员资格的开销可能会抵消不关闭文件的好处。

有没有更好的策略?

我没有嫁给 awk,所以欢迎使用其他脚本或编译语言的方法。


具体例子

数据生成(带条形码的 FASTQ)

以下生成的数据与我正在使用的数据类似。每个条目由 4 行组成,其中条形码是使用无歧义 DNA 字母表的 18 个字符的单词。

1024 个独特的条形码 | 100 万次阅读

cat /dev/urandom | tr -dc "ACGT" | fold -w 5 | \
awk '{ print "@batch."NR"_"[=10=]"AAAAAAAAAAAAA_ACGTAC length=1\nA\n+\nI" }' | \
head -n 4000000 > cells.1K.fastq

16384 个唯一条码 | 100 万次阅读

cat /dev/urandom | tr -dc "ACGT" | fold -w 7 | \
awk '{ print "@batch."NR"_"[=11=]"AAAAAAAAAAA_ACGTAC length=1\nA\n+\nI" }' | \
head -n 4000000 > cells.16K.fastq

awk 多路分解脚本

请注意,在这种情况下,我为每个唯一条形码写入 2 个文件。

demux.awk

#!/usr/bin/awk -f
BEGIN {
    if (length(outdir) == 0 || length(prefix) == 0) {
        print "Variables 'outdir' and 'prefix' must be defined!" > "/dev/stderr";
        exit 1;
    }
    print "[INFO] Initiating demuxing..." > "/dev/stderr";
}
{
    if (NR%4 == 1) {
        match(, /.*_([ACGT]{18})_([ACGTN]{6}).*/, bx);
        print bx[2] >> outdir"/"prefix"."bx[1]".umi";
    }
    print >> outdir"/"prefix"."bx[1]".fastq";

    if (NR%40000 == 0) {
        printf("[INFO] %d reads processed\n", NR/4) > "/dev/stderr";
    }
}
END {
    printf("[INFO] %d total reads processed\n", NR/4) > "/dev/stderr";
}

用法

awk -v outdir="/tmp/demux1K" -v prefix="batch" -f demux.awk cells.1K.fastq

或类似的 cells.16K.fastq.

假设您是唯一一个 运行 awk,您可以使用

验证打开文件的大致数量
lsof | grep "awk" | wc -l

观察到的行为

尽管文件大小相同,但具有 16K 个唯一条码的文件比只有 1K 个唯一条码的文件运行速度慢 10-20 倍。

没有看到任何示例 input/output 或您当前正在执行的脚本,这是非常猜测的,但如果您当前在字段 1 中有条形码并且正在做(假设 GNU awk,那么您没有自己的代码管理打开的文件):

awk '{print > }' file

那么如果管理打开的文件确实是您的问题,将其更改为:

sort file | '!=f{close(f};f=} {print > f}'

当然,以上是假设这些条形码值是什么,哪个字段保存它们,分隔字段的内容,输出顺序是否必须与原始顺序匹配,您的代码可能还做了什么随着输入的增加,速度变慢等等,因为您还没有向我们展示任何这些。

如果这不是您需要的全部内容,请编辑您的问题以包含缺少的 MCVE。


根据您的脚本更新问题以及输入是 4 行块的信息,我将通过在每条记录的前面添加键 "bx" 值并使用 NUL 来解决问题分隔 4 行块,然后使用 NUL 作为排序的记录分隔符和随后的 awk:

$ cat tst.sh
infile=""
outdir="${infile}_out"
prefix="foo"

mkdir -p "$outdir" || exit 1

awk -F'[_[:space:]]' -v OFS='\t' -v ORS= '
    NR%4 == 1 { print  OFS  OFS }
    { print [=12=] (NR%4 ? RS : "[=12=]") }
' "$infile" |
sort -z |
awk -v RS='[=12=]' -F'\t' -v outdir="$outdir" -v prefix="$prefix" '
BEGIN {
    if ( (outdir == "") || (prefix == "") ) {
        print "Variables 7outdir7 and 7prefix7 must be defined!" | "cat>&2"
        exit 1
    }
    print "[INFO] Initiating demuxing..." | "cat>&2"
    outBase = outdir "/" prefix "."
}
{
    bx1   = 
    bx2   = 
    fastq = 

    if ( bx1 != prevBx1 ) {
        close(umiOut)
        close(fastqOut)
        umiOut   = outBase bx1 ".umi"
        fastqOut = outBase bx1 ".fastq"
        prevBx1  = bx1
    }

    print bx2   > umiOut
    print fastq > fastqOut

    if (NR%10000 == 0) {
        printf "[INFO] %d reads processed\n", NR | "cat>&2"
    }
}
END {
    printf "[INFO] %d total reads processed\n", NR | "cat>&2"
}
'

当 运行 针对您在问题中描述的生成的输入文件时:

$ wc -l cells.*.fastq
4000000 cells.16K.fastq
4000000 cells.1K.fastq

结果是:

$ time ./tst.sh cells.1K.fastq 2>/dev/null

real    0m55.333s
user    0m56.750s
sys     0m1.277s

$ ls cells.1K.fastq_out | wc -l
2048

$ wc -l cells.1K.fastq_out/*.umi | tail -1
1000000 total

$ wc -l cells.1K.fastq_out/*.fastq | tail -1
4000000 total


$ time ./tst.sh cells.16K.fastq 2>/dev/null

real    1m6.815s
user    0m59.058s
sys     0m5.833s

$ ls cells.16K.fastq_out | wc -l
32768

$ wc -l cells.16K.fastq_out/*.umi | tail -1
1000000 total

$ wc -l cells.16K.fastq_out/*.fastq | tail -1
4000000 total