从多样本文本文件提取特定STR位点的高效方案咨询

背景与需求 #

现有122个STR分析结果文件，每个文件包含约80万条独特STR位点数据，单文件示例如下：

SAMPLE  CHROM   POS Allele_1    Allele_2    LENGTH
HG02035 chr1 230769616 (tcta)14 (tcta)16 4
HG02035 chr2 1489653 (aatg)8 (aatg)11 4
HG02035 chr2 68011947 (tcta)11 (tcta)11 4
HG02035 chr2 218014855 (ggaa)16 (ggaa)16 4
HG02035 chr3 45540739 (tcta)15 (tcta)16 43

另一文件示例：

SAMPLE  CHROM   POS Allele_1    Allele_2    LENGTH
HG02040 chr1 230769616 (tcta)15 (tcta)15 4
HG02040 chr2 1489653 (aatg)8 (aatg)8 4
HG02040 chr2 68011947 (tcta)10 (tcta)10 4
HG02040 chr2 218014855 (ggaa)21 (ggaa)21 4
HG02040 chr3 45540739 (tcta)17 (tcta)17 4

需要将所有样本中**同一STR位点（按CHROM+POS唯一标识）**的记录汇总到单个文件，输出示例（以chr1:230769616位点为例）：

HG02035 chr1 230769616 (tcta)14 (tcta)16 4
HG02040 chr1 230769616 (tcta)15 (tcta)15 4
HG02121 chr1 230769616 (tcta)2 (tcta)2 4
HG02131 chr1 230769616 (tcta)16 (tcta)16 4
HG02513 chr1 230769616 (tcta)14 (tcta)14 4

现有方案的问题 #

尝试使用以下awk命令实现需求：

awk '$1!="SAMPLE" {print $0 > $2"_"$3".locus.tsv"}' *.vcf

该命令逻辑正确，但因需要创建80万个文件，且每条记录都会触发一次文件打开/关闭操作，IO开销极大，导致处理耗时过长。

优化解决方案 #

方案：先排序再批量处理 #

核心思路是通过排序让同一位点的记录连续出现，再用awk仅在位点切换时才打开/关闭对应文件，将每个文件的IO操作从N次（N为该位点的样本数）降到1次，大幅减少IO开销。

执行命令：

# 合并所有文件、过滤表头、按位点排序后批量写入
cat *.vcf | grep -v "^SAMPLE" | sort -k2,2 -k3,3n | awk '{
    curr_key = $2 "_" $3
    if (curr_key != prev_key) {
        # 关闭上一个位点的文件（首次循环跳过）
        if (prev_key != "") close(prev_file)
        prev_file = curr_key ".locus.tsv"
        prev_key = curr_key
    }
    # 追加当前行到对应文件
    print >> prev_file
}'

命令解释 #

1. cat *.vcf：合并所有输入的STR结果文件
2. grep -v "^SAMPLE"：过滤掉所有文件的表头行（保留有效数据行）
3. sort -k2,2 -k3,3n：按染色体（第2列）和位点位置（第3列，数字排序）对数据行排序，确保同一位点的所有样本记录连续排列
4. awk处理逻辑：
   • 用curr_key存储当前位点的标识（CHROM_POS）
   • 当位点切换时，关闭上一个位点的文件，打开新位点的文件
   • 将当前行追加到对应位点的文件中，每个文件仅打开/关闭一次

额外优化建议 #

• 如果目标磁盘IO性能较差，可先将临时输出写入SSD分区，完成后再移动到目标目录
• 若服务器有多核CPU，可结合parallel对文件进行分组预处理，但上述单命令方案已能解决核心性能问题

内容的提问来源于stack exchange，提问作者Bách Nguyễn