Exploratory Data Analysis of Hsu et al
Gaja Matassa
Date: April 30, 2025
for (i in 1:length(params))
print(paste('Parameter:', names(params)[i], ' - Value:', params[[i]], '- Class:', class(params[[i]])))
## [1] "Parameter: InputFolder_WGBS - Value: ~/DataDir/7.EscapeesExploration/WGBS_Hsu/Output/ - Class: character"
## [1] "Parameter: InputFolder_bACEseq - Value: ~/DataDir/7.EscapeesExploration/bACEseq_Hsu/Output/ - Class: character"
## [1] "Parameter: OutputFolder - Value: - Class: character"InputFolder_WGBS <- params$InputFolder_WGBS
InputFolder_bACEseq <- params$InputFolder_bACEseq
# if (dir.exists(OutputFolder) == FALSE) {
# dir.create(OutputFolder, recursive=TRUE)
# }From Hsu et al. paper: A major function of TET1 is the oxidation of 5mC to 5hmC,39 where it can serve as a relatively stable epigenetic mark in non-dividing cells, or as an intermediate for active DNA demethylation. To examine 5hmC enrichment in hPGCLCs, we profiled genome-wide 5hmC levels at single-base resolution using bisulfite-assisted ACEseq (bACEseq).29 By pairing bACEseq with whole genome bisulfite sequencing (WGBS), true 5mC levels and sites can be inferred by subtracting bACEseq signals (5hmC only) from WGBS signals (5mC + 5hmC mixed readout), thereby resolving the epigenetic base ambiguity between 5mC and 5hmC. By profiling ∼40 million (M) common CpG sites between UCLA1 and UCLA2 hESCs and hPGCLCs, we identified a doubling of 5hmC levels in hPGCLCs coupled with a modest reduction of 5mC that can be partially explained by the presence of 5hmC (Figure 2A).
Therefore I decide for our downstream exploration to create a bsseq object with genomic loci shared by the two techniques WGBS and bACEseq and shared by UCLA1 and UCLA2 hESCs and hPGCLCs. It will have assays with methylation levels of 5mCs and 5hmCs separated.
bsseq_obj_WGBS <- readRDS(paste0(InputFolder_WGBS, "bsseq_obj_complete.rds"))
bsseq_obj <- readRDS(paste0(InputFolder_bACEseq, "bsseq_obj_complete.rds"))table(pData(bsseq_obj)$Type)##
## D4hPGCLC hESC hESC_KO
## 6 4 3table(pData(bsseq_obj_WGBS)$Type)##
## D4hPGCLC hESC hESC_KO
## 5 5 3bsseq_obj_WGBS <- methylSig::filter_loci_by_group_coverage(
bs = bsseq_obj_WGBS,
group_column = 'Type',
min_samples_per_group = c('D4hPGCLC' = 5, 'hESC' = 5, 'hESC_KO' = 0))
bsseq_obj <- methylSig::filter_loci_by_group_coverage(
bs = bsseq_obj,
group_column = 'Type',
min_samples_per_group = c('D4hPGCLC' = 6, 'hESC' = 4, 'hESC_KO' = 0))#uniqueInbACEseq <- GenomicRanges::setdiff(bsseq_obj, bsseq_obj_WGBS)
hits <- findOverlaps(granges(bsseq_obj_WGBS), granges(bsseq_obj), type = "equal")
bsseq_WGBS_common <- bsseq_obj_WGBS[queryHits(hits), ]
bsseq_bACE_common <- bsseq_obj[subjectHits(hits), ]Between the two techniques it doesn’t seem that we have the exact same samples. In particular, for UCLA1 hESCs only one replicate is shared and for UCLA2 hESCs there is one more replicate for WGBS than for bACEseq. Furthermore, for UCLA1 hPGCLCs there is one more replicate for bACEseq than WGBS. How to deal with it? Exclude the additional UCLA1 hPGCLC replicate and make a fake pairment of UCLA1_hESC_rep2 in WGBS and UCLA1_hESC_rep1 in bACEseq? Or just keep only the shared replicate for UCLA1 hESCs?
rownames(pData(bsseq_WGBS_common))## [1] "UCLA1_hESC_rep2" "UCLA1_hESC_rep3" "UCLA1_D4hPGCLC_rep2"
## [4] "UCLA1_D4hPGCLC_rep3" "UCLA2_hESC_rep1" "UCLA2_hESC_rep2"
## [7] "UCLA2_hESC_rep3" "UCLA2_D4hPGCLC_rep1" "UCLA2_D4hPGCLC_rep2"
## [10] "UCLA2_D4hPGCLC_rep3" "UCLA1_hESC_KO_CDKO1102" "UCLA1_hESC_KO_CDKO1105"
## [13] "UCLA1_hESC_KO_CDKO1312"rownames(pData(bsseq_bACE_common))## [1] "UCLA1_hESC_rep1" "UCLA1_hESC_rep3" "UCLA1_D4hPGCLC_rep1"
## [4] "UCLA1_D4hPGCLC_rep2" "UCLA1_D4hPGCLC_rep3" "UCLA2_hESC_rep2"
## [7] "UCLA2_hESC_rep3" "UCLA2_D4hPGCLC_rep1" "UCLA2_D4hPGCLC_rep2"
## [10] "UCLA2_D4hPGCLC_rep3" "UCLA1_hESC_KO_CDKO1102" "UCLA1_hESC_KO_CDKO1105"
## [13] "UCLA1_hESC_KO_CDKO1312"samples_to_remove <- c("UCLA1_hESC_rep1", "UCLA1_hESC_rep2", "UCLA2_hESC_rep1", "UCLA1_D4hPGCLC_rep1", "UCLA1_hESC_KO_CDKO1102", "UCLA1_hESC_KO_CDKO1105", "UCLA1_hESC_KO_CDKO1312")
bsseq_WGBS <- bsseq_WGBS_common[ , !rownames(pData(bsseq_WGBS_common)) %in% samples_to_remove]
rm(bsseq_WGBS_common)
bsseq_bACE <- bsseq_bACE_common[ , !rownames(pData(bsseq_bACE_common)) %in% samples_to_remove]
rm(bsseq_bACE_common)rownames(pData(bsseq_WGBS))## [1] "UCLA1_hESC_rep3" "UCLA1_D4hPGCLC_rep2" "UCLA1_D4hPGCLC_rep3"
## [4] "UCLA2_hESC_rep2" "UCLA2_hESC_rep3" "UCLA2_D4hPGCLC_rep1"
## [7] "UCLA2_D4hPGCLC_rep2" "UCLA2_D4hPGCLC_rep3"rownames(pData(bsseq_bACE))## [1] "UCLA1_hESC_rep3" "UCLA1_D4hPGCLC_rep2" "UCLA1_D4hPGCLC_rep3"
## [4] "UCLA2_hESC_rep2" "UCLA2_hESC_rep3" "UCLA2_D4hPGCLC_rep1"
## [7] "UCLA2_D4hPGCLC_rep2" "UCLA2_D4hPGCLC_rep3"Keeping only samples that are shared by both the techniques (without considering the 3 KO hESCs) we go from 10 samples to 8. Let’s proceed with them
granges(bsseq_WGBS)## GRanges object with 27815552 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 10468 *
## [2] chr1 10470 *
## [3] chr1 10483 *
## [4] chr1 10488 *
## [5] chr1 10492 *
## ... ... ... ...
## [27815548] chrKI270757.1 70888 *
## [27815549] chrKI270757.1 70920 *
## [27815550] chrKI270757.1 70921 *
## [27815551] chrKI270757.1 70981 *
## [27815552] chrKI270757.1 70982 *
## -------
## seqinfo: 176 sequences from an unspecified genome; no seqlengthsgranges(bsseq_bACE)## GRanges object with 27815552 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 10468 *
## [2] chr1 10470 *
## [3] chr1 10483 *
## [4] chr1 10488 *
## [5] chr1 10492 *
## ... ... ... ...
## [27815548] chrKI270757.1 70888 *
## [27815549] chrKI270757.1 70920 *
## [27815550] chrKI270757.1 70921 *
## [27815551] chrKI270757.1 70981 *
## [27815552] chrKI270757.1 70982 *
## -------
## seqinfo: 177 sequences from an unspecified genome; no seqlengthsif(!identical(granges(bsseq_WGBS), granges(bsseq_bACE))){
common_seqs <- intersect(seqlevels(bsseq_obj_WGBS), seqlevels(bsseq_obj))
seqlevels(bsseq_WGBS, pruning.mode="coarse") <- common_seqs
seqlevels(bsseq_bACE, pruning.mode="coarse") <- common_seqs}
rm(bsseq_obj_WGBS)
rm(bsseq_obj)if (identical(rownames(pData(bsseq_WGBS)), rownames(pData(bsseq_bACE))) & identical(granges(bsseq_WGBS), granges(bsseq_bACE))) {
bsseq_complete <- bsseq_WGBS
names(assays(bsseq_complete)) <- c("M_WGBS", "Cov_WGBS")
assays(bsseq_complete)$Meth_WGBS <- getMeth(BSseq = bsseq_WGBS, type = "raw")
names(assays(bsseq_complete)) <- c("M_WGBS", "Cov_WGBS", "Meth_WGBS")
assays(bsseq_complete)$M_bACEseq <- getCoverage(bsseq_bACE, type = 'M')
assays(bsseq_complete)$Cov_bACEseq <- getCoverage(bsseq_bACE, type = 'Cov')
assays(bsseq_complete)$Meth_bACEseq <- getMeth(BSseq = bsseq_bACE, type = "raw")
}assays(bsseq_complete)$Meth_hmC <- assays(bsseq_complete)$Meth_bACEseq
assays(bsseq_complete)$Meth_mC <- assays(bsseq_complete)$Meth_WGBS - assays(bsseq_complete)$Meth_bACEseq
summary(assays(bsseq_complete)$Meth_mC)## UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## Min. :-0.7000 Min. :-0.8333 Min. :-0.7500 Min. :-0.7500
## 1st Qu.: 0.6667 1st Qu.: 0.5556 1st Qu.: 0.5444 1st Qu.: 0.6667
## Median : 0.8235 Median : 0.7778 Median : 0.7500 Median : 0.8333
## Mean : 0.7234 Mean : 0.6819 Mean : 0.6579 Mean : 0.7334
## 3rd Qu.: 0.9231 3rd Qu.: 0.9167 3rd Qu.: 0.8750 3rd Qu.: 0.9412
## Max. : 1.0000 Max. : 1.0000 Max. : 1.0000 Max. : 1.0000
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3
## Min. :-0.8000 Min. :-0.9167 Min. :-1.0000 Min. :-0.9091
## 1st Qu.: 0.6807 1st Qu.: 0.5556 1st Qu.: 0.5476 1st Qu.: 0.5561
## Median : 0.8333 Median : 0.7500 Median : 0.7778 Median : 0.7778
## Mean : 0.7302 Mean : 0.6588 Mean : 0.6690 Mean : 0.6708
## 3rd Qu.: 0.9231 3rd Qu.: 0.8750 3rd Qu.: 0.9000 3rd Qu.: 0.9000
## Max. : 1.0000 Max. : 1.0000 Max. : 1.0000 Max. : 1.0000#cap values between 0 and 1 (to avoid small negatives)
assays(bsseq_complete)$Meth_mC[assays(bsseq_complete)$Meth_mC < 0] <- 0
assays(bsseq_complete)$Meth_mC[assays(bsseq_complete)$Meth_mC > 1] <- 1
summary(assays(bsseq_complete)$Meth_mC)## UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.6667 1st Qu.:0.5556 1st Qu.:0.5444 1st Qu.:0.6667
## Median :0.8235 Median :0.7778 Median :0.7500 Median :0.8333
## Mean :0.7241 Mean :0.6831 Mean :0.6588 Mean :0.7341
## 3rd Qu.:0.9231 3rd Qu.:0.9167 3rd Qu.:0.8750 3rd Qu.:0.9412
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.6807 1st Qu.:0.5556 1st Qu.:0.5476 1st Qu.:0.5561
## Median :0.8333 Median :0.7500 Median :0.7778 Median :0.7778
## Mean :0.7308 Mean :0.6598 Mean :0.6705 Mean :0.6718
## 3rd Qu.:0.9231 3rd Qu.:0.8750 3rd Qu.:0.9000 3rd Qu.:0.9000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000Violin plot showing global trend (and mean) of 5mCs + 5hmCs
violin_plot(assays(bsseq_complete)$Meth_WGBS, stat="mean", subset = 3000000)## Subsetting taking random CpGs## No id variables; using all as measure variables
Violin plot showing global trend (and mean) of only 5hmCs
violin_plot(assays(bsseq_complete)$Meth_hmC, stat="mean", subset = 3000000)## Subsetting taking random CpGs## No id variables; using all as measure variables
Violin plot showing global trend (and mean) of only 5mCs
violin_plot(assays(bsseq_complete)$Meth_mC, stat="mean", subset = 3000000)## Subsetting taking random CpGs## No id variables; using all as measure variables
saveRDS(bsseq_complete, "~/DataDir/7.EscapeesExploration/bsseq_Hsu_complete.rds")
saveRDS(bsseq_WGBS, "~/DataDir/7.EscapeesExploration/bsseq_WGBS.rds")
saveRDS(bsseq_bACE, "~/DataDir/7.EscapeesExploration/bsseq_bACE.rds")