【UK Biobank干货】教你如何在RAP云计算平台创建表型 ...
Research Analysis Platform（简称RAP，地址：https://ukbiobank.dnanexus.com/landing）是UKB与AMS和DNAnexus合作推出的云计算平台，功能强大，是未来数据分析和计算的趋势。本次推送介绍如何在RAP系统中创建表型数据（Phenotype）并进行Pre-GWAS质量控制，为后续的GWAS分析做准备。需使用到RAP平台的JupyterLab功能（Spark Cluster）和相应的Python代码。
【注意1】在依照本教程进行表型数据提取之前，应已确定好目标表型并进行了充分准备（例如在UKB官网查询数据的基本信息）。本教程以糖尿病为例。
【注意2】RAP使用前提：UK Biobank账号申请成功；Proposal申请成功；拥有AMS账号；数据使用费支付成功；数据使用协议签署成功。
【注意3】RAP有数据存储和使用费，但UKB直接release的raw data存储不收取费用，具体价格见（2022年8月价格）：
【注意4】RAP注册和使用初级教程见：
1、登录RAP系统，进入目标Project，点击上部菜单中的TOOLS，选择JupyterLab。
2、新建JupyterLab Environment：点击右上角“+New JupyterLab”图标。
3、设置JupyterLab Environment（Spark Cluster）：
【注意】2个地方需要进行设置：
（1）在"Project"项目中选择你自己的目标Project；
（2）在“Cluster Configuration”项目选择“Spark Cluster”【如果选的是默认的“Single Note”，后续运行代码时会报错，所以这里务必要注意】。
其他选项一般按默认值就可以。
【注意】页面底部提示了云计算收费标准（单位：英镑）。例如这里JupyterLab每4小时的使用费是4.032英镑。干完活儿务必记得End Environment，否则会一直计费。
点击"Start Environment"后，RAP需要花费一些时间进行初始化，耐心等待即可。
初始化完成 ，Status显示绿色的Ready，点击Name栏中对应的名称即可进入 JupyterLab Environment（新开窗口）。
4、设置Python环境：点击Notebook中的Python 3：
进入Python环境：
可以在左栏中将Untitled.ipynb文件重命名为demo1.ipynb：
5、【数据获取】我们首先导入dxdata、dxpy等packages用于获取和管理RAP上存储的UKB数据集。
【注意】我们从这里开始使用 Python代码，输入相应代码后记得运行。
# 导入packagesimport dxpy
import dxdata
import pandas as pd
import pyspark
6、自动发现分配的数据集ID并加载数据集：
# 自动发现分配的数据集ID并加载数据集dispensed_dataset_id = dxpy.find_one_data_object(typename='Dataset', name='app*.dataset', folder='/', name_mode='glob')['id']
dataset = dxdata.load_dataset(id=dispensed_dataset_id)
7、Spark初始化：
# Spark初始化（只运行一次，除非重启Kernel）sc = pyspark.SparkContext()
spark = pyspark.sql.SparkSession(sc)
8、访问participant entity（参与人信息）：
# 访问participant entityparticipant = dataset['participant']
9、根据Field Index、Instance Index、Array Index等信息选择参与人数据：
# 这个函数用来抓取一个字段ID列表的所有字段名（e.g. "pdef field_names_for_ids(field_ids):
from distutils.version import LooseVersion
fields = []
for field_id in field_ids:
field_id = 'p' + str(field_id)
fields += participant.find_fields(name_regex=r'^{}(_i\d+)?(_a\d+)?$'.format(field_id))
return sorted([field.name for field in fields], key=lambda n: LooseVersion(n))
field_ids = ['31', '21022', '22001', '22006', '22009', '22019', '22021', '22027',
'41202', '41204',
'20107', '2946', '1807',
'20110', '3526', '1845']
field_names = ['eid'] + field_names_for_ids(field_ids)
10、连接到Spark并将其转换为Pandas DataFrame
# 连接到Spark，并将其转换为Pandas DataFramedf = participant.retrieve_fields(names=field_names, engine=dxdata.connect())
pdf = df.toPandas()
11、【从这里开始，我们可以在Pandas DataFrame中进行数据交互】可视化表格：
# 显示表格的前6行pdf.head(6)
12、提取目标表型（以糖尿病为例）的相关UKB数据编码：
假设我们感兴趣的表型是糖尿病，根据UKB主页中的ICD-10诊断数据（UKB Data-Field 41202，见 https://biobank.ndph.ox.ac.uk/ukb/field.cgi?id=41202），得知其子编码为E10-E14 Diabetes mellitus，包括E10、E11、E13和E14：
# 利用关键词查找fields的function
def field_codes_by_keyword(field_name, keyword):
return dict([(k,v) for (k,v) in participant[field_name].coding.codes.items() if keyword.lower() in v.lower()])
# 确定糖尿病相关ICD-code（E10*、E11*、E13*、E14*）
diabetes_icd_codes = list({**field_codes_by_keyword("p41202", "e10."), **field_codes_by_keyword("p41202", "e11."), **field_codes_by_keyword("p41202", "e13."), **field_codes_by_keyword("p41202", "e14.")})
diabetes_icd_codes
运行后得到的diabetes_icd_codes即为糖尿病相关子编码。
13、构建表型变量（在本例中：has_diabetes_icd10变量）：
【定义标准】如果一个人的疾病记录中包含糖尿病相关编码，则定义表型变量has_diabetes_icd10为1，否则为0。
# 根据住院信息获取疾病情况def icd10_codes(row):
main_icd10_codes = row['p41202'] or []
secondary_icd10_codes = row['p41204'] or []
return list( set(main_icd10_codes+secondary_icd10_codes) )
pdf['icd10_codes'] = pdf.apply(icd10_codes, axis=1)
# 如果是糖尿病相关ICD-10子编码，记录为“1”
def has_diabetes_icd10(row):
return 0 if set(row['icd10_codes']).isdisjoint(diabetes_icd_codes) else 1
pdf['has_diabetes_icd10'] = pdf.apply(has_diabetes_icd10, axis=1)
# 显示新表格的前6行
pdf.head(6)
拉到最右，可以看到has_diabetes_icd10被成功定义：
14、质量控制（QC）：
【QC标准（可根据研究需要灵活修改）】
保留所报告性别与遗传性别相同的参与者；
保留英国白人族裔；
保留没有性染色体非整倍体的参与者；
过滤掉那些异质性和缺失率的离群者，或区分出过度的亲属关系，例如有10个或更多的三级亲属。
# 质量控制QCpdf_qced = pdf[
(pdf['p31'] == pdf['p22001']) & # Filter in sex and genetic sex are the same
(pdf['p22006']==1) & # in_white_british_ancestry_subset
(pdf['p22019'].isnull()) & # Not Sex chromosome aneuploidy
(pdf['p22021']!=10) & # Not Ten or more third-degree relatives identified (not 'excess_relatives')
(pdf['p22027'].isnull()) # Not het_missing_outliers
]
# 查看通过QC的样本量
pdf_qced.shape[0]
这里显示共有406,947个样本通过了QC。
15、可视化表型变量分布情况：
# 查看我们得到的表型变量分布情况pdf_qced['has_diabetes_icd10'].hist(bins=30)
糖尿病作为二元表型变量，病例/对照（Case/Control）的值分别为1/0。
16、整理得到的Pandas DataFrame并保存到Project中，以便用于后续GWAS分析：
# 整理数据# 重命名变量以便于理解
import re
pdf_qced = pdf_qced.rename(columns=lambda x: re.sub('p22009_a','pc',x))
pdf_qced = pdf_qced.rename(columns={'eid':'IID', 'p31': 'sex', 'p21022': 'age',
'p22006': 'ethnic_group',
'p22019': 'sex_chromosome_aneuploidy',
'p22021': 'kinship_to_other_participants',
'p22027': 'outliers_for_heterozygosity_or_missing'})
# 增加FID列
pdf_qced['FID'] = pdf_qced['IID']
# 将ICD10 code由list变为string
pdf_qced['icd10_codes'] = pdf_qced['icd10_codes'].apply(lambda x: 'NA' if x is None else ",".join([str(i) for i in x]))
# 生成QC后的表型数据table
pdf_phenotype = pdf_qced[['FID', 'IID', 'sex', 'age',
'has_diabetes_icd10',
'icd10_codes', 'has_diabetes_icd10',
'ethnic_group', 'sex_chromosome_aneuploidy',
'kinship_to_other_participants',
'outliers_for_heterozygosity_or_missing'] + [f'pc{i}' for i in range(1, 41)]]
# 显示新表格的前6行
pdf_phenotype.head(6)
17、关联基因数据（用时较长，耐心等待）：
# 关联Microarray基因数据path_to_family_file = '/mnt/project/Bulk/Genotype Results/Genotype calls/ukb22418_c1_b0_v2.fam'
plink_fam_df = pd.read_csv(path_to_family_file, delimiter='\s', dtype='object',
names = ['FID','IID','Father ID','Mother ID', 'sex', 'Pheno'], engine='python')
# 显示新表格的前6行
plink_fam_df.head(6)
18、将表型变量与基因变量关联：
# 将表型变量与基因变量关联并生成新的表型DateFramediabetes_microarray_df = pdf_phenotype.join(plink_fam_df.set_index('IID'), on='IID', rsuffix='_fam', how='inner')
# 删除*.fam中不用的列
diabetes_microarray_df.drop(
columns=['FID_fam','Father ID','Mother ID', 'sex_fam', 'Pheno'], axis=1, inplace=True, errors='ignore'
)
19、将表型数据保存为*.csv格式文件：
# 将表型数据保存为*.csv格式文件diabetes_microarray_df.to_csv('diabetes_microarray_df.phe', sep='\t', na_rep='NA', index=False, quoting=3)
pdf_phenotype.to_csv('diabetes.phe', sep='\t', na_rep='NA', index=False, quoting=3)
成功后左栏中会显示新生成的csv文件。
20、将我们生成的表型csv文件上传到RAP中的新建文件夹中（这里新建demo_diabetes文件夹作为演示）：
%%bash# 将我们生成的表型csv文件上传到RAP中
dx upload diabetes_microarray_df.phe -p --path /demo_diabetes/ --brief
dx upload diabetes.phe -p --path /demo_diabetes/ --brief
上传成功会显示相应的task number（本例中是file-GK5P90QJg95JjV7p6YPqy3vz和file-GK5P918Jg954Z7Ff6B2zgvQG）。
21、至此，我们成功提取并构建了经过质量控制的糖尿病表型数据。
返回RAP界面，进入相应的Project，可以看到刚才新建的demo_diabetes文件夹，里面存放了diabetes_microarray_df.phe和diabetes.phe表型文件：
点击文件可浏览概况：
表型文件下载到电脑后可以用R等统计软件读取。
【再次提醒】数据提取成功千万别忘了关JupyterLab，不然会一直扣费 ！！
更多UKB攻略/干货，请继续关注「遗传社科研究」后续推送！