ixxmu
commented
4 months ago 孟德尔随机化：下载UKbiobank数据本地处理 by 医学和生信笔记
今天学习下如何下载UKBiobank的数据并进行整理，把它变成TwoSampleMR能够使用的格式。
数据下载
打开ukbiobank的官网：https://www.nealelab.is/uk-biobank
打开后往下拉，找到results can be found here，点击它：
点击之后即可来到数据界面，首先是README这个sheet，这里会有一些介绍，每个文件包含哪些内容，每个文件的列名是什么意思等，都有详细的说明，大家一定要看！
为了方便以后使用，这个文件是可以免费下载的，建议大家直接下来，后面就可以用EXCEL随时打开查看了。
点击下面的Manifest 201807就可以来到数据详情和下载界面了，UKbiobank这个免费的数据已经好久不更新了，可以看到数据就留在2018年了。
数据详情界面有每个数据的phenotype description、phenotype code、数据下载链接等信息。
可以看到每个表型（phenotype）都有6个数据，分别是inrt和raw各有3种（分为female、male、both_sexes）。
把下面的滚动条往右边拉到最后，就可以看到每个数据的下载链接了，你鼠标放上去就会有复制链接的提示，直接复制在浏览器或者迅雷等软件中打开即可下载数据了。
每个表型都有6个数据，那么我们应该是用哪个呢？官方的github[1]介绍说选择irnt作为接下来的分析：
Overall, the results are largely consistent regardless of the choiceof IRNT(inverse rank-normal transformed) or raw untransformedphenotypes. Since IRNT does appear to provide a marginal boost to theh2g results, especially in terms of significance, we choose to treatthe IRNT version as the primary analysis for continuous phenotypes.(Results for the raw, untransformed versions will still appear in thecomplete results file though.)
所以我们也不纠结了，直接用irnt的数据分析，如果你的数据需要分开性别，那么你就单独下载某个性别的数据，要么就是下载both_sexes的。
数据读取及整理
我这里选择了fat这个表型的数据（上面有展示）。下载好之后直接用vroom读取即可：
fat_bothsex <- vroom::vroom("./ukbiobank/100004_irnt.gwas.imputed_v3.both_sexes.tsv.bgz")
dim(fat_bothsex)
## [1] 13791467 11
colnames(fat_bothsex)
## [1] "variant" "minor_allele" "minor_AF"
## [4] "low_confidence_variant" "n_complete_samples" "AC"
## [7] "ytx" "beta" "se"
## [10] "tstat" "pval"
这个数据也是很大的，1000多万行，11列。
# 看下前6行
head(fat_bothsex)
## # A tibble: 6 × 11
## variant minor_allele minor_AF low_confidence_variant n_complete_samples AC
## <chr> <chr> <dbl> <lgl> <dbl> <dbl>
## 1 1:1579… T 0 TRUE 51453 0
## 2 1:6948… A 2.04e-5 TRUE 51453 2.10e0
## 3 1:6956… C 1.80e-4 TRUE 51453 1.85e1
## 4 1:1398… T 2.03e-5 TRUE 51453 2.09e0
## 5 1:6927… C 1.12e-1 FALSE 51453 1.15e4
## 6 1:6937… G 1.17e-1 FALSE 51453 1.20e4
## # ℹ 5 more variables: ytx <dbl>, beta <dbl>, se <dbl>, tstat <dbl>, pval <dbl>
这个数据的每一列是什么意思，都在前面说过的README部分有详细的介绍，大家一定要去看！
比如第一列是variant，它的意思是chr:pos:ref:alt，alt是effect allele，需要使用这个和注释文件合并。这些信息都在README里面：
这个数据是没有rsid的，但是你仔细读过README之后就会发现，在variants.tsv.bgz这个文件中是有rsid的，而且这个文件也有chr:pos:ref:alt这些信息：
所以把这个文件下载下来，和我们的fat_bothsex文件合并一下，就有rsid了。variants.tsv.bgz这个文件就在Manifest 201807的第11行，找到下载链接直接下载即可：
下载好之后我们把它读取进来：
variants_anno <- vroom::vroom("./ukbiobank/variants.tsv.bgz")
dim(variants_anno)
## [1] 13791467 25
发现这个数据和fat_bothsex这个数据的行数是一样的。其实在README的描述中人家就说的很清楚了，这个顺序就是一样的，你可以通过variant这一列进行合并，也可以直接复制粘贴到一起：
查看前6行：
variants_anno[1:6,]
## # A tibble: 6 × 25
## variant chr pos ref alt rsid varid consequence consequence_category
## <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1:15791… 1 15791 C T rs54… 1:15… splice_reg… missense
## 2 1:69487… 1 69487 G A rs56… 1:69… missense_v… missense
## 3 1:69569… 1 69569 T C rs25… 1:69… missense_v… missense
## 4 1:13985… 1 139853 C T rs53… 1:13… splice_reg… missense
## 5 1:69279… 1 692794 CA C 1:69… 1:69… upstream_g… non_coding
## 6 1:69373… 1 693731 A G rs12… 1:69… upstream_g… non_coding
## # ℹ 16 more variables: info <dbl>, call_rate <dbl>, AC <dbl>, AF <dbl>,
## # minor_allele <chr>, minor_AF <dbl>, p_hwe <dbl>, n_called <dbl>,
## # n_not_called <dbl>, n_hom_ref <dbl>, n_het <dbl>, n_hom_var <dbl>,
## # n_non_ref <dbl>, r_heterozygosity <dbl>, r_het_hom_var <dbl>,
## # r_expected_het_frequency <dbl>
看看是不是完全一样：
# 肯定是完全一致
identical(variants_anno$variant, fat_bothsex$variant)
## [1] TRUE
所以我们直接合并一下即可：
fat_bothsex <- cbind.data.frame(variants_anno[,2:6], fat_bothsex)
head(fat_bothsex)
## chr pos ref alt rsid variant minor_allele minor_AF
## 1 1 15791 C T rs547522712 1:15791:C:T T 0.00000e+00
## 2 1 69487 G A rs568226429 1:69487:G:A A 2.03879e-05
## 3 1 69569 T C rs2531267 1:69569:T:C C 1.80062e-04
## 4 1 139853 C T rs533633326 1:139853:C:T T 2.03117e-05
## 5 1 692794 CA C 1:692794_CA_C 1:692794:CA:C C 1.12102e-01
## 6 1 693731 A G rs12238997 1:693731:A:G G 1.16788e-01
## low_confidence_variant n_complete_samples AC ytx beta
## 1 TRUE 51453 0.00000 0.00000 NaN
## 2 TRUE 51453 2.09804 -1.99658 -0.9567220
## 3 TRUE 51453 18.52940 -2.86950 -0.1222020
## 4 TRUE 51453 2.09020 -2.00107 -0.9583010
## 5 FALSE 51453 11535.90000 -162.96100 -0.0233942
## 6 FALSE 51453 12018.20000 -105.14000 -0.0126221
## se tstat pval
## 1 NaN NaN NaN
## 2 0.6939280 -1.378710 0.1679920
## 3 0.2382830 -0.512845 0.6080620
## 4 0.6939330 -1.380970 0.1672940
## 5 0.0107131 -2.183700 0.0289888
## 6 0.0101812 -1.239750 0.2150740
MR分析
首先是加载R包：
library(TwoSampleMR)
现在fat_bothsex这个数据是一个GWAS的原始数据，你可以把它用作是结局相关的数据，也可以把它用作是暴露相关的数据。
如果你要把这个数据用作是结局相关的数据，那就是直接进行以下代码即可：
# 注意列名一一对应
outcome_data <- format_data(dat = fat_bothsex,
type = "outcome", # 类型是outcome
snp_col = "rsid",
beta_col = "beta",
se_col = "se",
eaf_col = "minor_AF",
effect_allele_col = "alt",
other_allele_col = "ref",
pval_col = "pval",
samplesize_col = "n_complete_samples",
#ncase_col = "ncase_col",
#ncontrol_col = "ncontrol_col",
#gene_col = "nearest_genes",
chr_col = "chr",
pos_col = "pos"
)
现在这个outcome_data就是一个结局数据了，如果你已经有了暴露数据，那么就可以直接进行MR分析了。
如果你要把这个数据用作是暴露数据，那么首先需要根据P值筛选一下：
# 先去掉NaN
df1 <- fat_bothsex[!is.nan(fat_bothsex$pval),]
# 根据P值筛选,我这个数据没有小于5e-8的，所以选择了5e-6作为演示
df1 <- subset(df1, pval < 5e-6)
然后再格式化数据：
# 注意列名一一对应
exposure_data <- format_data(dat = df1,
type = "exposure", # 类型是exposure
snp_col = "rsid",
beta_col = "beta",
se_col = "se",
eaf_col = "minor_AF",
effect_allele_col = "alt",
other_allele_col = "ref",
pval_col = "pval",
samplesize_col = "n_complete_samples",
#ncase_col = "ncase_col",
#ncontrol_col = "ncontrol_col",
#gene_col = "nearest_genes",
chr_col = "chr",
pos_col = "pos"
)
然后再运行下clump（这个过程需要联网的，如果你想完全本地化运行，可参考下载IEU OPEN GWAS数据本地处理）：
# 这一步需要联网
exposure_data <- clump_data(exposure_data,
clump_kb=10000, clump_r2=0.001)
现在这个exposure_data就是一个暴露数据了，如果你已经有了结局数据，那么就可以直接进行MR分析了。
其中还有一些问题并没有说，比如一个ref有多个alt，或者minor_allele有多个等问题，以后再介绍，大家遇到了也可以自己搜索下。
搞定！easy！
参考资料
github: https://nealelab.github.io/UKBB_ldsc/select_topline.html
https://mp.weixin.qq.com/s/CCq9iIsti5wpCyItYAhmNQ