【孟德尔随机化】下载Pan-Biobank 数据并作为SMR分析
【孟德尔随机化】下载Pan-Biobank 数据并作为SMR分析
大名鼎鼎的uk biobank【UK Biobank — Neale lab】不必再介绍，目前该数据库已经将跨种族的gwas数据也放到了主页——
目前的药靶分析不再局限于某一种族人群的数据，因跨种族的数据以其样本量巨大，为药靶分析锦上添花~
英国生物库收集了 50 万个基因和表型信息配对的个体，对常见疾病和性状的遗传病因学研究具有极大的价值。然而，该数据集的大多数全基因组分析仅使用欧洲血统的个体。分析一个更具包容性和多样性的数据集能提高分析能力和发现潜力。在这里，我们介绍了对 7,228 个表型进行的多血统分析，涉及 6 个大陆血统组，共计 16,131 项全基因组关联研究。我们在文章发表前向公众免费发布了这些汇总统计数据。
如何下载对应表型的数据？
点击谷歌链接以后，
就可以crt+F 查询感兴趣的表型，拖到表格后面，附有相应的下载链接：
欣喜地下载完数据，正准备分析，秃然发现没有rsid……怎么办呢？几千万行的数据用R包在本地转换显然是不现实的！
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.
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函数啦
# # 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数据也整理出相同格式的一列
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一下
##match
full_liver <-liver %>% mutate(.,rsid=newdat$rsid[match(.$varid,newdat$varid)])
head(full_liver)
选择SMR分析所需列并整理
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()
对数转换
range(Liver_gwas$p) #查看数值范围
Liver_gwas$p <- 10^(liver_gwas$p) ##对数转换
range(Liver_gwas$p) #再次查看数值范围
保存数据
View(head(Liver_gwas))
write.table(Liver_gwas,file = "gwas.txt", quote =FALSE,row.names = F)
到这里就结束了对pan-biobank一个表型gwas数据的整理，以此类推…… - END -