Skip to content

Commit

Permalink
Converted the rest of the old base plotting code to the tidyverse to …
Browse files Browse the repository at this point in the history 
…make it simpler and programmatic
  • Loading branch information
@ben-laufer
ben-laufer committed Feb 11, 2020
1 parent 8e6e8b1 commit effd82a
Showing 1 changed file with 155 additions and 181 deletions.
336 changes: 155 additions & 181 deletions GATplots.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,12 @@
# Takes the output of GAT for hg38 CpG and genic annotation enrichments and makes a ggplot2
# Ben Laufer

# Install -----------------------------------------------------------------
# Packages ----------------------------------------------------------------

if (!requireNamespace(c("tidyverse", "reshape2"), quietly = TRUE))
install.packages(c("tidyverse", "reshape2"), repos = "https://cloud.r-project.org")
library(tidyverse)
library(reshape2)
if (!requireNamespace(c("tidyverse"), quietly = TRUE))
install.packages("tidyverse", repos = "https://cloud.r-project.org")
suppressPackageStartupMessages(library(tidyverse))
options(readr.num_columns = 0)

# Functions ---------------------------------------------------------------

Expand All @@ -27,106 +27,89 @@ gg_color_hue <- function(n) {

# CpG ---------------------------------------------------------------------

# Load
GAT <- read.delim("Gat_CpG_hyper_hypo_results.tsv")

# Nice names
GAT <- GAT %>%
dplyr::as_tibble() %>%
cat("\n","Tidying CpG annotations...")
GAT <- readr::read_tsv("Gat_CpG_hyper_hypo_results.tsv") %>%
dplyr::mutate(annotation = dplyr::recode_factor(annotation,
"hg38_cpg_inter" = "Open Sea",
"hg38_cpg_shelves" = "CpG Shelf",
"hg38_cpg_islands" = "CpG Island",
"hg38_cpg_shores" = "CpG Shore")
)

GAT.1 <- GAT

# Fix Fold Change
index <- GAT$fold < 1
GAT$fold[index] <- -(1/GAT$fold[index])

# Order
GAT <- dcast(GAT, annotation ~ track, value.var = "fold")
order <- c("CpG Island",
"CpG Shore",
"CpG Shelf",
"Open Sea")
GAT$annotation <- factor(GAT$annotation,levels = order)
GAT <- GAT[match(order, GAT$annotation),]
row.names(GAT) <- c(1:nrow(GAT))
GAT$annotation <- factor(GAT$annotation, levels = rev(order))

# Colors for Custom Order
n = nrow(GAT)
cols = gg_color_hue(n)
names(cols) <- order

# Melt
GAT <- melt(GAT)

# Significance
GAT.1 <- dcast(GAT.1, annotation ~ track, value.var = "qvalue")
GAT.1$annotation <- factor(GAT.1$annotation, levels = order)
GAT.1 <- GAT.1[match(order, GAT.1$annotation),]
row.names(GAT.1) <- c(1:nrow(GAT.1))
GAT.1$annotation <- factor(GAT.1$annotation, levels = rev(order))
GAT.1 <- melt(GAT.1)
GAT.1$signif <- ifelse(GAT.1$value < 0.05,1,0)
GAT.1 <- GAT.1[-3]
GAT <- merge(GAT, GAT.1, by = c("annotation", "variable"))

# Plot
CpG <- ggplot(data = GAT, aes(annotation,
y = value,
fill = annotation)
) +
geom_bar(stat = "identity") +
coord_flip() +
labs(y = "Fold Enrichment",
x = element_blank()) +
theme_classic() +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 16),
strip.text = element_text(size = 16),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14)
) %>%
dplyr::mutate(annotation = forcats::fct_relevel(annotation,
"Open Sea",
"CpG Shelf",
"CpG Shore",
"CpG Island")
) %>%
dplyr::mutate(fold = dplyr::case_when(fold < 1 ~ -1/fold,
fold >= 1 ~ fold)
) %>%
dplyr::mutate(signif = dplyr::case_when(qvalue <= 0.05 ~ 1,
qvalue > 0.05 ~ 0)
) %>%
dplyr::select(annotation,
track,
fold,
signif)
cat("Done.")

cat("\n", "Choosing colors...")
cols = nlevels(GAT$annotation) %>%
gg_color_hue() %>%
rev()

names(cols) <- levels(GAT$annotation)
order <- rev(levels(GAT$annotation))
cat("Done.")

cat("\n", "Plotting...")
(ggplot(data = GAT,
aes(annotation,
y = fold,
fill = annotation)
) +
scale_y_continuous(limits = c(-4,4),
breaks = c(-4,-3,-2,-1,0,1,2,3,4)
) +
scale_fill_manual(values = cols,
breaks = order,
name = "Annotation") +
facet_grid(~variable) +
geom_hline(yintercept = 0) +
geom_text(data = GAT[(GAT$signif == 1 & GAT$value > 0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = 0.5,
nudge_x = -0.05) +
geom_text(data = GAT[(GAT$signif == 1 & GAT$value < 0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = -0.5,
nudge_x = -0.05)

ggsave("Gat_CpG_hyper_hypo.pdf",
plot = CpG,
device = NULL,
width = 11,
height = 6)
geom_bar(stat = "identity") +
coord_flip() +
labs(y = "Fold Enrichment",
x = element_blank()) +
theme_classic() +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 16),
strip.text = element_text(size = 16),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14)
) +
scale_y_continuous(limits = c(-4,4),
breaks = c(-4,-3,-2,-1,0,1,2,3,4)
) +
scale_fill_manual(values = cols,
breaks = order,
name = "Annotation") +
facet_grid(~track) +
geom_hline(yintercept = 0) +
geom_text(data = GAT[(GAT$signif == 1 & GAT$fold > 0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = 0.5,
nudge_x = -0.05) +
geom_text(data = GAT[(GAT$signif == 1 & GAT$fold < 0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = -0.5,
nudge_x = -0.05)
) %>%
ggsave("Gat_CpG_hyper_hypo.pdf",
plot = .,
width = 11,
height = 6)
cat("Done.")

# Genic ------------------------------------------------------------------

# Load
GAT2 <- read.delim("Gat_genic_hyper_hypo_results.tsv")

# Nice names
GAT2 <- GAT2 %>%
dplyr::as_tibble() %>%
cat("\n","\n", "Tidying genic annotations...")
GAT2 <- readr::read_tsv("Gat_genic_hyper_hypo_results.tsv") %>%
dplyr::mutate(annotation = dplyr::recode_factor(annotation,
"enhancers_merged.bed" = "Enhancers",
"intergenic_merged.bed" = "Intergenic",
Expand All @@ -137,89 +120,80 @@ GAT2 <- GAT2 %>%
"exons_merged.bed" = "Exon",
"fiveUTRs_merged.bed" = "5' UTR",
"promoters_merged.bed" = "Promoter")
)

GAT2.1 <- GAT2

# Fix Fold Change
index1 <- GAT2$fold < 1
GAT2$fold[index1] <- -(1/GAT2$fold[index1])

# Order
GAT2 <- dcast(GAT2, annotation ~ track, value.var = "fold")
order <- c("Enhancers",
"1-5 kb Upstream",
"Promoter",
"5' UTR",
"Exon",
"Exon/Intron Boundaries",
"Intron",
"3' UTR",
"Intergenic")
GAT2$annotation <- factor(GAT2$annotation, levels = order)
GAT2 <- GAT2[match(order, GAT2$annotation),]
row.names(GAT2) <- c(1:nrow(GAT2))
GAT2$annotation <- factor(GAT2$annotation, levels = rev(order))

# Colors for Custom Order
n = nrow(GAT2)
cols = gg_color_hue(n)
names(cols) <- order

# Melt
GAT2 <- melt(GAT2)

# Significance
GAT2.1 <- dcast(GAT2.1, annotation ~ track, value.var = "qvalue")
GAT2.1$annotation <- factor(GAT2.1$annotation,levels = order)
GAT2.1 <- GAT2.1[match(order, GAT2.1$annotation),]
row.names(GAT2.1) <- c(1:nrow(GAT2.1))
GAT2.1$annotation <- factor(GAT2.1$annotation, levels = rev(order))
GAT2.1 <- melt(GAT2.1)
GAT2.1$signif <- ifelse(GAT2.1$value < 0.05,1,0)
GAT2.1 <- GAT2.1[-3]
GAT2 <- merge(GAT2, GAT2.1, by = c("annotation", "variable"))

# Plot
gene <- ggplot(data = GAT2, aes(x = annotation,
y = value,
fill = annotation)
) +
geom_bar(stat = "identity") +
coord_flip() +
labs(y = "Fold Enrichment",
x = element_blank()
) +
theme_classic() +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 16),
strip.text = element_text(size = 16),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14)
) %>%
dplyr::mutate(annotation = forcats::fct_relevel(annotation,
"Intergenic",
"3' UTR",
"Intron",
"Exon/Intron Boundaries",
"Exon",
"5' UTR",
"Promoter",
"1-5 kb Upstream",
"Enhancers")
) %>%
dplyr::mutate(fold = dplyr::case_when(fold < 1 ~ -1/fold,
fold >= 1 ~ fold)
) %>%
dplyr::mutate(signif = dplyr::case_when(qvalue <= 0.05 ~ 1,
qvalue > 0.05 ~ 0)
) %>%
dplyr::select(annotation,
track,
fold,
signif)
cat("Done.")

cat("\n", "Choosing colors...")
cols = nlevels(GAT2$annotation) %>%
gg_color_hue() %>%
rev()

names(cols) <- levels(GAT2$annotation)
order <- rev(levels(GAT2$annotation))
cat("Done.")

cat("\n", "Plotting...")
(ggplot(data = GAT2,
aes(x = annotation,
y = fold,
fill = annotation)
) +
scale_y_continuous(limits = c(-3.5,3.5),
breaks = c(-3,-2,-1,0,1,2,3)
) +
scale_fill_manual(values = cols,
breaks = order,
name = "Annotation") +
facet_grid(~variable) +
geom_hline(yintercept = 0) +
geom_text(data = GAT2[(GAT2$signif ==1 & GAT2$value >0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = 0.5,
nudge_x = -0.09) +
geom_text(data = GAT2[(GAT2$signif ==1 & GAT2$value <0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = -0.5,
nudge_x = -0.09)

ggsave("Gat_genic_hyper_hypo.pdf",
plot = gene,
device = NULL,
width = 11,
height = 6)
geom_bar(stat = "identity") +
coord_flip() +
labs(y = "Fold Enrichment",
x = element_blank()
) +
theme_classic() +
theme(axis.text = element_text(size = 16),
axis.title = element_text(size = 16),
strip.text = element_text(size = 16),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14)
) +
scale_y_continuous(limits = c(-3.5,3.5),
breaks = c(-3,-2,-1,0,1,2,3)
) +
scale_fill_manual(values = cols,
breaks = order,
name = "Annotation") +
facet_grid(~track) +
geom_hline(yintercept = 0) +
geom_text(data = GAT2[(GAT2$signif ==1 & GAT2$fold >0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = 0.5,
nudge_x = -0.09) +
geom_text(data = GAT2[(GAT2$signif ==1 & GAT2$fold <0), ],
label = "*",
size = 8,
show.legend = FALSE,
nudge_y = -0.5,
nudge_x = -0.09)
) %>%
ggsave("Gat_genic_hyper_hypo.pdf",
plot = .,
width = 11,
height = 6)
cat("Done.", "\n", "\n")

0 comments on commit effd82a

Please sign in to comment.