PheWAS(全表型组关联分析)----PheWAS实战(三)
2020-01-10 本文已影响0人
bcl_hx
介绍完PheWAS与GWAS的概念及其对比以及PheWAS的R包后，相信大家对对PheWAS原理以及相关R包有了初步认知。下面我以我做的一个项目的数据为例跑PheWAS。
1. 原始数据介绍
#基因型文件
第一行行名代表不同的个体
第一列rs代表SNP
B,H,D代表基因型(实际使用时候应该转化为0,1,2)
file
#表型文件
id列:10001.......代表不同的表型编号(你可以根据自己的表型编号)。
列名:BXD1.......代表不同的个体。
61.4,54.1.......代表表型值。
file
2.数据格式转化
我们要求转成的格式前面已经介绍过了，下面我直接放R中的代码。
#1.BXD_Geno文件处理
#1.1genotype BHD重编码为0，1，2
#加载读取exccel的包
library(readxl)
#读取BXD_Geno文件
BXD_Geno<-read_excel("/Users/bcl/Desktop/PheWAS2.0/genotype/BHD重编码/BXD_Geno.xlsx")
#把B替换为0，H替换为1，D替换为2
BXD_Geno[BXD_Geno=="B"]=0
BXD_Geno[BXD_Geno=="H"]=1
BXD_Geno[BXD_Geno=="D"]=2
#转置
BXD_Geno<-t(BXD_Geno)
df<-head(BXD_Geno,5)
#生成文件写入一个新的csv文件备用(第一个数去掉改为id，方便后面筛选用)
write.table(BXD_Geno,file="/Users/bcl/Desktop/PheWAS2.0/genotype/BHD重编码/BXD_Geno_recode.csv",quote=FALSE,sep=",",col.names = FALSE)
#1.2select_the_common_part(表型文件和基因型文件中都有的个体)
#加载包
library(dplyr)
#读取文件
df1<-read_excel("/Users/bcl/Desktop/PheWAS2.0/genotype/select_the_common_part/individual_with_phenotype.xlsx")
df2<-read.csv("/Users/bcl/Desktop/PheWAS2.0/genotype/BHD重编码/BXD_Geno_recode.csv")
#semi_join以df1为筛选标准，保留df2中与df1相同的部分。
df1<-tbl_df(df1)
df2<-tbl_df(df2)
select_BXD_Geno<-semi_join(df2,df1,by="id")
df<-head(select_BXD_Geno,5)
#输出为csv(选出的基因型数据，第一轮筛选基因型，根据的是pheno和geno共有的个体)
write.table(select_BXD_Geno,file="/Users/bcl/Desktop/PheWAS2.0/genotype/select_the_common_part/select_BXD_Geno.csv",quote=FALSE,sep=",",row.names=FALSE)
#2.BXD_phenotype文件处理
#读取文件
BXD_phenotype<-read_excel("/Users/bcl/Desktop/PheWAS2.0/phenotype/BXD_phenotype.xlsx",sheet=1)
select_BXD_Geno<-read.csv("/Users/bcl/Desktop/PheWAS2.0/genotype/select_the_common_part/select_BXD_Geno.csv",check.names = FALSE)
#提取BXD_phenotype中有价值的信息(ID列和各个个体表型值列)
BXD_phenotype_extract<-BXD_phenotype
#转置
BXD_phenotype_t<-t(BXD_phenotype_extract)
df<-head(BXD_phenotype_t,5)
#添加一列[id.....每个小鼠个体名称(BXD???)]
new_first_col<-rownames(BXD_phenotype_t)
BXD_phenotype_t<-cbind(new_first_col,BXD_phenotype_t)
df<-head(BXD_phenotype_t,5)
#列名命名为id，表型编号(100001......)，并删除第一行[用第一行数据去命名]
colnames(BXD_phenotype_t)<-BXD_phenotype_t[1,]
BXD_phenotype_t<-BXD_phenotype_t[-1,]
#行名命名为1：总行数
rownames(BXD_phenotype_t)<-c(1:nrow(BXD_phenotype_t))
df<-head(BXD_phenotype_t,5)
#以phenotype和genotype共有的个体为筛选标准，选出phenotype文件中符合的个体。
df3<-tbl_df(BXD_phenotype_t)
df4<-tbl_df(select_BXD_Geno[,1:2])
select_BXD_pheno<-semi_join(df3,df4,by="id")
df<-head(select_BXD_pheno,5)
#输出
write.table(select_BXD_pheno,file="/Users/bcl/Desktop/PheWAS2.0/phenotype/select_BXD_pheno.csv",quote=FALSE,sep=",",row.names=FALSE,col.names = TRUE)
#3.排序
#select_BXD_geno和select_BXD_pheno的文件按照其中一个文件的个体信息，对另外一个文件的个体信息进行排序(以select_BXD_geno第一列为标准)
select_BXD_Geno<-as.data.frame(select_BXD_Geno)
df<-head(select_BXD_Geno,5)
select_BXD_pheno<-as.data.frame(select_BXD_pheno)
df<-head(select_BXD_pheno,5)
for(i in 1:nrow(select_BXD_Geno)){ #genotype文件
for(j in 1:nrow(select_BXD_pheno))
{
if(select_BXD_Geno[i,1]==select_BXD_pheno[j,1]) #如果i和j的个体同
{
select_BXD_pheno[c(i,j),]=select_BXD_pheno[c(j,i),] #select_BXD_pheno的i行与j行互换
next; #不找了，下个循环
}
}
}
df<-head(select_BXD_pheno,5)
write.table(select_BXD_pheno,file="/Users/bcl/Desktop/PheWAS2.0/phenotype/reorder_select_BXD_pheno.csv",quote=FALSE,sep=",",row.names=FALSE)
#select_BXD_geno和reorder_select_BXD_pheno的第一列重新编号为1.....个体数(每个个体代表什么在re_select_geno文件中)
select_BXD_Geno$id<-c(1:nrow(select_BXD_Geno))
df<-head(select_BXD_Geno,5)
write.table(select_BXD_Geno,file="/Users/bcl/Desktop/PheWAS2.0/genotype/select_the_common_part/select_BXD_Geno_num.csv",quote=FALSE,sep=",",row.names=FALSE)
select_BXD_pheno$id<-c(1:nrow(select_BXD_pheno))
write.table(select_BXD_pheno,file="/Users/bcl/Desktop/PheWAS2.0/phenotype/reorder_select_BXD_pheno_num.csv",quote=FALSE,sep=",",row.names=FALSE)
df<-head(select_BXD_P,5)
4.运行PheWAS
#4.run the PheWAS
#加载包
library(PheWAS)
library(stringr)
#读取文件
re_select_BXD_Geno_num<-read.csv(file = "/Users/bcl/Desktop/PheWAS2.0/genotype/select_the_common_part/select_BXD_Geno_num.csv",header = TRUE,check.names = FALSE)
phenotypes<-read.csv(file = "/Users/bcl/Desktop/PheWAS2.0/phenotype/reorder_select_BXD_pheno_num.csv",header = TRUE,check.names = FALSE)
p_value_min_row<-data.frame(matrix(NA,ncol(re_select_BXD_Geno_num)*5,15))
k=1
for(i in 2:ncol(re_select_BXD_Geno_num)){
genotypes<-re_select_BXD_Geno_num[,c(1,i)] #选择Geno中的第1列和第i列(即每次选择一个snp)
results=phewas(phenotypes,genotypes,cores = 14) #得到p值结果
#index<-which(results$p==min(results$p),arr.ind = TRUE) #挑选p值最小的，然后返回索引。
#选出最小的5个，先升序然后输出最小5个所在的行
results1<-results[order(results[,7]),]
for(j in 1:5){
p_value_min_row[k,]<-results1[j,] #循环写入一个数据框
k=k+1
}
}
colnames(p_value_min_row)<-colnames(results) #赋列名
write.table(p_value_min_row,file="/Users/bcl/Desktop/PheWAS2.0/p_value_min_row.csv",quote=FALSE,sep=",",row.names=FALSE)
5.曼哈顿图绘制
#5.绘制曼哈顿图(按照gwas的方法做)[写文章挑选有代表性的结果绘图即可]
#要4列:第一列SNP，第二列CHR，第三列BP，第四列p(注:phewas结果中只有第二列CHR没有，要补一列)
library(qqman)
new_col_in_one<-genotypes[,2]
new_col_in_one[1:nrow(results)]<-colnames(new_col_in_one)#注:这里面的results为第4步phewas函数产生的一个结果
colnames(new_col_in_one)<-c("SNP")
#添加第二列CHR
new_col_CHR<-re_select_BXD_Geno[,2]
new_col_CHR[1:1663,1]<-1
colnames(new_col_CHR)<-c("CHR")
#设置图的属性
par(cex=0.8) #设置点的大小
color_set<-c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
genotypes<-re_select_BXD_Geno[,c(1,2)] #选择Geno中的第1列和第i列(即每次选择一个snp)
results=phewas(select_BXD_P,genotypes) #得到p值
#把前面的一列并进来，合为四列
results_for_manhattan<-cbind(results$snp,new_col_CHR,results$phenotype,results$p)
colnames(results_for_manhattan)<-c("SNP","CHR","BP","P")
results_for_manhattan$BP<-as.numeric(as.character(results_for_manhattan$BP)) #BP列转化为数值型
#run manhattan
manhattan(results_for_manhattan,main="PheWAS Manhattan Plot",xlim=c(10000,13000),genomewideline=2.0,col=color_set,xlab="") #调横坐标位置，因为我的表型编号值大概在10000-13000左右，导致做出来的图颜色并没有我预想的颜色，这里可以调整编号的方式实现上面设置的颜色。
file
本文由博客一文多发平台 OpenWrite 发布！
