大名鼎鼎的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
里这一句提示~
代码语言:javascript复制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函数啦
代码语言: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数据也整理出相同格式的一列
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 -