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_buildGRCh38_pavalue_removedld.tsv",header = F,sep = "\t")
fnbmd<-read.table("fn-bmd-gwas_grch38_pvalue_removedld.tsv",header = F,sep="\t")
lada<-read.table("lada_gwas_grch38_pvalue_removedld.tsv",header = F,sep="\t")
#t1d<-read.table("t1d_gwas_buildGRCh38_pavalue.tsv",header = F,sep = "\t")
#fnbmd<-read.table("fn-bmd-gwas_grch38_pvalue.tsv",header = F,sep="\t")
#lada<-read.table("lada_gwas_grch38_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/
