Skip to content

01 · 数据读取与标准化

对应模块R/adapters/read_dispatch.Rread_tsv.Rread_instrument_table.RR/02_standardize.R

作用

任意来源、任意格式的 GWAS 数据读进来,归一到一套内部标准 schema,供下游所有方法使用。这是"适配层 + 标准 schema"解耦的核心:换数据源不用动方法引擎。

内部标准 schema

SNP, chr, pos, effect_allele, other_allele, eaf, beta, se, pval, samplesize, ncase, ncontrol

对齐 GWAS-SSF 与 TwoSampleMR 习惯。读进来先转成 raw(标准列名 + 原始编码),再由 standardize_dataset() 做编码归一与清洗。

读取分发(read_gwas

统一入口 read_gwas(unit, cfg, role, snp_whitelist, window, pos_whitelist)format 分发到三个适配器:

format适配器典型来源
tsvread_tsv_rawFinnGen、GWAS-SSF、通用扁平 sumstats
instrument_tableread_instrument_table_rawUKB-PPP protein_info.csv 式预算工具表
vcfread_vcf_rawGWAS-VCF(M3 stub)

列名自动探测(detect_columns

不用手工写死列名。config.yamlcolumn_aliases 内置了各源的常见别名,探测时:

  1. column_map 显式指定了某列 → 用它(且校验该列真的在表头里,否则报错);
  2. 否则在别名字典里找匹配(大小写不敏感)。

例如效应等位的别名:[effect_allele, ea, EA, A1, alt, ALT, ALLELE1, hm_effect_allele] —— FinnGen 的 alt、UKB-PPP 的 ALLELE1、GWAS-SSF 的 effect_allele 都能自动认出。认不出会明确报错告诉你缺哪列。

内存安全:流式预过滤(read_tsv_filtered

结局 GWAS 常有上千万行,直接读会爆内存。read_tsv.R 的关键设计是读之前先用 shell 流式过滤,只把需要的行喂给 fread

  • snp_whitelist(MR 只需工具 SNP):grep -Fw -f 白名单 只捞出这些 rsID 的行;
  • window(chr,start,end)(coloc 只需区域):awk 只取该区间的行;
  • pos_whitelist(无 rsID 的大文件,按位置匹配):awkchr:pos 集合过滤。

峰值内存仅几百 MB。.gz 文件用 zcat 解压串流。

一个踩过的坑

rsID 预过滤时,表头与数据各用一条独立的 zcat 管道{ zcat file | head -1; zcat file | grep ... })。若共用一个管道,head 预读缓冲块后退出会触发 SIGPIPE,曾导致部分文件读到 0 行。

工具表专属处理(read_instrument_table

预算工具表一张表里多个蛋白、一行一个已选 pQTL,额外支持两种变换:

  • 行过滤 filter:如 { "cis.trans": "cis" } 只取 cis 行;
  • 复合 ID 拆列 variant_id_split:把 Variant.IDpattern: "chr:pos:OA:EA" 拆成四列。pattern 显式声明冒号各段含义(A0=OA, A1=EA),拆出的列优先于别名探测。

标准化变换(standardize_dataset

读进 raw 后统一处理:

1. p 值编码

pval_encoding: raw 原样;neglog1010^(-pval) 还原(FinnGen mlogp、UKB-PPP log10p 都是 -log10 尺度)。为防 p=0 破坏后续 log/qnorm,用 .Machine$double.xmin 兜底,p>1 截到 1。

2. 效应量归一

内部一律用 beta(= logOR)effect_type: OR | HR 或工具表标了 OR 列时,取 log(beta)beta 原样。

3. 清洗(并追踪剩余数)

剔除:缺失 SNP/beta/se/pvalse ≤ 0、非 ^[ACGT]+$ 的非法等位、EA == OA 的单态记录、重复 SNP。每步丢多少行都记进日志(STROBE-MR 追踪)。

4. 匹配键(position 模式)

无 rsID 的数据集(如 Mahajan T2D),用 chr:pos:排序等位 建键:等位取 pmin/pmax 排序后拼接,使键与效应等位方向无关(方向留给谐化对齐),比固定 NEA:EA 顺序更稳健;且小写以匹配 TwoSampleMR::format_datatolower

对应 config 字段

yaml
defaults:
  match_by: "rsid"           # rsid(默认,build-safe)| position
  require_rsid: true         # 强制两侧都有 rsID
column_aliases: {...}        # 列名别名字典(自动探测用)
exposures/outcomes:
  - format: "tsv"            # tsv | instrument_table | vcf
    pval_encoding: "raw"     # raw | neglog10
    effect_type: "beta"      # beta | OR | HR
    column_map: {...}        # 显式列名映射(留空=自动探测)
    filter: {...}            # 工具表行过滤
    variant_id_split: {...}  # 复合 ID 拆列

校验(01_config.R

  • build 一致性match_by=rsid 下要求暴露/结局同一基因组版本,否则报错(除非开 liftover);
  • rsID 占比validate_rsid() 要求 ≥50% 的 SNP 形如 rs[0-9]+,否则报错(position 模式跳过)。

下一步:选出的标准化暴露数据进入 工具变量选择

个人科研与运维文档 · 内容持续修订