您好,登录后才能下订单哦!
密码登录
登录注册
点击 登录注册 即表示同意《亿速云用户服务条款》
这篇文章主要介绍“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多效基因实例分析”的内容就介绍到这里了,感谢大家的阅读。如果想了解更多行业相关的知识,可以关注亿速云行业资讯频道,小编每天都会为大家更新不同的知识点。
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。