使用 Hail 输出 PLINK 文件：一步到位

在现代基因组学研究中，数据的处理和转换是至关重要的一环。PLINK 文件（包括 .bed, .bim, .fam 三个文件）是遗传学研究中的标准文件格式之一，经常用于 GWAS（全基因组关联研究）等分析。对于研究人员来说，将原始数据转换为 PLINK 文件格式是分析过程中的关键步骤。而 Hail，作为一个处理大规模基因数据的高效工具，提供了方便的方法来导出这些文件。

前面，我们介绍了如何掌握 Hail，这里，我们将介绍如何使用 Hail 的 export_plink 功能，轻松导出 PLINK 文件，并深入探讨如何在数据预处理和分析过程中使用它。

使用 Hail 导出 PLINK 文件

Hail 提供了一个非常便利的函数 export_plink()，能够将 Hail MatrixTable 格式的数据导出为 PLINK 所需的三种文件格式：.bed、.bim 和 .fam。这些文件是 PLINK 软件包用于处理和分析基因型数据的标准文件格式。

1. 载入数据

首先，我们需要加载基因型数据，通常这些数据存储在 VCF 格式的文件中。使用 Hail 的 import_vcf() 函数，可以轻松导入 VCF 文件。

import hail as hl

mt = hl.import_vcf('file://path_to_vcf_data.vcf', force_bgz=True, reference_genome='GRCh38')

2. 数据注释

为了确保导出的 PLINK 文件包含完整的信息，我们通常需要对数据进行注释。例如，标注样本的性别、家族信息等。这可以通过 annotate_cols() 函数完成。

mt = mt.annotate_cols(
    pat_id='0',         # 父代 ID
    mat_id='0',         # 母代 ID
    is_female=True,     # 是否为女性
    pheno=-9            # 表型（-9 表示缺失值）
)

在这个例子中，我们为每个样本添加了父母 ID、性别和表型信息。

3. 导出为 PLINK 文件

现在，数据已经准备好，可以使用 export_plink() 函数将其导出为 PLINK 所需的 .bed, .bim, .fam 文件。

hl.export_plink(
    mt, 'output/example', 
    fam_id=mt.s,       # 样本的 ID
    ind_id=mt.s,       # 样本的 ID
    pat_id=mt.pat_id,  # 父代 ID
    mat_id=mt.mat_id,  # 母代 ID
    is_female=mt.is_female, # 性别信息
    pheno=mt.pheno     # 表型信息
)

参数解析： - fam_id: 样本的家族 ID，通常使用样本的 ID（即 mt.s）。 - ind_id: 样本的个体 ID，这里使用 mt.s 作为样本 ID。 - pat_id 和 mat_id: 分别代表父母的 ID，通常可以使用默认值 0 表示没有父母信息，或者根据数据 - is_female: 样本的性别信息，Hail 会根据性别表达为 1（男性）或 2（女性）。这里通过注释 is_female 字段来指示性别。 - pheno: 表型数据，可以是布尔值（如是否患病）或者是数值型数据（如身高、体重等）。

导出的文件会包含以下三部分：

• output/example.bed: 存储基因型数据（二进制格式）。
• output/example.bim: 包含变异信息，如染色体、位置、参考和变异等。
• output/example.fam: 包含样本信息，如家族 ID、个体 ID、父母 ID、性别和表型。

进阶功能：批量导出与错误处理

当处理大量基因数据时，可能会涉及多个数据批次。在这种情况下，使用批处理来导出 PLINK 文件变得尤为重要。例如，我们可以使用 Python 脚本批量处理多个数据集，并通过适当的异常处理机制确保过程顺利进行。

import subprocess

# 批量导出并上传
batch_files = ['batch1.vcf', 'batch2.vcf', 'batch3.vcf']  # 示例批次
for batch_file in batch_files:
    # 导入 VCF 文件并处理
    mt = hl.import_vcf(f'file://{batch_file}', force_bgz=True, reference_genome='GRCh38')
    mt = mt.annotate_cols(
        pat_id='0',
        mat_id='0',
        is_female=True,
        pheno=-9
    )
    
    # 导出 PLINK 文件
    hl.export_plink(mt, f'output/{batch_file.split(".")[0]}', fam_id=mt.s, ind_id=mt.s, pat_id=mt.pat_id, mat_id=mt.mat_id, is_female=mt.is_female, pheno=mt.pheno)

    # 上传到远程存储
    subprocess.run(['dx', 'upload', f'output/{batch_file.split(".")[0]}.bed'])
    subprocess.run(['dx', 'upload', f'output/{batch_file.split(".")[0]}.bim'])
    subprocess.run(['dx', 'upload', f'output/{batch_file.split(".")[0]}.fam'])