Fig1b

Contents

• Introduction
• FHS
• Hazard Rate Lookup Table
• Proportion Variation Explained Lookup Table
• Merge files for joint panel

Introduction

FHS

### plot with PVE
setwd("~/dynamicHRpaper/")
fig_path="./Figs/Fig1b/"
f3=readRDS("data/fhs_prs_with_pc_corrrection.rds")
df=readRDS("data/fh_full.rds")

dat=data.frame(df %>%group_by(round(AGE1,0)) %>%summarise(length(AGE1)))
dat=dat[c(11:50),]


#hazards=readRDS("output/hazards_.rds")
hazards=data.frame(readRDS("output/hazards_fh_EQ_withPC_corrections.rds"))
hazards=data.frame(hazards)

hazards=data.frame(hazards)
colnames(hazards)=c("prs.HR","prs.OR","prs.r2",
                    
                    "ldl.HR","ldl.OR","ldl.r2",
                    
                    "smoke.HR","smoke.OR","smoke.r2",
                    "age.HR","age.OR","age.r2",
                    "sbp.HR","sbp.OR","sbp.r2",
                    #"dm.HR","dm.OR","dm.r2",
                    "chol.HR","chol.OR","chol.r2",
                    "bmi.HR","bmi.OR","bmi.r2",
                    "sex.HR","sex.OR","sex.r2",
                    
                    "hdl.HR","hdl.OR","hdl.r2")




hrmat=hazards[,grep("HR",colnames(hazards))];hrmat$age=dat$round.AGE1..0.

ormat=hazards[,grep("OR",colnames(hazards))];ormat$age=dat$round.AGE1..0.
r2mat=hazards[,grep("r2",colnames(hazards))];r2mat$age=dat$round.AGE1..0.

colnames(hrmat)=colnames(ormat)=colnames(r2mat)=c("Polygenic risk score","LDL Cholesterol","Smoking","Age","Systolic blood pressure","Total cholesterol","Body Mass Index","Male Sex","HDL cholesterol","age")

#####

setwd("~/dynamicHRpaper/")

upperbound=data.frame(readRDS("output/hazards_rate_upperbound_fh_EQ_withPC_corrections.rds"))
lowerbound=data.frame(readRDS("output/hazards_rate_lowerbound_fh_EQ_withPC_corrections.rds"))

upperbound$age=dat$round.AGE1..0.
u=melt(upperbound,id.vars="age")

lowerbound$age=dat$round.AGE1..0.
l=melt(lowerbound,id.vars="age")


##### Bring together
haz=melt(hrmat,id.vars = "age")
ors=melt(ormat,id.vars = "age")
r2s=melt(r2mat,id.vars = "age")

haz$low.se = l$value
haz$upper.se = u$value


s=reshape(haz, direction = "wide",
idvar = "age", timevar = "variable")


r=reshape(r2s, direction = "wide",
idvar = "age", timevar = "variable")
write.csv(s,"output/hrindex_fhs_eq_with_pc.csv",row.names = F)

write.csv(r,"output/r2index_fhs_eq_with_pc.csv",row.names = F)

r2s$variable=ors$variable=haz$variable <- factor(haz$variable, levels = c("Total cholesterol","Male Sex","Smoking","Polygenic risk score","LDL Cholesterol","Systolic blood pressure","Body Mass Index","Age","HDL cholesterol"))

hr_fhs = ggplot(subset(
  haz,
  variable %in% c(
    "Polygenic risk score",
    "Total cholesterol",
    "Systolic blood pressure",
    "Pooled cohort equations",
    #"Age",
    "Smoking",
    "Male Sex",
    "HDL cholesterol",
    "Smoking")), aes(age, value, colour = variable,fill=variable)) +
  #geom_point(size = 1) + 
  stat_smooth(method = "loess") + labs(y = "Hazard Ratio of Coronary Artery Disease", x ="Age of Risk Assessment", colour = "Risk Factor") + scale_color_manual(values = col_p,drop = TRUE,breaks = c("Diabetes mellitus","Male Sex","Smoking","Polygenic risk score","LDL Cholesterol","Systolic blood pressure",
"Body Mass Index","Pooled cohort equations","Age","HDL cholesterol",
"Total cholesterol"),guide="none")+
  scale_fill_manual(values=col_p, drop=TRUE, breaks=c("Diabetes mellitus","Male Sex","Smoking","Polygenic risk score","LDL Cholesterol","Systolic blood pressure","Body Mass Index","Pooled cohort equations","Age","HDL cholesterol","Total cholesterol"),guide="none")


r2_fhs = ggplot(subset(
  r2s,
  variable %in% c(
    "Polygenic risk score",
    "Total cholesterol",
    "Systolic blood pressure",
    "Pooled cohort equations",
    #"Age",
    "Smoking",
    "Male Sex",
    "HDL cholesterol",
    "Smoking")
), aes(age, value, colour = variable,fill=variable))+
  #geom_point(size = 1) +
  stat_smooth(method = "loess") + labs(y = "Proportion of Variation Explained", x ="Age of Risk Assessment", colour = "Risk Factor") + scale_color_manual(values = col_p,drop = TRUE,breaks = c("Diabetes mellitus","Male Sex","Smoking","Polygenic risk score","LDL Cholesterol","Systolic blood pressure","Body Mass Index","Pooled cohort equations","Age","HDL cholesterol","Total cholesterol"))+
  scale_fill_manual(values=col_p, drop=TRUE, breaks=c("Diabetes mellitus","Male Sex","Smoking","Polygenic risk score","LDL Cholesterol","Systolic blood pressure","Body Mass Index","Pooled cohort equations","Age","HDL cholesterol","Total cholesterol"),guide="none")

h1=hr_fhs+theme_classic(base_size = 20)+tems
# 
# ggsave(h1,file="hr_fhs_nob.tiff",dpi = 300,width = 6,height = 6)

r2=r2_fhs+theme_classic(base_size = 20)+theme(legend.position = "none")+tems
#saveRDS(r2,"output/r2fhs_withPC.rds")
h1=h1+theme_classic(base_size = 20)+theme(legend.position = "none")+tems
#saveRDS(h1,"output/h1hs_withPC.rds")

#ggsave(h1,file=paste0(fig_path,"hr_fhs.tiff"),dpi = 300,width = 6,height = 6)
#ggsave(r2,file=paste0(fig_path,"r2_fhs.tiff"),dpi = 300,width = 6,height = 6.1)

ggplotly(h1)

ggplotly(r2)

Hazard Rate Lookup Table

library(DT)
m=cbind(haz[,c(1:2)],round(haz[,c(3:5)],2))
colnames(m)=colnames(haz)
DT::datatable(m)

Proportion Variation Explained Lookup Table

library(DT)
m=cbind(r2s[,c(1:2)],round(r2s[,c(3)],2))
colnames(m)=colnames(r2s)
DT::datatable(m)

Merge files for joint panel

#hall=readRDS("output/hall.rds")
#r2u=readRDS("output/r2.rds")
setwd("~/Dropbox/rootfilesforbackup/dynamichr/")
hall=readRDS("output/hall_withPC.rds")
r2u=readRDS("output/r2_withPC.rds")
h1=readRDS("output/h1hs_withPC.rds")
r2fhs=readRDS("output/r2fhs_withPC.rds")
r2_fhs=r2fhs+ggtitle("FOS")+theme_classic(base_size=15)+labs(x="Age of Risk Assessment, years")+theme(legend.text = element_text(size = 15), legend.title = element_text(size = 15), axis.text.x = element_text(size = 15),axis.text.y = element_text(size = 15))+guides(color=guide_legend(override.aes=list(fill=NA)))

r2u=r2u+ggtitle("UKB")+theme_classic(base_size=15)+labs(x="Age of Risk Assessment, years")+theme(legend.text = element_text(size = 15), legend.title = element_text(size = 15), axis.text.x = element_text(size = 15),axis.text.y = element_text(size = 15))+guides(color=guide_legend(override.aes=list(fill=NA)))

h1=h1+ggtitle("FOS")+theme_classic(base_size=15)+ylim(0,5)+labs(x="Age of Risk Assessment, years",y="Hazard Ratio (95% CI) of CAD")+theme(legend.text = element_text(size = 15), legend.title = element_text(size = 15), axis.text.x = element_text(size = 15),axis.text.y = element_text(size = 15))+guides(color=guide_legend(override.aes=list(fill=NA)))

hall = hall + ggtitle("UKB") + theme_classic(base_size=15) + ylim(0,5) + labs(x="Age of Risk Assessment, years") + theme(legend.text = element_text(size = 15), legend.title = element_text(size = 15), axis.text.x = element_text(size = 15),axis.text.y = element_text(size = 15))+guides(color=guide_legend(override.aes=list(fill=NA)))


print1=ggarrange(h1,hall+rremove("ylab"), ncol=2, common.legend = TRUE, legend="bottom",labels = c("A.","B."),heights=c(3.5,3.5))

print2=ggarrange(r2_fhs,r2u+rremove("ylab"), ncol=2, common.legend = TRUE, legend="bottom",labels = c("A.","B."))


#ggsave(print1,filename = "Figs/Fig2/Fig1combined_eq_new_pc.pdf",dpi = 600,width = 12,height=8)

#ggsave(print2,filename = "Figs/Fig2/Fig2combined_eq_new_pc.pdf",dpi = 600,width = 12,height=8)