读取公共数据库中GWAS结果的vcf格式文件
这个数据库中的vcf文件是gwas分析之后的,所以可直接生成GWAS数据文件而不需表型数据。index文件是vcf文件的索引。原作文章链接:https://mp.weixin.qq.com/s/_CvOR5Lckjafzvsm7xYAxw。以上是花了3rmb买完整理后的。同时也希望大家支持原创,毕竟连奶茶都不到。这里以https://gwas.mrcieu.ac.uk/数据库为例。
·
这里以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 <- 462933
write.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=":")))#根据对应的注释内容,计数一个number
beta <- as.numeric(unlist(strsplit(as.character(arGt$GT),split=":"))[seq(1,dim(arGt)[1]*number,number)]) #提取每个SNP的beta值
个人觉得这样改一下更好,适应更多场景
更多推荐
已为社区贡献1条内容
所有评论(0)