From 0a19158d3896c85f3ffb433fc469619ac54f5bcd Mon Sep 17 00:00:00 2001 From: Hongbo Liu Date: Tue, 29 Mar 2022 01:36:42 -0400 Subject: [PATCH] Add files via upload --- Figure_Code/Code.For.Figure3.R | 33 ++++++++++++++++++ Figure_Code/Code.For.Figure4.R | 62 ++++++++++++++++++++++++++++++++++ Figure_Code/Code.For.Figure5.R | 58 +++++++++++++++++++++++++++++++ Figure_Code/Code.For.Figure6.R | 35 +++++++++++++++++++ Figure_Code/Code.For.Figure7.R | 62 ++++++++++++++++++++++++++++++++++ 5 files changed, 250 insertions(+) create mode 100644 Figure_Code/Code.For.Figure3.R create mode 100644 Figure_Code/Code.For.Figure4.R create mode 100644 Figure_Code/Code.For.Figure5.R create mode 100644 Figure_Code/Code.For.Figure6.R create mode 100644 Figure_Code/Code.For.Figure7.R diff --git a/Figure_Code/Code.For.Figure3.R b/Figure_Code/Code.For.Figure3.R new file mode 100644 index 0000000..1608a78 --- /dev/null +++ b/Figure_Code/Code.For.Figure3.R @@ -0,0 +1,33 @@ +########################################################################### +##### Code for Figure.3 ##### +##### Hongbo Liu (hongbo919@gmail.com) ##### +########################################################################### + + +########################################################################### +##### Code for Figure.3b ##### +########################################################################### +##### Note: Manhattan plot of human kidney meQTL data (n=443). +library(qqman) + +meQTL <- read.table("meQTL.forManhattanPlot.txt", header=T, sep="\t"); colnames(meQTL) <- c("CHR","POS","SNP","PVAL") +jpeg('Figure.3c.jpg', width = 20, height = 5, units = 'in', res = 500, quality = 100) +manhattan(meQTL, chr="CHR", bp="POS", snp="SNP", p="PVAL", col = c("deepskyblue4", "goldenrod3"), suggestiveline = F, genomewideline = F, cex = 0.1, main = "Manhattan Plot") +dev.off() + + +########################################################################### +##### Code for Figure.3e ##### +########################################################################### +##### Note: Enrichment of tissue-specific meQTL CpGs in tissue enhancers +library(ggplot2) +Tissue.Specific.meQTL.CpG.2HMM <- read.table("Tissue.Specific.meQTL.CpG.2HMM.txt", header=T, sep="\t"); + +pdf(Figure.3d.pdf",width=5,height=3) +ggplot(data = Tissue.Specific.meQTL.CpG.2HMM, aes(x=Enhancer, y=meQTL, fill=OddsRatio)) + + geom_tile(color = "white", size = 1)+ + scale_fill_gradient2(low = "#3C5488FF", high = "#DC0000FF", mid = "white", + midpoint = 1, limit = c(0.5,1.52), space = "Lab", + name="") + + theme_minimal() +dev.off() diff --git a/Figure_Code/Code.For.Figure4.R b/Figure_Code/Code.For.Figure4.R new file mode 100644 index 0000000..b2dc2ee --- /dev/null +++ b/Figure_Code/Code.For.Figure4.R @@ -0,0 +1,62 @@ +########################################################################### +##### Code for Figure.4 ##### +##### Hongbo Liu (hongbo919@gmail.com) ##### +########################################################################### + + +########################################################################### +##### Code for Figure.4b ##### +########################################################################### +##### Note: Estimated proportion of heritability (h2med/h2g) mediated by DNA methylation and gene expression in 414 kidneys across 34 GWAS traits. +library(dplyr) +library(ggplot2) +library(Rmisc) +library(ggrepel) + +read.table("h2med.Estimate_over_h2.Combined.meQTL.eQTL.txt", header=T, sep="\t")-> h2med_over_h2 +read.table("h2med.Estimate_over_h2_SE.Combined.meQTL.eQTL.txt", header=T, row.names = 1, sep="\t")-> h2med_over_h2_SE + +pdf(file='Figure.4b.pdf', height = 4, width = 4) +ggplot(h2med_over_h2, aes(y=meQTL.h2med_over_h2, x=eQTL.h2med_over_h2)) + theme(legend.position="none") + +geom_pointrange(aes(xmax = eQTL.h2med_over_h2 + eQTL.h2med_over_h2.SE, xmin = eQTL.h2med_over_h2 - eQTL.h2med_over_h2.SE, color = TraitType)) + +geom_pointrange(aes(ymax = meQTL.h2med_over_h2 + meQTL.h2med_over_h2.SE, ymin = meQTL.h2med_over_h2 - meQTL.h2med_over_h2.SE, color = TraitType)) + +geom_point(aes(color = TraitType, size = 3)) + scale_color_manual(values=c("gray57", "darkgreen", "gold3", "dodgerblue", "magenta")) + geom_abline(intercept = 0, slope = 1, linetype="dashed", color = "gray") + +scale_x_continuous(limits = c(-0.08,0.7), breaks=seq(0, 0.6, by = 0.2)) + +scale_y_continuous(limits = c(-0.08,0.7),breaks=seq(0, 0.6, by = 0.2)) + +theme_bw(base_size = 12) + +theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),legend.position = c(2, 2)) +dev.off() + + +########################################################################### +##### Code for Figure.4e ##### +########################################################################### +##### Note: Enrichment of kidney methylation-mediated heritability for kidney function GWAS traits (x-axis) in different ChromHMM regulatory elements (y-axis). +##### White to red indicates h2med enrichment (nominal two-sided p < 0.05 calculated by MESC). +##### Asterisk indicates h2med enrichment passing FDR q < 0.05 (accounting for 374 tests for 11 chromatin state CpG sets and 34 GWAS traits). +library(ggplot2) +library(gplots) +h2med_enrichment <- read.table("Kidney.meQTL.Heritability.HMM.h2med_enrichment.txt",header=T, sep="\t") +Kidney.CpGset.h2med.HMM_State <- read.table("Kidney.CpGset.h2med.HMM_State.txt",header=T, sep="\t") + +# A 4 x 4 matrix of characters (stars) +labs <- h2med_enrichment +for (i in 1:nrow(labs)) { + for (j in 1:ncol(labs)) { + labs[i,j] <- "" + if (subset(Kidney.CpGset.h2med.HMM_State, Gene_category == rownames(labs)[i] & GWAS == colnames(labs)[j])$h2med_enrichment_qvalue < 0.05 & h2med_enrichment[i,j] > 0){ + labs[i,j] <- "*" + } + } +} + +my_palette <- colorRampPalette(c("gray", "white", "red"))(n = 74) +pdf("Figure.4e.pdf", width = 7, height = 8) +x<- as.matrix(h2med_enrichment); +heatmap <- heatmap.2(x, scale = "none", col=my_palette, trace = "none", density.info = "none", Rowv=FALSE, Colv=FALSE, cellnote=labs, notecol="black", notecex=2, margins = c(20, 20)) +dev.off() + + + + + diff --git a/Figure_Code/Code.For.Figure5.R b/Figure_Code/Code.For.Figure5.R new file mode 100644 index 0000000..436ce13 --- /dev/null +++ b/Figure_Code/Code.For.Figure5.R @@ -0,0 +1,58 @@ +########################################################################### +##### Code for Figure.5 ##### +##### Hongbo Liu (hongbo919@gmail.com) ##### +########################################################################### + + +########################################################################### +##### Code for Figure.5b ##### +########################################################################### +##### Note: Enrichment of kidney methylation-mediated heritability for kidney function GWAS traits (x-axis) in kidney cell type-specific accessible regions (y-axis). +##### White to red indicates h2med enrichment (nominal p < 0.05 calculated by MESC). +##### Asterisk indicates h2med enrichment passing FDR q < 0.05 (accounting for 408 tests for 12 cell type CpG sets and 34 GWAS traits). +library(ggplot2) +library(gplots) +h2med_enrichment <- read.table("Kidney.meQTL.Heritability.snATAC_Cluster.h2med_enrichment.txt",header=T, sep="\t") +Kidney.CpGset.h2med.snATAC_Cluster <- read.table("Kidney.CpGset.h2med.snATAC_Cluster.txt",header=T, sep="\t") + +# A 4 x 4 matrix of characters (stars) +labs <- h2med_enrichment +for (i in 1:nrow(labs)) { + for (j in 1:ncol(labs)) { + labs[i,j] <- "" + if (subset(Kidney.CpGset.h2med.snATAC_Cluster, Gene_category == rownames(labs)[i] & GWAS == colnames(labs)[j])$h2med_enrichment_qvalue < 0.05 & h2med_enrichment[i,j] > 0){ + labs[i,j] <- "*" + } + } +} + +my_palette <- colorRampPalette(c("gray", "white", "red"))(n = 74) +pdf("Figure.5e.pdf", width = 7, height = 8) +x<- as.matrix(h2med_enrichment); +heatmap <- heatmap.2(x, scale = "none", col=my_palette, trace = "none", density.info = "none", Rowv=FALSE, Colv=FALSE, cellnote=labs, notecol="black", notecex=2, margins = c(20, 20)) +dev.off() + + + +########################################################################### +##### Code for Figure.5c ##### +########################################################################### +##### Note: Single cell GWAS trait enrichment in human kidney cells by gchromVAR. +library(ggplot2) +library(gplots) +UKBB2KidATAC.snATAC_Cluster.Mean <- read.table("UKBB2KidATAC.snATAC_Cluster.Mean.tsv", header=T, sep=",") + +my_palette <- colorRampPalette(c("blue", "white", "red"))(n = 256) +library("gplots") +pdf("Figure.5c.pdf", width = 7, height = 7) +x<- as.matrix(UKBB2KidATAC.snATAC_Cluster.Mean); +heatmap <- heatmap.2(x, scale = "none", Rowv = FALSE, Colv = FALSE, col=my_palette, trace = "none", density.info = "none", margins = c(7, 7)) +dev.off() + + + + + + + + diff --git a/Figure_Code/Code.For.Figure6.R b/Figure_Code/Code.For.Figure6.R new file mode 100644 index 0000000..d8a3d85 --- /dev/null +++ b/Figure_Code/Code.For.Figure6.R @@ -0,0 +1,35 @@ +########################################################################### +##### Code for Figure.6 ##### +##### Hongbo Liu (hongbo919@gmail.com) ##### +########################################################################### + + +########################################################################### +##### Code for Figure.6c ##### +########################################################################### +##### Note: Manhattan plot highlighting 330 genes with evidence of multiple traits colocalization. +library(ggman) +GWAS <- read.table("GWAS.forManhattanPlot.txt", header=T, sep="\t") +HighlightSNPs <- read.table("HighlightSNPs.txt", header=T, sep="\t") +HighlightGMESNPs <- read.table("HighlightGMESNPs.txt", header=T, sep="\t") + +GWAS[which(GWAS$PVAL < 1e-300),4] = 1e-300 + +p1 <- ggman(GWAS, sigLine = 7.30103, snp = "RSID", bp = "POS", chrom = "CHR", pvalue = "PVAL", relative.positions = TRUE) + theme(text = element_text(size = 20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +p2 <- ggmanHighlightGroup(p1, highlightDfm = HighlightSNPs, snp = "RSID", group = "Group", size = 2, alpha = 1, stroke = 0, legend.title = "Coloc Type") +p3 <- ggmanHighlight(p2, highlight = as.character(HighlightGMESNPs$RSID), size = 2, stroke = 0) + +jpeg('Figure.6c.jpg', width = 15, height = 8, units = 'in', res = 500, quality = 100) +ggmanLabel(p3, labelDfm = Gene_Highlit, snp = "RSID", label = "Gene", type = "text", colour = "black", fontface = 'italic', size = 3.5, + nudge_y = 50, + max.iter = 5000, + vjust=1, + direction='y', + nudge_x=0.1, + segment.color = "grey50", + segment.alpha = 0.5, + segment.size = 0.2) +dev.off() + + + diff --git a/Figure_Code/Code.For.Figure7.R b/Figure_Code/Code.For.Figure7.R new file mode 100644 index 0000000..8b06274 --- /dev/null +++ b/Figure_Code/Code.For.Figure7.R @@ -0,0 +1,62 @@ +########################################################################### +##### Code for Figure.7 ##### +##### Hongbo Liu (hongbo919@gmail.com) ##### +########################################################################### + + +########################################################################### +##### Code for Figure.7a ##### +########################################################################### +##### http://locuszoom.org/genform.php?type=yourdata + + +########################################################################### +##### Code for Figure.7b ##### +########################################################################### +##### Note: Genotype (rs2252281, x-axis) and normalized CpG methylation (cg15971010, y-axis) association in human kidneys (n=443) +library(ggplot2) +library(ggsci) +library(ggpubr) +S443.Methyl.Genotype <- read.table("Chr17.rs2252281.cg15971010.S443.Methyl.Genotype.infor.txt", header=T, sep="\t") +meQTL <- read.table("Chr17.rs2252281.cg15971010.S443.meQTL.infor.txt", header=T, sep="\t") +S443.Methyl.Genotype <- na.omit(S443.Methyl.Genotype) + +pdf(file='Figure.7b.pdf',width=3,height=3.8) +ggboxplot(S443.Methyl.Genotype, x = "rs2252281", y = "cg15971010", fill = "rs2252281", title = paste ("meQTL\nbeta = ", round(meQTL$beta,2), "\np = ", signif(meQTL$pvalue,2), sep = ""), add = "jitter", add.params = list(size = 0.1, jitter = 0.1, color = "black"), legend = "none") + scale_color_npg() + scale_fill_npg() + theme(plot.title = element_text(size=12)) +dev.off() + + +########################################################################### +##### Code for Figure.7c ##### +########################################################################### +##### Note: Genotype (rs2252281 x-axis) and normalized gene expression (SLC47A1, y-axis) association in human kidney tubule samples (n=356). +library(ggplot2) +library(ggsci) +library(ggpubr) +S356.Gex.Genotype <- read.table("Gex.Genotype.infor.txt", header=T, sep="\t") +eQTL <- read.table("eQTL.infor.txt", header=T, sep="\t") +S356.Gex.Genotype <- na.omit(S356.Gex.Genotype) + +pdf(file='Figure.7b.pdf',width=3,height=3.8) +ggboxplot(S356.Gex.Genotype, x = "rs2252281", y = "SLC47A1", fill = "rs2252281", title = paste ("eQTL\nbeta = ", round(eQTL$beta,2), "\np = ", signif(eQTL$pvalue,2), sep = ""), add = "jitter", add.params = list(size = 0.1, jitter = 0.1, color = "black"), legend = "none") + scale_color_npg() + scale_fill_npg() + theme(plot.title = element_text(size=12)) +dev.off() + + +########################################################################### +##### Code for Figure.7d ##### +########################################################################### +##### Note: Manhattan plot of phenome-wide association study of predicted loss-of-function (pLOF) variants in SLC47A1 in UKBB. +library(dplyr) +library(PheWAS) +load("SLC47A1_pLOF_burden_PheWAS.RData") +pdf(file='Figure.7d.pdf',width=7,height=3.2) +phewasManhattan(SLC47A1_pLOF_burden_PheWAS, Phenotypes, OR.direction = , annotate.angle= 35, y.axis.interval = 1, sort.by.category.value=T, annotate.list = annoteLable.pLOF, base.labels =T, annotate.size = 0, title="pLOF burden PheWAS manhattan plot") +dev.off() + + + + + + + +