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多效基因實例分析”的內容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業相關的知識，可以關注億速云行業資訊頻道，小編每天都會為大家更新不同的知識點。