读取公共数据库中GWAS结果的vcf格式文件

这里以https://gwas.mrcieu.ac.uk/数据库为例
这个数据库中的vcf文件是gwas分析之后的，所以可直接生成GWAS数据文件而不需表型数据。index文件是vcf文件的索引。

library(vcfR)arData <- vcfR::read.vcfR('/home/zhang123/cs/ukb-b-.vcf.gz') #读取vcf文件str(arData) #查看vcf数据的结构head(arData@meta,12)#查看注释信息head(arData@fix)head(arData@gt)# 数据处理arFix <- as.data.frame(arData@fix[,(1:5)]) #提取arData@fix的前5列arGt <- as.data.frame(arData@gt[,2]) #提取arData@gt的第一列colnames(arGt) <- "GT" #对arGt的列名重命名beta <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(1,dim(arGt)[1]*5,5)]) #提取每个SNP的beta值se <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(2,dim(arGt)[1]*5,5)]) #提取每个SNP的se值eaf <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(4,dim(arGt)[1]*5,5)]) #提取每个SNP的eaf值pval <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(3,dim(arGt)[1]*5,5)]) #提取每个SNP的-log10后的P值mydata <- data.frame(beta = beta, se = se, pval = pval, maf = eaf ) #合并数据mydata <- cbind(arFix,mydata)mydata$pval <- 10^(-mydata$pval) #转化为原始P值head(mydata)mydata$N <- 462933write.table(mydata, "/home/zhang123/cs/ukb-b-.txt", quote = FALSE, row.names = FALSE)

原作文章链接：https://mp.weixin.qq.com/s/_CvOR5Lckjafzvsm7xYAxw
以上是花了3rmb买完整理后的。同时也希望大家支持原创，毕竟连奶茶都不到。

收到优化代码评论，按需存取：

number <- length(unlist(strsplit(as.character(arData@gt[1,1]),split=":")))#根据对应的注释内容，计数一个numberbeta <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(1,dim(arGt)[1]*number,number)]) #提取每个SNP的beta值个人觉得这样改一下更好，适应更多场景