GWAS cFDR多效基因实例分析

这篇文章主要介绍"GWAS cFDR多效基因实例分析"的相关知识，小编通过实际案例向大家展示操作过程，操作方法简单快捷，实用性强，希望这篇"GWAS cFDR多效基因实例分析"文章能帮助大家解决问题。

GWAS cFDR 分析

library(dplyr)library(ggplot2)library(RColorBrewer)#library(GWAScFDR) #library(condFDR)library(cfdr)mycol=c(brewer.pal(7,"Set2"),rev(brewer.pal(8,"Set1")),brewer.pal(8,"Set3"),brewer.pal(12,"Paired"),brewer.pal(8,"Dark2"),brewer.pal(7,"Accent"))stratifiedQQplot = function(p1,p2,                            xlab=_expression(nominal-log[10]~~(q[expected])),                            ylab=_expression(empirical-log[10]~~(p[observed])),t="Stratified Q-Q Plot"){        library(ggplot2)        #p1 = p1[p1!=0]        #p2 = p2[p2!=0]                dat = NULL        for(cutoff in c(1,0.1,0.01,0.001,0.0001)){                p = p1[p2<=cutoff]                x = -log10(seq(from = 0,to = 1,length.out = length(p)+1))[-1]                print(max(x))                y = sort(-log10(p),decreasing = T)                dat = rbind(dat,data.frame(x=x,y=y,cutoff))        }        dat$cutoff = factor(dat$cutoff)        p = ggplot(dat,aes(x=x,y=y,fill=cutoff,colour=cutoff)) +                geom_line(size=1.2) +                geom_abline(intercept=0,slope=1) +                labs(title = t) +                labs(x = xlab) +                labs(y = ylab) +                ylim(0,log10(length(p1)))+                scale_fill_manual(values=mycol)+                theme_bw()+ theme(                        panel.grid=element_blank(),                        axis.text.x=element_text(colour="black", hjust=1),                        axis.text.y=element_text(colour="black"),                        panel.border=element_rect(colour = "black"),                        legend.key = element_blank(),                        plot.title = element_text(hjust = 0.5))}######################################################################################setwd("/share/nas1/huangls/project/zx-20210914-383-gwas_cfdr") t1d<-read.table("t1d_gwas_buildGRCh48_pavalue_removedld.tsv",header = F,sep = "\t") fnbmd<-read.table("fn-bmd-gwas_grch48_pvalue_removedld.tsv",header = F,sep="\t") lada<-read.table("lada_gwas_grch48_pvalue_removedld.tsv",header = F,sep="\t")  #t1d<-read.table("t1d_gwas_buildGRCh48_pavalue.tsv",header = F,sep = "\t") #fnbmd<-read.table("fn-bmd-gwas_grch48_pvalue.tsv",header = F,sep="\t") #lada<-read.table("lada_gwas_grch48_pvalue.tsv",header = F,sep="\t")  colnames(t1d)<-c("chr","pos","ref","alt","t1d.p","variant_id","AF","id") colnames(fnbmd)<-c("chr","pos","fnbmd.p","id") colnames(lada)<-c("chr","pos","lada.p","id")  aa=inner_join(t1d,fnbmd,by="id") df<-inner_join(aa,lada,by="id") #df$t1d.p[df$t1d.p==0]=1.5e-300 #df$fnbmd.p[df$fnbmd.p==0]=1.5e-300 #df$lada.p[df$lada.p==0]=1.5e-300 # #cfdr https://annahutch.github.io/fcfdr/articles/t1d_app.html#   af=df$AF df$maf=ifelse(af>0.5,1-af,af)  df<-df[nchar(df$ref)==1 & nchar(df$alt)==1,  ] df<-df[!(df$fnbmd.p==1 | df$lada.p==1 | df$t1d.p==1),  ]  df=df[,c("chr.x"  , "pos.x"  , "ref"  ,   "alt" ,  "id","variant_id","AF", "fnbmd.p"  ,"t1d.p"  , "lada.p" )]  write.table(df,"all.pvalue.tsv",quote = F,row.names = F,sep = "\t")  #xx=read.table("all.pvalue.tsv",header = T,sep = "\t") library(RColorBrewer) palette(brewer.pal(8,"Set1"))################################################################## # https://www.biorxiv.org/content/10.1101/2020.12.04.411710v2.full.pdf  pp=c(1,0.1,0.01,0.001,0.0001) PP<- c("fnbmd.t1d.cfdr","fnbmd.lada.cfdr")  TT<- c("FNBMD|T1D","FNBMD|LADA")  dd=NULL for(j in 1:length(PP)){         t=unlist(strsplit(PP[j],"[.]"))         # p1=df[,paste0(t[1],".p")]         # p2=df[,paste0(t[2],".p")]         p1=paste0(t[1],".p")         p2=paste0(t[2],".p")       ###################QQQQ##################################         titl=paste0(toupper(t[1]),"|",toupper(t[2]))         p = stratifiedQQplot(df[,p1],df[,p2],                              ylab =bquote(Nominal~~-log[10](P[.(titl)])) ,                              xlab =bquote(Empirical~~-log[10](q[.(titl)])) ,                              t=titl)         png(filename=paste0(t[1],".",t[2],".qq.png"), height=5*300, width=6*300, res=300, units="px")         print(p)         dev.off()         pdf(file=paste0(t[1],".",t[2],".qq.pdf"), height=5, width=6)         print(p)         dev.off()                  titl=paste0(toupper(t[2]),"|",toupper(t[1]))         p = stratifiedQQplot(df[,p2],df[,p1],                              ylab =bquote(Nominal~~-log[10](P[.(titl)])) ,                              xlab =bquote(Empirical~~-log[10](q[.(titl)])) ,                              t=titl)         png(filename=paste0(t[2],".",t[1],".qq.png"), height=5*300, width=6*300, res=300, units="px")         print(p)         dev.off()         pdf(file=paste0(t[2],".",t[1],".qq.pdf"), height=5, width=6)         print(p)         dev.off()###########################cFDR##########################################         cf1=cfdr::cfdr(df[,p1],df[,p2])         cf2=cfdr::cfdr(df[,p2],df[,p1])         #cf1=GWAScFDR::cFDR(df[,p1],df[,p2])         #cf2=GWAScFDR::cFDR(df[,p2],df[,p1])         cf=data.frame(cFDR1=cf1,cFDR2=cf2)         cf$ccFDR=apply(cf, 1, max)         #ccf=condFDR::ccFDR(df,p1,p2,p_threshold  = 0.1,mc.cores = 8)         res=data.frame(df,cf)         res=plyr::rename(res,c("cFDR1"=paste0(toupper(t[1]),"|",toupper(t[2]),".cFDR"),"cFDR2"=paste0(toupper(t[2]),"|",toupper(t[1]),".cFDR"),"ccFDR"=paste0(toupper(t[1]),"|",toupper(t[2]),".ccFDR")))         write.table(res,paste0(t[1],".",t[2],".cfdr.all.res.tsv"),quote = F,row.names = F,sep = "\t") }# merge all data   fnbmd.lada=read.table("fnbmd.lada.cfdr.all.res.tsv",head=T,sep="\t",check.names = F) fnbmd.t1d=read.table("fnbmd.t1d.cfdr.all.res.tsv",head=T,sep="\t",check.names = F)  fnbmd.t1d=fnbmd.t1d[,c("id","FNBMD|T1D.cFDR",  "T1D|FNBMD.cFDR",  "FNBMD|T1D.ccFDR")] fnbmd.cfdr.res=inner_join(fnbmd.lada,fnbmd.t1d,by="id",check.names = F) head(fnbmd.cfdr.res)write.table(fnbmd.cfdr.res,"cfdr.all.res.tsv",quote = F,row.names = F,sep = "\t")#做ANNOVAR注释之后合并aa=read.table("cfdr.hg38_multianno.txt",head=T,sep="\t",check.names = F)bb=read.table("cfdr.all.res.tsv",head=T,sep="\t",check.names = F) aa$id=paste(aa$Chr,aa$Start,sep = "-")aa=aa[,c("id","avsnp150","Func.refGene"  ,      "Gene.refGene" ,   "GeneDetail.refGene", "ExonicFunc.refGene", "AAChange.refGene", "cytoBand")]cfdr.ann=inner_join(bb,aa,by="id",check.names = F)write.table(cfdr.ann,"cfdr.all.anno.res.tsv",quote = F,row.names = F,sep = "\t")

ANNOVAR注释:

perl /share/work/biosoft/annovar/2019Oct24/table_annovar.pl \ --buildver hg38 \ --outfile cfdr \ --remove \ --protocol refGene,cytoBand,avsnp150\ --operation g,r,f \ --nastring . \ cfdr.var.ann.avinput /share/work/database/ref/Homo_sapiens/humandb/hg38/

关于"GWAS cFDR多效基因实例分析"的内容就介绍到这里了，感谢大家的阅读。如果想了解更多行业相关的知识，可以关注行业资讯频道，小编每天都会为大家更新不同的知识点。