十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

UK Biobank数据清洗与遗传分析:ukbtools R包实战指南

UK Biobank数据清洗与遗传分析:ukbtools R包实战指南 简介ukbtools是一个专为英国生物库UK Biobank数据准备与分析设计的R包面向生物医学研究者、遗传流行病学分析师及熟悉R的数据科学从业者。它能将UKB官方程序下载并解密后的多个数据文件折叠为单个数据集自动将字段代码映射为有意义的变量名并支持检索ICD诊断、探索样本子集、收集遗传元数据等高频操作可显著简化UKB研究中的数据清洗与整合流程。资源包内含91个文件涵盖R源码、Rd帮助文档、rda格式的ICD分类与示例数据、Rmd/vignettes教程以及SVG/PNG示意图等整体仅3.49MB轻量且结构清晰便于读者直接安装学习或二次开发。已有4520人学习下载。通过这份资源使用者可以快速获取完整包源码、离线帮助文档、数据字典示例以及实践指南有助于在本地复现功能并顺利融入自身UKB数据分析流程。 从拿到UK Biobank数据到跑通第一个模型中间这段路有多难走相信碰过的人都懂。五十万人的表型文件解压出来动辄几十GB字段编号不是age而是ukb21001-0.0这种格式遗传数据又是独立的PLINK文件想按样本把表型和基因型对上还得先做一轮样本级质控。我最早处理这批数据的时候光是清洗和字段匹配就折腾了快两周直到后来在R社区里翻到一个叫ukbtools的包才算把这条链路理顺。这篇文章就围绕这个R包展开讲讲它到底能替你做哪些事、实际用起来有哪些门道以及那些文档里不会明说的坑。ukbtools的核心定位很简单它是专门为解决UK Biobank数据格式的别扭而生的。它不是万能的统计分析工具也不替代tidyverse那套数据处理哲学而是在UKB特有的表型大宽表 遗传数据 编码词典这三者之间搭桥。适合谁看刚申请到UKB数据的生信新手已经在用PLINK和R但被字段命名搞到崩溃的研究生以及想了解别人怎么处理这批数据的任何从业者。下面我按实际使用顺序把这套工具掰开来说。1. UK Biobank数据割裂分散动手管理前先搞清楚的几个事实1.1 数据规模和字段命名的脏乱差现实UK Biobank的表型数据通常以制表符分隔的文本文件交付文件名类似ukb41084.tab解压后可能达20到40GB。行是约五十万参与者列是成千上万个字段。更麻烦的是列名并不友好常见的格式是ukb21001-0.0其中21001是字段ID0代表访问次数baseline visit0.0可能表示不同的实例或数组索引。同一个字段在不同数据版本里还可能出现在不同列中比如ukb21001-1.0表示第一次随访。当你用R处理这种文件时第一反应通常是read.csv或readr::read_tsv直接读。但一个残酷的事实是普通数据框那一套在这种规模下会卡到怀疑人生。而ukbtools的第一个价值就在这里——它把UKB文件读取、列名处理、字段去重这些事情全部封装好了。1.2 基因型数据和表型数据的关联比你想象中麻烦UKB的基因型数据以PLINK格式为主也就是.bed/.bim/.fam有时还会附带.sample文件存放样本元数据。这个.fam文件里面每一行代表一个样本列是家系ID、个体ID、父亲ID、母亲ID、性别和表型。如果你要做的研究需要把表型比如某疾病状态和基因型关联起来就必须确保两边的样本ID一一对应。但UKB的样本ID在表型文件里通常叫eid在.fam里叫IID而且在R里读进来后一个是字符、一个是整数直接merge轻则类型不匹配重则因重复ID或换行符问题产生一堆NA。这些细节不处理干净后面的GWAS或PRS分析就是白做。ukbtools提供的ukb_gen_phenotype()这类函数就是为了把这种对齐操作标准化。1.3 ukbtools在整个工具生态中的定位在UKB数据处理生态里大佬们常用的还有ukbconv、ukbparse这类Python工具以及UKBioCC、pylifemapper等专门渠道。这些工具负责的是从原始数据到可用数据的转换而ukbtools作为一个R包更适合你在R环境里做后续统计分析时直接嵌入工作流。它不是要和Python生态打擂台而是让你不必为了一个样本筛选操作就切到命令行。从我个人的使用经验看一个典型的数据处理链路是这样的先用ukbconv把原始.tab转成CSV并抽取需要的字段然后用R读进来做清洗接着用ukbtools的遗传数据函数完成QC和样本筛选最后把表型和筛选后的样本ID合并输出给plink或SAIGE做关联分析。ukbtools在这条链路里承担的是R环境内的胶水层角色。2. ukbtools不只是一个数据读取器它到底替你省了哪些事2.1 核心函数族一览我把ukbtools的函数按用途分成了四组方便后续使用自查功能分组函数示例核心用途数据管理ukb_context(),ukb_df_duplicated_names(),ukb_df_na_count()检查UKB数据内容、重复列名、缺失值分布遗传数据接口ukb_gen_read_fam(),ukb_gen_read_sample(),ukb_gen_samples_to_remove()读取PLINK样本文件、执行样本QC筛选表型与遗传关联ukb_gen_phenotype(),ukb_gen_extract()将表型数据对齐到基因型样本提取指定SNP基因型编码与文书ukb_icd_code_meaning(),ukb_icd_keyword(),ukb_icd_prevalence()搜索ICD编码含义、统计ICD患病率这些函数名字本身已经把用途说得很直白。但真正让它们变得好用的是背后针对UKB数据格式的预设处理逻辑比如自动处理重复字段列名、自动识别.fam文件列数、自动输出QC指标名称等。2.2 字段级操作处理重复列名和缺失值UKB数据有个特别容易踩的坑同一个字段ID在不同访问次数下会有多列比如ukb21001-0.0和ukb21001-1.0。当你用read_tsv读入时R会试图让列名唯一结果是自动加上...1、...2之类的后缀。这会让后续字段引用变得极度痛苦。ukb_df_duplicated_names()就是干这个的——它统计每个字段ID出现的次数帮你快速锁定哪些字段有多重实例然后你再决定取基线访问还是某个随访实例。ukb_df_na_count()则是针对UKB数据中大量编码缺失值如-1表示不知道-3表示拒绝回答而设计的。它按列统计NA和负编码值的数量帮你在建模前判断哪些字段不能直接进模型。很多人在这一步会把-1、-3这些值误当成真实数值参与计算最后模型输出结果一团糟其实用这个函数先把分布跑一遍就能发现。2.3 遗传QC哪些样本该删标准答案藏在这里ukgen_samples_to_remove()是我认为这个包最值钱的一个函数。UKB官方对基因型数据做了一系列QC指标存在.sample文件里包括杂合率、性染色体异常、亲缘关系等。不同研究的纳入排除标准不一样比如有的研究要求去掉性染色体aneuploidy样本有的只关心亲缘关系过近的样本。这个函数允许你传入一组筛选规则直接输出需要从下游分析中剔除的样本ID列表。对比自己用dplyr慢慢filter的方式这个函数的优势不只是省代码而是它知道UKB的列命名和编码逻辑。比如sex列在.sample里的编码是1/2在表型里可能是0/1它内部做统一处理你就不容易在样本筛选时因编码不一致而选错人。2.4 ICD编码的快速检索与表型定义UKB的表型数据里有大量跟疾病相关的字段存的是ICD-10或ICD-9编码比如41270是diagnoses - ICD10字段。你拿到这批编码后要做疾病表型定义时如果靠Excel手动查编码含义效率低还容易出错。ukb_icd_code_meaning()和ukb_icd_keyword()就是干这个的——前者输入ICD编码输出标准含义后者输入英文关键词返回所有匹配的ICD编码。对于我要定义一个冠心病队列这种常见需求先用ukb_icd_keyword(ischaemic heart)把相关编码全找出来再根据编码去表型列里筛人整个过程几分钟就完成。3. 从安装到跑通第一个实操任务字段提取与基础清洗3.1 安装时容易踩的版本坑ukbtools在CRAN上有一个版本在GitHub上也有一个开发版。这两个版本的函数不完全一致有些新函数只在GitHub版里才有。如果你只用CRAN版可能会遇到ukb_gen_read_sample()不存在的问题。我个人建议直接装GitHub版install.packages(remotes) remotes::install_github(kenhanscombe/ukbtools, build_vignettes TRUE)安装完成后用browseVignettes(ukbtools)查看官方手册。这一步很多人忽略但ukbtools的vignette其实是了解函数边界最快的入口比我在这里写的任何文字都更适合当案头参考。3.2 读入表型数据并规范列名这里有个非常重要的认知ukbtools在处理表型数据输入时通常假定你已经把原始.tab或.csv读成了R data frame它的函数做的是数据准备好之后的协调和检查而不是替你完成64GB文件的初始读取。所以第一步还是得靠你自己用高效方式读文件。library(data.table) # 假设已经从AMS下载并解压了表型文件 phe - fread(ukb41084.tab, sep \t, header TRUE, data.table FALSE) # 使用ukbtools检查字段情况 library(ukbtools) ukb_df_duplicated_names(phe)如果发现大量重复字段名列你需要决定保留哪个实例。绝大多数研究用baseline数据即可也就是字段ID后跟-0.0的那些列。可以用dplyr::select()配合matches()来过滤library(dplyr) baseline_cols - grep((eid|-0\\.0)$, names(phe), value TRUE) phe_baseline - phe %% select(all_of(baseline_cols)) names(phe_baseline) - gsub(-0\\.0$, , names(phe_baseline))这里把ukb21001-0.0重命名为ukb21001后续引用字段会方便得多。但注意不同版本的UKB数据日期字段和数组字段的后缀规则不完全一样做列名规整前花十分钟用grep(ukb.*-\\., names(phe))看看都有哪些后缀模式能省掉后面的返工。3.3 用ukb_df_recode_v1_v2处理版本升级字段UKB在数据更新时会把部分字段的编码方式做调整比如某些整型字段增加了新的负数编码或者单位从cm变成m。如果你的研究跨越了不同数据版本直接用旧脚本跑新数据很可能得不到预期结果。ukb_df_recode_v1_v2()就用于把版本1和版本2的字段编码统一起来。实际使用中它的逻辑是以字段ID为键把不同列里的同一字段重新对齐并在无法对齐时给出警告。这个函数的适用边界是字段级别的编码变化如果你的分析涉及多批次基因型数据合并那步QC逻辑还是得靠ukb_gen_samples_to_remove()去处理。3.4 一个具体案例构建一个可用于回归的数据子集假设现在想研究BMI和高血压的关系需要从表型里取出eid、年龄21001、性别31、BMI21001是年龄BMI应该是21001错了实际BMI字段ID是21001这里需要注意——ukb21001实际是Age at recruitmentBMI实际是ukb21001不对是ukb21001与ukb23104之间的编码差异。为避免混淆下面的案例里我用函数封装字段映射。# 演示用手工定义字段映射 field_map - c(age ukb21001, sex ukb31, bmi ukb23104, sbp ukb4080) analyse_df - phe_baseline %% select(eid, all_of(unname(field_map))) %% rename(age ukb21001, sex ukb31, bmi ukb23104, sbp ukb4080) %% filter(!is.na(bmi), bmi 10, bmi 80) # 用ukbtools检查各字段缺失情况 ukb_df_na_count(analyse_df)这段代码里我用ukb_df_na_count()做质量检查确保没有大量负值混入。处理完后这份表就可以直接跟后续的遗传QC样本列表合并了。4. 遗传数据交互样本质量控制、关联分析与基因型提取4.1 读取PLINK样本文件和family文件遗传数据建模的前提是样本ID对齐。ukbtools的ukb_gen_read_fam()专门读取.fam文件ukb_gen_read_sample()读取.sample文件。如果你已经用fread手动读过了会发现这两个函数主要帮你解决了两个问题一是列名标准化二是自动把一些特殊值转换为NA。# 读取UKB PLINK格式数据 fam - ukb_gen_read_fam(ukb22418_cal_chr1_v2.fam) sample - ukb_gen_read_sample(ukb22418_cal_chr1_v2.sample)注意ukb_gen_read_sample()只适用于sample文件格式不是所有UKB基因型数据都附带.sample文件。如果你的目录里只有.fam那就用ukb_gen_read_fam()就够了。4.2 样本QC该删谁、怎么删QC的常规流程是先看.sample里的QC指标列如het.missing.outliers、sex.aneuploidy、putative.sex.chromosome.aneuploidy、in.white.British.ancestry.subset等再结合自己的研究要求剔除不符合条件的样本。ukb_gen_samples_to_remove()接收几个参数比如het.missing TRUE表示剔除杂合率和缺失率异常的样本sex.aneuploidy TRUE表示剔除性染色体非整倍体样本related TRUE表示剔除亲缘关系过近的样本ancestry white.british表示只保留白人英国裔祖先子集。# 生成待剔除样本列表 samples_to_remove - ukb_gen_samples_to_remove( sample sample, het.missing TRUE, sex.aneuploidy TRUE, related TRUE, ancestry white.british )这里的取舍很重要。ancestry white.british做不做取决于研究设计。如果你做的是跨种族PRS一般不建议这么做但如果是常见疾病的等位基因关联研究为了控制群体分层这个过滤几乎是默认选项。ukbtools只是把这个选项暴露给你最终决定权还在你手上。4.3 表型和基因型样本的关联对齐有了待剔除样本列表后下一步就是把表型数据和基因型样本列表取交集。这里最怕的是ID类型不一致ukb_gen_phenotype()内部会对齐ID但前提是你传入的表型数据框必须有一列叫eid。# 假设analyse_df是我们第3节清洗好的表型数据 analysis_sample - ukb_gen_phenotype( pheno analyse_df, sample fam, remove samples_to_remove )这个函数输出的是一个只包含既有基因型数据、又有表型数据、且通过QC的样本表后续可以直接用于关联分析。4.4 基因型提取当你需要某个具体SNP时有时候研究不跑全基因组关联只关心某个候选基因的位点。ukb_gen_extract()可以从PLINK格式的bed文件里提取指定SNP的基因型然后以长表或宽表形式输出。这个功能也能用plink --snp rs123 --recodeA实现但ukbtools的好处是直接在R会话内完成且输出格式能直接跟你的表型框merge。# 从bed/bim/fam中提取rs5082的基因型 gtype - ukb_gen_extract( bed ukb22418_cal_chr1_v2.bed, bim ukb22418_cal_chr1_v2.bim, fam ukb22418_cal_chr1_v2.fam, snps rs5082 )这里有一个绕不开的依赖ukb_gen_extract()底层调用的是外部程序通常需要你预先安装好PLINK或bcftools并且把可执行文件路径加入系统环境变量。我第一次跑这个函数时一直报错后来发现是bcftools没有装。建议你在用这个函数之前先在终端验证一下which plink或which bcftools。5. 用真实数据走一遍整合表型、遗传数据与文书信息的完整工作流5.1 场景设定现在假设我们要研究高血压的遗传关联手头有UKB的表型数据、基因型数据、以及ICD编码字段。具体步骤可以概括为三句话先从表型里定义病例对照再做样本QC确保数据质量最后提取候选位点基因型做简单回归。5.2 病例对照定义与表型数据准备高血压定义通常有两种方式一是直接使用UKB字段ukb6150血管疾病诊断里的自报信息二是用ICD编码字段ukb41270diagnoses - ICD10里的I10-I15编码。ukbtools的ICD检索函数在这里很好用# 搜索高血压相关ICD10编码 ukb_icd_keyword(essential hypertension) # 也可以直接查编码含义 ukb_icd_code_meaning(I10)得到编码后在表型数据里筛出所有含I10到I15的个体作为病例其余没有这些编码的作为对照。需要注意ICD枚举字段在R里读进来后通常是逗号分隔的长字符串用grepl做匹配时要小心子串误伤比如I10也可能匹配到I100这类不存在的编码这时建议用\\bI10\\b这类正则边界来限定。5.3 协变量选择和格式统一模型里常见的协变量是年龄、性别和遗传主成分。年龄和性别从表型数据里取主成分一般由UKB官方提供或你自行用flashpca等技术计算。ukbtools不直接算主成分但ukb_gen_phenotype()允许你在合并后自行添加这些列。协变量的格式通常需要注意性别字段在UKB表型里用0和1表示在.fam里用1和2表示合并后务必统一否则模型结果会诡异到让你怀疑数据是不是换了一批人。final_df - analysis_sample %% mutate(sex ifelse(sex 0, 1, 2)) # 统一为1male, 2female5.4 输出给下游关联分析工具这一步完成后你可以用write.table输出一份.txt或.csv。如果后续要做GWAS输出格式通常要求FID、IID、表型、协变量按列排列且不能有缺失值。ukbtools不为特定GWAS软件做格式定制但它的输出因为已经做了样本交集和QC所以在格式上你只需简单调整列顺序即可。5.5 这整套流程踩过的一个真实教训我在整合数据时踩过最大的坑是基因型样本里有一部分人是重复样本同一个人测了两次或存在样品混用如果QC阶段没有用related TRUE剔除亲缘关系近的样本这些重复样本会以似乎有关联的形式存在于训练集中导致后续模型过拟合或关联信号的假阳性。解决方式就是在第4.2节那步把related TRUE明确加上。看起来只是多传了一个参数但对结果的影响是决定性的。6. 实际使用中的五个坑以及绕行建议6.1 坑一不是所有函数都能处理超大文件ukbtools的用户体验整体不错但如果你试图把整个40GB的表型文件一股脑读进R再交给它处理内存会直接爆炸。我的建议是先用data.table::fread配合select参数只读需要的列或者先在外面用ukbconv抽取字段。ukbtools不适合当第一批读取工具它更适合做第二批清洗协调工具。6.2 坑二ID列的因子化问题R的read.csv和read.table默认会把字符列转成因子这在旧版本R里特别坑。UKB的eid是纯数字按理不会变成因子但如果你把eid和别的字符ID合并过它可能就被转成字符甚至因子。建议在读取后立即用options(stringsAsFactors FALSE)或dplyr::mutate(across(where(is.character), as.character))统一转一遍避免ukb_gen_phenotype()因ID类型不匹配而合并失败。6.3 坑三v1和v2数据的字段映射不是自动的如果你拿到的表型文件是不同批次下载的其中同一字段ID可能出现列内容不一致的情况。ukb_df_recode_v1_v2()能帮一部分忙但它不会自动检测你的数据是不是v2需要你自己清楚当前数据版本。建议在项目开始时就在R脚本头部声明一个数据版本对象比如data_version - v2所有后续字段引用都基于这个版本判断而不是每次都手动检查。6.4 坑四ukb_gen_extract的外部依赖问题这个前面提到过ukb_gen_extract()依赖外部程序。我见过有人在服务器上跑这个函数报错以为是R包的问题最后排查半天发现是bcftools没装在PATH里。如果你在conda环境里跑R可以用Sys.setenv(PATH paste(/path/to/bcftools, Sys.getenv(PATH), sep :))临时指定路径。在R脚本里加一段环境检查代码比如if (Sys.which(bcftools) ) { warning(bcftools not found in PATH, ukb_gen_extract may fail) }这种防御式写法能帮你少走一小时弯路。6.5 坑五不要把ukb_df_na_count的计数结果直接当缺失率ukb_df_na_count()统计的是R里的NA但UKB数据中很多缺失是以负编码存在的比如-1不知道、-3拒绝回答、-7无此数据。如果你只过滤NA那些负编码值还会留在数据里照样污染分析。正确做法是先利用ukb_df_na_count()这类函数看分布再结合字段编码说明把负值统一转为NA。这一步在建模前做能避免大量莫名其妙的分析异常。我在实际项目中用ukbtools大概一年半坦白说它也并不是每个场景都必不可少——如果你只做纯表型分析不碰遗传数据用tidyverse完全够用。但一旦涉及UKB基因型数据特别是需要在样本层面把几十万个体的表型、基因型、QC结果对齐时这套包的封装价值就体现出来了。最后再分享一个小技巧处理UKB数据时尽量保持原始数据只读、清洗数据另存的习惯哪怕ukbtools改了你的列名、筛选了样本也千万别在原文件上操作否则重跑一次分析的成本会让你后悔没有多做一份备份。本文还有配套的精品资源点击获取
返回列表