使用AWK和Perl进行生物信息学数据重构以支持PCA分析

一位用户在处理用于主成分分析（PCA）的数据时遇到了挑战。他们需要将大约450个文件中的数据整理成一个矩阵格式。每个文件包含类似以下结构的数据行（例如，样本 HG00097，单倍型 1，染色体 1 和 X）：

1
2
3
4
5
6
7
8

HG00097#1_chr1  D1000   15
HG00097#1_chr1  D1001   196
HG00097#1_chr1  D1002   48
HG00097#1_chr1  D1003   55
HG00097#1_chr1  D1012   38
HG00097#1_chr1  D1018   2
HG00097#1_chrX  D966856 1
HG00097#1_chrX  D974    1

其中，第一列是样本标识符，第二列是“D”值，第三列是该值在特定染色体中出现的次数。

目标是将所有文件合并到一个矩阵中，其中：

• 所有个体样本标识符作为行名。
• 所有唯一的“D”值作为列标题。
• 单元格包含对应的出现次数，缺失值用0填充。

基于提供的示例，期望的输出矩阵（以制表符或逗号分隔）应如下所示（此处以表格形式展示以便阅读）：

用户最初尝试使用R的pivot_wider()函数，但由于数据集巨大（约450万行）而失败，因此寻求基于命令行工具（如AWK）的解决方案。

社区提供了多种解决方案：

1. Perl 解决方案 一种简洁的Perl单行命令可以生成制表符分隔的矩阵：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

perl -lane '
  $f{$F[0]}->{$F[1]} = $F[2];
  $d{$F[1]} = 0;
  END {
    $, = "\t";
    @d = sort keys %d;
    print "", @d;
    for $f (sort keys %f) {
      print $f, map {$f{$f}->{$_} // 0} @d;
    }
  }' -- file1 file2 ...

这个脚本读取所有文件，构建一个两层嵌套的哈希结构来存储计数，并收集所有唯一的“D”值。在END块中，它先输出列标题（所有排序后的“D”值），然后遍历每个样本标识符，输出其对应的计数，缺失值用0补全。输出可以通过管道传递给mlr（Miller）等工具进行美观打印。

2. Miller (mlr) 解决方案 Miller是专门为数据操作设计的工具，可以更直接地处理这种“重塑”操作：

1

cat file* | mlr --n2c reshape -s 2,3 then unsparsify --fill-with 0 then label id

该命令执行以下操作：

• --n2c: 将输入视为没有表头的、以空格分隔的字段，并输出CSV格式。
• reshape -s 2,3: 将第2个字段（“D”值）的值作为新的列名，并将第3个字段的值放入对应的列下。
• unsparsify --fill-with 0: 确保所有记录（行）都有相同的字段（列），缺失的列用0填充。
• label id: 将第一个字段的标签命名为“id”。

回答者指出，在测试中处理包含9000个唯一“D”值的模拟数据（约450万行）时，该命令运行时间约为26秒，内存消耗约4.6 GB。

3. AWK 与 sort 组合的解决方案 一个纯AWK（结合sort）的解决方案旨在更节省内存。它分两步处理：首先使用sort获取所有唯一的“D”值作为列顺序，然后使用AWK逐文件处理数据，只将当前文件的计数和全局的列顺序保留在内存中。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

awk '
    BEGIN {
        OFS = "\t"
        while ( (getline < ARGV[1]) > 0 ) {
            haps[++numHaps] = $2
            vals[$2] = $2
        }
        ARGV[1] = ""
        sample = "Sample"
    }
    FNR == 1 {
        prt()
        sample = $1
    }
    { vals[$2] = $3 }
    END { prt() }

    function prt(       hapNr) {
        printf "%s", sample
        for ( hapNr=1; hapNr<=numHaps; hapNr++ ) {
            printf "%s%s", OFS, vals[haps[hapNr]]
            vals[haps[hapNr]] = 0
        }
        print ""
    }
' <(sort -u -k2,2 "$@") "$@"

这个脚本假设每个输入文件对应一个样本标识符（FNR == 1时切换）。prt函数负责打印每一行。

4. 算法思路讨论 有回答者从算法角度分析了问题本质：这是一个大规模数据透视（pivot）问题。关键在于高效地确定矩阵的行数（唯一样本标识符数量）和列数（唯一“D”值数量），然后以紧凑格式（如二维数组）存储和填充数据。虽然数据集行数多，但如果唯一值数量可控，现代计算机的内存和计算能力足以应对。建议可以使用C++、Python或R（通过流式或分块处理）来实现。

讨论与性能对比 在后续讨论中，用户澄清总行数4.5M是指合并后作为R函数输入的数据量，每个文件大约有5000行。用户测试发现Perl版本速度稍慢。回答者比较了Perl、AWK和Miller的速度，指出在特定测试数据集（约50万行）上，Perl（不含mlr后处理）比AWK版本更快（25秒 vs 90秒），而mlr本身在处理大数据集时可能因内存不足而崩溃。具体的性能表现高度依赖于唯一“D”值的数量和样本数量。

此外，用户补充了关于数据排序和性染色体（X/Y）分布的额外信息，这些信息可能影响最终矩阵的行顺序和完整性，但在核心的数据重构逻辑上，上述解决方案仍然适用。