用AWK和Perl处理大规模生物信息数据的实用技巧

问题描述

我需要将约450个文件中的生物信息学数据转换为矩阵格式以进行PCA分析。每个文件包含约5000行数据，总数据量约450万行。之前尝试使用R的pivot_wider()函数，但由于数据集过大而失败。

数据格式示例：

1
2
3


HG00097#1_chr1  D1000   15
HG00097#1_chr1  D1001   196
HG00097#1_chrX  D966856 1

每行包含三个字段：

唯一标识符（格式：[样本]#[单体型]_[染色体]）
D值（如D1000）
该值在特定染色体中出现的次数

目标输出

需要创建一个矩阵：

行名：所有样本标识符
列名：所有唯一的D值
单元格：对应样本和D值的出现次数（缺失值用0填充）

示例输出：

1
2
3
4


|                | D1000 | D1001 | D1002 | ... | D966856 |
|----------------|-------|-------|-------|-----|---------|
| HG00097#1_chr1 | 15    | 196   | 48    | ... | 0       |
| HG00097#1_chrX | 0     | 0     | 0     | ... | 1       |

解决方案

1. 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 file3 ...

特点：

使用哈希表存储数据
自动处理缺失值（用0填充）
输出为制表符分隔格式
内存消耗约4.6GB（处理450万行数据时）

2. AWK+sort解决方案

 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
27
28


#!/usr/bin/env bash

awk '
    BEGIN {
        OFS = "\t"
        while ( (getline < ARGV[1]) > 0 ) {
            haps[++numHaps] = $2
            vals[$2] = $2
        }
        ARGV[1] = ""
        sample = "Sample"
    }
    sample != $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 "$@") "$@"

AWK getline使用注意事项：

getline应谨慎使用，仅在特定场景下适用
安全的使用模式：if/while ( (getline var < file) > 0)
避免影响内置变量（$0, NF等）
使用后及时关闭文件：close(file)

3. Miller工具解决方案

1

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

特点：

专门设计用于数据处理
简洁的命令行语法
支持多种输出格式（CSV、TSV等）
处理450万行数据约需26秒

性能比较

工具	内存消耗	处理时间	适用场景
Perl	~4.6GB	中等	需要复杂数据处理
AWK	较低	快	内存受限环境
Miller	~4.6GB	26秒	简单数据转换

数据特性说明

文件组织：
- 单体型#1和#2的数据分别存储在不同文件中
- 每个文件内按D值从小到大排序
- 染色体按1-22、X、Y顺序排列
样本类型：
- XX样本：单体型#1和#2都以X结尾
- XY样本：单体型#1以Y结尾，单体型#2以X结尾
D值范围：
- 约330,303个唯一的D值
- 需要为每个缺失组合填充0值

实用建议

内存管理：处理大规模数据时，AWK方案通常内存效率更高
排序需求：如需要特定排序（如数值排序D值），可在Perl或AWK脚本中添加相应排序逻辑
输出格式：根据后续分析工具需求选择合适的输出格式（TSV、CSV等）
错误处理：确保处理过程中能正确处理文件读取错误和边界情况

这些方案均已在实际大规模生物信息学数据上测试验证，可根据具体计算环境选择最合适的工具。