【孟德尔随机化】下载Pan-Biobank 数据并作为SMR分析

2023-11-30 12:26:22 浏览数 (2)

大名鼎鼎的uk biobank【UK Biobank — Neale lab】不必再介绍,目前该数据库已经将跨种族的gwas数据也放到了主页——

目前的药靶分析不再局限于某一种族人群的数据,因跨种族的数据以其样本量巨大,为药靶分析锦上添花~

英国生物库收集了 50 万个基因和表型信息配对的个体,对常见疾病和性状的遗传病因学研究具有极大的价值。然而,该数据集的大多数全基因组分析仅使用欧洲血统的个体。分析一个更具包容性和多样性的数据集能提高分析能力和发现潜力。在这里,我们介绍了对 7,228 个表型进行的多血统分析,涉及 6 个大陆血统组,共计 16,131 项全基因组关联研究。我们在文章发表前向公众免费发布了这些汇总统计数据。

如何下载对应表型的数据?

点击谷歌链接以后,

就可以crt F 查询感兴趣的表型,拖到表格后面,附有相应的下载链接:

欣喜地下载完数据,正准备分析,秃然发现没有rsid……怎么办呢?几千万行的数据用R包在本地转换显然是不现实的!

代码语言:javascript复制
lidat<- read_delim(file = "../Primary_outcome/categorical-20002-both_sexes-1158.tsv.bgz",delim="t")
view(head(lidat))

还有一些缺失值也是需要过滤的。再看看,发现p值竟然有负值!

这个时候大家需要回头看看README manifest里这一句提示~

Please note: the p-values in the summary statistic flat files are negative log10 p-values. P-value columns have been renamed with a "neglog10" suffix to indicate this.

代码语言:javascript复制
gwas$p <- 10^(gwas$p) ##对数转换
注意这里不能用exp函数进行简单的转换,因为exp并不以10为指数

结局了p值的问题,细心的小伙伴已经发现了端倪,之前在下载页面我打了两个箭头,再认真看看:

  • The variant manifest contains detailed information on each variant (download on Amazon AWS, tbi).

其实rsid的信息已经统一放在这里了,点击链接即可下载。

然后呢?就是很简单的match函数啦

代码语言:javascript复制
# # rsid文件
# snpdat <-read_delim(file = "../Primary_outcome/full_variant_qc_metrics.txt.bgz",delim = "t") 
# head(snpdat)
# newdat <- snpdat[,1:6]
# head(newdat)
# saveRDS(newdat,file="snpdat.rds")

match之前要看看newdat的varid列的样子,需要将gwas数据也整理出相同格式的一列

代码语言:javascript复制
liver <- lidat %>% as.data.frame() %>% drop_na() %>% 
          unite(.,"tmp",c("chr","pos"),sep=":",remove=T) %>% 
           unite(.,"varid",c("tmp","ref","alt"),sep="_",remove=F) 
head(liver)

match一下

代码语言:javascript复制
##match
full_liver <-liver %>% mutate(.,rsid=newdat$rsid[match(.$varid,newdat$varid)])
head(full_liver)

选择SMR分析所需列并整理

代码语言:javascript复制
liver_gwas <- full_liver[,c("rsid","alt","ref","af_cases_EUR","beta_EUR","se_EUR","pval_EUR")]
liver_gwas$n <- "1000"
colnames(liver_gwas) <- c("SNP","A1","A2","freq","b","se","p","n")
liver_gwas[1:3,]
Liver_gwas <-liver_gwas %>%  drop_na()

对数转换

代码语言:javascript复制
range(Liver_gwas$p) #查看数值范围
Liver_gwas$p <- 10^(liver_gwas$p) ##对数转换
range(Liver_gwas$p) #再次查看数值范围

保存数据

代码语言:javascript复制
View(head(Liver_gwas))
write.table(Liver_gwas,file = "gwas.txt", quote =FALSE,row.names = F)

到这里就结束了对pan-biobank一个表型gwas数据的整理,以此类推…… - END -

0 人点赞