hail

Hail 是一个开源、分布式、专门面向大规模基因组学数据的计算框架，由 Broad Institute 主导开发，用 Python + Spark 实现。它把基因组学中常见的 VCF、PLINK、BGEN、gVCF 等格式抽象成 MatrixTable / VariantDataset（VDS）这样的高维稀疏矩阵，允许用户在 TB 级数据上像操作 NumPy/Pandas 一样写几行代码就完成 GWAS、罕见变异分析、质量控制、药物靶点发现等任务。

国际上部分大的人群队列项目采用了 Hail 作为分析框架，如 All of Us Research Program(百万人)、UK Biobank (50万人)、gnomAD、Trans-biobank analysis with 676,000 individuals elucidates the association of polygenic risk scores of complex traits with human lifespan(67万)、Hematopoietic mosaic chromosomal alterations increase the risk for diverse types of infection(76万)

项目地址 https://github.com/hail-is/hail

官方文档 Hail 0.2

学习资源

https://github.com/UK-Biobank/UKB-RAP-Notebooks-Genomics/blob/main/Notebooks/G208_GWAS_in_Hail.ipynb

https://documentation.dnanexus.com/science/using-hail-to-analyze-genomic-data

Hail-GWAS教程笔记

格式介绍¶

MatrixTable（MT）¶

MatrixTable 是 Hail 的核心数据结构，把 变异矩阵 抽象为一个 行 × 列 × 条目 的三维张量：

行键：(locus, alleles) —— 每个变异位点
列：样本
条目：基因型、深度、概率等字段

物理存储

目录结构：xxx.mt/ 下包含 .parquet 格式的列式存储和 metadata.json。
支持列/行随机采样、快速过滤、分布式线性代数运算。

典型用途

QC、GWAS、PCA、群体遗传学分析、机器学习。

Variant Dataset（VDS）¶

VDS 是 Hail 专为 gVCF → 联合调用设计的双层结构：

reference_data：MatrixTable，把大片参考序列压缩成区间行（REF block）。
variant_data：MatrixTable，仅保存真正的变异位点。

物理存储

目录结构：xxx.vds/ 下有两个子目录

xxx.vds/
  ├── reference_data.mt/
  └── variant_data.mt/

VDS 磁盘占用比 VCF/MT 小 5–20×。

典型用途

大规模 gVCF 合并、增量追加样本、低频变异检测。

vcf 和 MT/VDS 比较¶

VCF 适合交换与归档，Hail 的 MT/VDS 才是高效、可扩展的分析引擎。

维度	VCF	MatrixTable / VDS
磁盘	文本 + 重复 REF 行 → 巨大	列式压缩 + 区间合并 → 小 5–20×
速度	行式解析，随机访问慢	列式 + 索引，Spark 并行，秒级过滤
增量合并	必须重跑全部样本	VDS 支持 `combine_variant_datasets` 增量合并
内存	全文件解压	仅加载需要的列/行
功能	基本查看/过滤	内置 PCA、线性回归、机器学习、群体统计
容错	单节点	分布式 Spark 容错
格式多样性	易出 header/编码问题	Hail 统一类型系统，自动类型检查

gvcf 合并¶

主要使用函数 hl.vds.new_combiner()，其选项说明如下

hl.vds.new_combiner(
    output_path: str,                         # 最终 VDS 输出目录
    temp_path: str,                           # 临时目录（生命周期需自行清理）
    save_path: str | None = None,             # 断点续跑计划文件 *.json
    gvcf_paths: list[str] | None = None,      # 输入 gVCF 路径列表
    vds_paths: list[str] | None = None,       # 输入 VDS 路径列表（用于二次合并）
    intervals: list[hl.Interval] | None = None,  # 手动区间列表
    import_interval_size: int | None = None,  # 按 bp 大小自动切区间
    use_genome_default_intervals: bool = False,  # 使用内置 WGS 分区
    use_exome_default_intervals: bool = False,   # 使用内置 WES 分区
    reference_genome: hl.ReferenceGenome = 'default',
    contig_recoding: dict[str, str] | None = None,  # 染色体名映射
    gvcf_external_header: str | None = None,        # 外部 header 文件
    gvcf_sample_names: list[str] | None = None,     # 强制指定样本名
    gvcf_info_to_keep: list[str] | None = None,     # 要保留的 INFO 字段
    gvcf_reference_entry_fields_to_keep: list[str] | None = None,
    call_fields: list[str] = ['PGT'],
    branch_factor: int = 100,           # 合并树分支系数
    target_records: int = 24_000,       # 每个分区目标记录数
    gvcf_batch_size: int = 50,          # 每批同时处理的 gVCF 数量
    force: bool = False,                # 覆盖已存在输出
)

从 gVCF 构建全基因组 VDS

hail_joint_gvcf.py

import hail as hl
import glob

hl.init()

irgsp_ref = hl.ReferenceGenome.from_fasta_file(
    name='IRGSP_1_0',
    fasta_file="IRGSP-1.0_genome.fasta",
    index_file="IRGSP-1.0_genome.fasta.fai",
)

intervals_tbl = hl.import_locus_intervals('IRGSP-1.0_genome.1based.bed', reference_genome='IRGSP_1_0')

intervals = intervals_tbl.interval.collect()

gvcfs = sorted(
    glob.glob('/data/gvcf/*.g.vcf.gz')
)

combiner = hl.vds.new_combiner(
    output_path='/data/output/dataset.vds',
    temp_path='/data/output/tmp/',
    save_path='/data/output/combiner_plan.json',  # 可断点续跑
    gvcf_paths=gvcfs,
    reference_genome=irgsp_ref,
    intervals=intervals,
    gvcf_batch_size=20,
    branch_factor=50,
)

combiner.run()

运行 python hail_joint_gvcf.py，如果样本数较多可能出现 out of memory 的报错，需要提高 driver-memory 的大小（默认只有1-2G），即 PYSPARK_SUBMIT_ARGS="--driver-memory 20g pyspark-shell" python hail_joint_gvcf.py。

中断后只需 hl.vds.new_combiner(save_path='/data/output/combiner_plan.json').run() 即可继续。

在 vds 中追加 gvcf

hail_append_gvcf.md

import hail as hl
hl.init()

# 1) 已有的参考基因组
irgsp_ref = hl.ReferenceGenome.from_fasta_file(
    name='IRGSP_1_0',
    fasta_file="IRGSP-1.0_genome.fasta",
    index_file="IRGSP-1.0_genome.fasta.fai",
)

# 2) 读取旧 VDS
old_vds = hl.vds.read_vds('/data/output/dataset.vds')

# 3) 新增加的 gVCF 列表
new_gvcfs = [
    '/data/gvcf2/sample1.g.vcf.gz',
    '/data/gvcf2/sample2.g.vcf.gz',
    '/data/gvcf2/sample3.g.vcf.gz',
    '/data/gvcf2/sample4.g.vcf.gz',
]

# 4) 把新的 gVCF 组合成一个新的 VDS（临时目录）
tmp_vds = '/data/output/new_samples.vds'

intervals_tbl = hl.import_locus_intervals('IRGSP-1.0_genome.1based.bed', reference_genome='IRGSP_1_0')

intervals = intervals_tbl.interval.collect()


combiner = hl.vds.new_combiner(
    output_path=tmp_vds,
    temp_path='/data/output/tmp/',
    gvcf_paths=new_gvcfs,
    reference_genome=irgsp_ref,
    intervals=intervals,
)
combiner.run()

new_vds = hl.vds.read_vds(tmp_vds)

# 5) 按样本维度合并（行维度自动对齐）
merged_vds = hl.vds.combiner.combine_variant_datasets(
    [old_vds, new_vds]
)

# 6) 写回新的总 VDS（或覆盖旧目录）
merged_vds.write('/data/output/dataset2.vds', overwrite=True)

import hail as hl
hl.init()

irgsp_ref = hl.ReferenceGenome.from_fasta_file(
    name='IRGSP_1_0',
    fasta_file="IRGSP-1.0_genome.fasta",
    index_file="IRGSP-1.0_genome.fasta.fai",
)

intervals_tbl = hl.import_locus_intervals('IRGSP-1.0_genome.1based.bed', reference_genome='IRGSP_1_0')

intervals = intervals_tbl.interval.collect()

vds_list = sorted(
    glob.glob('/data/vds/*.dataset.vds')
)

combiner = hl.vds.new_combiner(
    output_path=tmp_vds,
    temp_path='/data/output/tmp/',
    vds_paths=vds_list,
    reference_genome=irgsp_ref,
    intervals=intervals,
)
combiner.run()

本文阅读量次
本站总访问量次