First, basic parameters are set up in this RMarkdown, such as loading dependencies, setting paths and setting up a uniform plot structure.
knitr::opts_chunk$set(warning=FALSE, cache = TRUE, message = FALSE)
require(tidyverse)
require(edgeR)
require(readxl)
require(readr)
require(RColorBrewer)
require(DT)
source("HelperFunctions.R")
# plot style
theme_point<-theme_bw()+theme(strip.background = element_blank())
theme_bar<-theme_bw()+theme(axis.text.x = element_text(angle=90, vjust=0.5, hjust=1),strip.background = element_blank(),axis.ticks.x = element_blank(),axis.line.x = element_blank())
theme_boxplot<-theme_bw()+theme(axis.text.x = element_text(angle=90, vjust=0.5, hjust=1),strip.background = element_blank(),axis.ticks.x = element_blank(),axis.title.x = element_blank(),axis.line.x = element_blank(),legend.position = "none")
color_panel<-c("#e35d6a","#5bb75b","#428bca","#e87810","#23496b","#ffbf00","#cc2028","#039748","pink","gray","darkgray")
color_panel1 <- c("#039748","#039748","#ffbf00","#ffbf00","#e35d6a","#e35d6a","#5bb75b","#5bb75b","#428bca","#428bca","#23496b","#23496b","#cc2028","#cc2028","#e87810")
color_panel2 <- brewer.pal(6, name = "Paired")
names(color_panel2) <- c("Biomatrica", "EDTA","RNA Streck","Roche","DNA Streck", "PAXgene")
makeRes <- function(fit, contrast){
lrt <- glmLRT(fit, contrast=contr)
dt1 <- datatable(as.data.frame(topTags(lrt)))
dt2 <- datatable(as.data.frame(summary(decideTests(lrt))))
plot <- plotMD(lrt)
return(list(dt1, dt2, plot))
}
# Global variables
data_path <- "../data/"
readcutoff <- 15
sample_annotation<-read_excel("./annotation_RVP1808.xlsx")
Differential methylation of promotor regions was performed with EdgeR (https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf), and FRY/KEGG according to https://www.bioconductor.org/packages/devel/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html
The cut-off was set at 15 reads per region.
if (file.exists("fitted_DMR_full.RData")){
load("fitted_DMR_full.RData")
} else {
if (file.exists("bismark2DGE_object_full.Rdata")){
load("bismark2DGE_object_full.Rdata")
} else {
files<-list.files(data_path,recursive=TRUE)
files<-files[grep("18.cov$", files)]
yall <- readBismark2DGE(paste0(data_path,sample_annotation$cov_filename), sample.names=sample_annotation$UniqueID)
save(yall, file = "bismark2DGE_object_full.Rdata")
}
yall <- yall[yall$genes$Chr %in% c(1:21, "X", "Y"), ]
# Select sites in promotor regions
TSS <- nearestTSS(yall$genes$Chr, yall$genes$Locus, species="Hs")
yall$genes$EntrezID <- TSS$gene_id
yall$genes$Symbol <- TSS$symbol
yall$genes$Strand <- TSS$strand
yall$genes$Distance <- TSS$distance
yall$genes$Width <- TSS$width
InPromoter <- yall$genes$Distance >= -1000 & yall$genes$Distance <= 2000
yIP <- yall[InPromoter,,keep.lib.sizes=FALSE]
yall <- rowsum(yIP, yIP$genes$EntrezID, reorder=FALSE)
yall$genes$EntrezID <- NULL
# Normalisation
Methylation <- gl(2,1,ncol(yall), labels=c("Me","Un"))
Me <- yall$counts[, Methylation=="Me"]
Un <- yall$counts[, Methylation=="Un"]
Coverage <- Me + Un
HasCoverage <- rowSums(Coverage >= 15) == 45
HasBoth <- rowSums(Me) > 0 & rowSums(Un) > 0
y <- yall[HasCoverage & HasBoth,, keep.lib.sizes=FALSE]
# Correct for library size
TotalLibSizepr <- 0.5*y$samples$lib.size[Methylation=="Me"] + 0.5*y$samples$lib.size[Methylation=="Un"]
y$samples$lib.size <- rep(TotalLibSizepr, each=2)
# design matrix
designSL <- model.matrix(~0+TimeLapse:Tube, data=sample_annotation)
design <- modelMatrixMeth(designSL)
colnames(design) <- str_replace(colnames(design), ":", "_")
colnames(design) <- str_replace(colnames(design), " ", "_")
y <- estimateDisp(y, design=design, trend="none")
fit <- glmFit(y, design)
}
In the first differential methylation analysis, we compare each time every tube, at timepoint 0.
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
contr <- makeContrasts(contr = TimeLapseT0_TubeBiomatrica - TimeLapseT0_TubePaxgene , levels=design)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
contr <- makeContrasts(contr = TimeLapseT0_TubeDNA_Streck - TimeLapseT0_TubePaxgene , levels=design)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
In the second differential methylation analysis, we compare within tubes: the comparison is now between timepoint 0 and timepoint 72.
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
# DNAStreck T0 vs T72
contr <- makeContrasts(contr = TimeLapseT72_TubeDNA_Streck - TimeLapseT0_TubeDNA_Streck , levels=design)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
# EDTA T0 vs T72
contr <- makeContrasts(contr = TimeLapseT72_TubeEDTA - TimeLapseT0_TubeEDTA , levels=design)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
require(GO.db)
require(org.Hs.eg.db)
cyt.go <- c("GO:0006915", "GO:0071887", "GO:0030101")
term <- select(GO.db, keys=cyt.go, columns="TERM")
datatable(term)
tmp <- org.Hs.egGO2ALLEGS
Rkeys(tmp) <- cyt.go
cyt.go.genes <- as.list(tmp)
out <- fry(y, index=cyt.go.genes, design=design, contrast=contr)
datatable(as.data.frame(out))
res <- glmLRT(fit, contrast=contr)
index <- rownames(fit) %in% cyt.go.genes[[1]]
barcodeplot(res$table$logFC, index=index, labels=c("DOWN","UP"), main=cyt.go[1])
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
datatable(as.data.frame(summary(decideTests(lrt)), col.names = c("direction", "condition", "freq")))
save_csv = data.frame(condition = rownames(contr), val = contr) %>% filter(contr != 0)
save_csv$condition <- str_replace(save_csv$condition, "TimeLapse", "")
save_csv$condition <- str_replace(save_csv$condition, "Tube", "")
save_csv <- save_csv %>% separate(col = condition, into = c("TimeLapse", "Tube"), extra = "merge")
char1 <- save_csv$TimeLapse[1]
char2 <- save_csv$Tube[1]
char3 <- save_csv$TimeLapse[2]
char4 <- save_csv$Tube[2]
write.csv(x = as.data.frame(topTags(lrt, n = 100)) %>% mutate(comparison = paste0(char1, "_", char2, "_", char3, "_", char4)), file = paste0("./tables/", char1, "_", char2, "_", char3, "_", char4, "_diffMeth.csv"), row.names=TRUE)
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] org.Hs.eg.db_3.10.0 GO.db_3.10.0 AnnotationDbi_1.48.0
## [4] IRanges_2.20.2 S4Vectors_0.24.3 Biobase_2.46.0
## [7] BiocGenerics_0.32.0 DT_0.13 RColorBrewer_1.1-2
## [10] readxl_1.3.1 edgeR_3.28.1 limma_3.42.2
## [13] forcats_0.5.0 stringr_1.4.0 dplyr_0.8.5
## [16] purrr_0.3.3 readr_1.3.1 tidyr_1.0.2
## [19] tibble_3.0.0 ggplot2_3.3.0 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.1 bit64_0.9-7 jsonlite_1.6.1 modelr_0.1.6
## [5] assertthat_0.2.1 blob_1.2.1 cellranger_1.1.0 yaml_2.2.1
## [9] pillar_1.4.3 RSQLite_2.2.0 backports_1.1.5 lattice_0.20-38
## [13] glue_1.3.2 digest_0.6.25 rvest_0.3.5 colorspace_1.4-1
## [17] htmltools_0.4.0 pkgconfig_2.0.3 broom_0.5.5 haven_2.2.0
## [21] scales_1.1.0 generics_0.0.2 ellipsis_0.3.0 withr_2.1.2
## [25] cli_2.0.2 magrittr_1.5 crayon_1.3.4 memoise_1.1.0
## [29] evaluate_0.14 fs_1.4.0 fansi_0.4.1 nlme_3.1-144
## [33] xml2_1.3.0 tools_3.6.3 hms_0.5.3 lifecycle_0.2.0
## [37] munsell_0.5.0 reprex_0.3.0 locfit_1.5-9.4 compiler_3.6.3
## [41] rlang_0.4.5 grid_3.6.3 rstudioapi_0.11 htmlwidgets_1.5.1
## [45] crosstalk_1.1.0.1 rmarkdown_2.1 gtable_0.3.0 codetools_0.2-16
## [49] DBI_1.1.0 R6_2.4.1 lubridate_1.7.4 knitr_1.28
## [53] bit_1.1-15.2 stringi_1.4.6 Rcpp_1.0.4 vctrs_0.2.4
## [57] dbplyr_1.4.2 tidyselect_1.0.0 xfun_0.12