Exploratory Data Analysis of Hsu et al bACEseq
Gaja Matassa
Date: April 29, 2025
for (i in 1:length(params))
print(paste('Parameter:', names(params)[i], ' - Value:', params[[i]], '- Class:', class(params[[i]])))
## [1] "Parameter: InputFolder - Value: ~/DataDir/7.EscapeesExploration/bACEseq_Hsu - Class: character"
## [1] "Parameter: OutputFolder - Value: ~/DataDir/7.EscapeesExploration/bACEseq_Hsu/Output/ - Class: character"InputFolder <- params$InputFolder
OutputFolder <- params$OutputFolder
if (dir.exists(OutputFolder) == FALSE) {
dir.create(OutputFolder, recursive=TRUE)
}Colors for figures
# sample_colors <- c("3_EG2_SL" = "#ff7f00", "3_EG9_SL" = "#ff7f00", "3_EG31_SL" = "#ff7f00",
# "2S" = "#dd1c77", "3S" = "#dd1c77", "5S" = "#dd1c77", "4S" = "#dd1c77",
# "2M" = "#377eb8", "3M" = "#377eb8", "5M" = "#377eb8", "4M" = "#377eb8",
# "2P" = "#4daf4a", "3P" = "#4daf4a", "5P" = "#4daf4a")
#
cellTypes_colors <- c(hESC_KO ="#dd1c77", D4hPGCLC ="#4daf4a", hESC = "#ff7f00")1. Loading data
BedGraph files in input are already merged (from individual Cytosines in a CpG to per-CpG metrics) and filtered (threshold: 10 reads). Here bedIntersect was already used to keep only on-bait CpGs.
bdg_files <- list.files(path = InputFolder, pattern = '*sorted.bedGraph', full.names = TRUE) #output from our preprocessingCreation of metadata from file names
sample_anno <- data.frame("SampleGEO" = sapply(strsplit(basename(bdg_files), split = "-|_"), function(x) x[1]),
"Line" = c(sapply(strsplit(basename(bdg_files), split = "-|_")[1:10], function(x) x[11]), sapply(strsplit(basename(bdg_files), split = "-|_")[11:13], function(x) x[2])),
"TET1KO" = c(rep(FALSE,10), rep(TRUE,3)),
"Type" = c(sapply(strsplit(basename(bdg_files), split = "-|_")[1:10], function(x) x[13]),
paste0(sapply(strsplit(basename(bdg_files), split = "-|_")[11:13], function(x) x[5]), "_KO")),
"Replicate" = c(sapply(strsplit(basename(bdg_files), split = "-|_")[1:2], function(x) x[15]),
sapply(strsplit(basename(bdg_files), split = "-|_")[3:5], function(x) x[14]),
sapply(strsplit(basename(bdg_files), split = "-|_")[6:7], function(x) x[15]),
sapply(strsplit(basename(bdg_files), split = "-|_")[8:10], function(x) x[14]),
sapply(strsplit(basename(bdg_files), split = "-|_")[11:13], function(x) x[4])),
#"Sequencing" = rep("WGBS", length(bdg_files)),
#"Day" = c(sapply(strsplit(basename(bdg_files), split = "-|_")[1:10], function(x) x[14]), NA, NA, NA)
"FileName" = basename(bdg_files))
sample_anno$SampleID <- paste(sample_anno$Line, sample_anno$Type, sample_anno$Replicate, sep = "_")
if(!any(duplicated(sample_anno$SampleID))){
rownames(sample_anno) <- sample_anno$SampleID
}else{stop("Duplicated sample names!")}Match order
bdg_files<-bdg_files[order(match(basename(bdg_files), sample_anno$FileName))]2. Reading bedgraph files with read.bismark from
bsseq package
2.1 Generation of BSseq object
strandCollapse = FALSE because MethylDackel has an
option to destrand data when methylation calls are made so that the
output is already destranded.
if(identical(basename(bdg_files), sample_anno$FileName)){
bsseq_obj = bsseq::read.bismark(
files = bdg_files,
colData = sample_anno,
rmZeroCov = FALSE,
strandCollapse = FALSE #false if MethylDackel is used because MethylDackel has an option to destrand data when methylation calls are made so that the output is already destranded.
)
} else{stop("Sample names don't match between Bedgraph file and metadata")}- The number of samples in the BSSeq object is 13.
- The number of CpGs in the BSSeq object is 54853237.
2.2 Exploration of BSseq object
bsseq::pData(bsseq_obj)## DataFrame with 13 rows and 7 columns
## SampleGEO Line TET1KO Type
## <character> <character> <logical> <character>
## UCLA1_hESC_rep1 GSM6759276 UCLA1 FALSE hESC
## UCLA1_hESC_rep3 GSM6759277 UCLA1 FALSE hESC
## UCLA1_D4hPGCLC_rep1 GSM6759278 UCLA1 FALSE D4hPGCLC
## UCLA1_D4hPGCLC_rep2 GSM6759279 UCLA1 FALSE D4hPGCLC
## UCLA1_D4hPGCLC_rep3 GSM6759280 UCLA1 FALSE D4hPGCLC
## ... ... ... ... ...
## UCLA2_D4hPGCLC_rep2 GSM6759284 UCLA2 FALSE D4hPGCLC
## UCLA2_D4hPGCLC_rep3 GSM6759285 UCLA2 FALSE D4hPGCLC
## UCLA1_hESC_KO_CDKO1102 GSM6759286 UCLA1 TRUE hESC_KO
## UCLA1_hESC_KO_CDKO1105 GSM6759287 UCLA1 TRUE hESC_KO
## UCLA1_hESC_KO_CDKO1312 GSM6759288 UCLA1 TRUE hESC_KO
## Replicate FileName
## <character> <character>
## UCLA1_hESC_rep1 rep1 GSM6759276_FMH-011-E..
## UCLA1_hESC_rep3 rep3 GSM6759277_FMH-011-E..
## UCLA1_D4hPGCLC_rep1 rep1 GSM6759278_FMH-011-E..
## UCLA1_D4hPGCLC_rep2 rep2 GSM6759279_FMH-011-E..
## UCLA1_D4hPGCLC_rep3 rep3 GSM6759280_FMH-011-E..
## ... ... ...
## UCLA2_D4hPGCLC_rep2 rep2 GSM6759284_FMH-011-E..
## UCLA2_D4hPGCLC_rep3 rep3 GSM6759285_FMH-011-E..
## UCLA1_hESC_KO_CDKO1102 CDKO1102 GSM6759286_UCLA1-hTE..
## UCLA1_hESC_KO_CDKO1105 CDKO1105 GSM6759287_UCLA1-hTE..
## UCLA1_hESC_KO_CDKO1312 CDKO1312 GSM6759288_UCLA1-hTE..
## SampleID
## <character>
## UCLA1_hESC_rep1 UCLA1_hESC_rep1
## UCLA1_hESC_rep3 UCLA1_hESC_rep3
## UCLA1_D4hPGCLC_rep1 UCLA1_D4hPGCLC_rep1
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep2
## UCLA1_D4hPGCLC_rep3 UCLA1_D4hPGCLC_rep3
## ... ...
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep2
## UCLA2_D4hPGCLC_rep3 UCLA2_D4hPGCLC_rep3
## UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1102
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1105
## UCLA1_hESC_KO_CDKO1312 UCLA1_hESC_KO_CDKO1312bsseq package assumes that the following data has been
extracted from alignments:
- genomic positions, including chromosome and location, for
methylation loci.The genomic positions are stored in the bsseq object as
a
GRangesobject.GRangesare general genomic regions; they represent a single base methylation loci as an interval of width 1 (which may seem a bit strange, but they said there are good reasons for this).
head(GenomicRanges::granges(bsseq_obj))## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 10468 *
## [2] chr1 10469 *
## [3] chr1 10470 *
## [4] chr1 10471 *
## [5] chr1 10483 *
## [6] chr1 10484 *
## -------
## seqinfo: 177 sequences from an unspecified genome; no seqlengthsMethylation calls are from these chromosomes:
levels(bsseq_obj@rowRanges@seqnames@values)## [1] "chr1" "chr2" "chr3" "chr4"
## [5] "chr5" "chr6" "chr7" "chr8"
## [9] "chr9" "chr10" "chr11" "chr12"
## [13] "chr13" "chr14" "chr15" "chr16"
## [17] "chr17" "chr18" "chr19" "chr20"
## [21] "chr21" "chr22" "chrX" "chrY"
## [25] "chrMT" "chrGL000008.2" "chrGL000009.2" "chrGL000194.1"
## [29] "chrGL000195.1" "chrGL000205.2" "chrGL000208.1" "chrGL000213.1"
## [33] "chrGL000214.1" "chrGL000216.2" "chrGL000218.1" "chrGL000219.1"
## [37] "chrGL000220.1" "chrGL000221.1" "chrGL000224.1" "chrGL000225.1"
## [41] "chrGL000226.1" "chrKI270302.1" "chrKI270303.1" "chrKI270304.1"
## [45] "chrKI270305.1" "chrKI270310.1" "chrKI270311.1" "chrKI270312.1"
## [49] "chrKI270315.1" "chrKI270316.1" "chrKI270317.1" "chrKI270320.1"
## [53] "chrKI270322.1" "chrKI270330.1" "chrKI270333.1" "chrKI270336.1"
## [57] "chrKI270337.1" "chrKI270362.1" "chrKI270363.1" "chrKI270366.1"
## [61] "chrKI270378.1" "chrKI270382.1" "chrKI270383.1" "chrKI270384.1"
## [65] "chrKI270385.1" "chrKI270386.1" "chrKI270387.1" "chrKI270388.1"
## [69] "chrKI270389.1" "chrKI270390.1" "chrKI270391.1" "chrKI270392.1"
## [73] "chrKI270393.1" "chrKI270394.1" "chrKI270395.1" "chrKI270411.1"
## [77] "chrKI270412.1" "chrKI270414.1" "chrKI270417.1" "chrKI270418.1"
## [81] "chrKI270419.1" "chrKI270420.1" "chrKI270422.1" "chrKI270423.1"
## [85] "chrKI270424.1" "chrKI270425.1" "chrKI270429.1" "chrKI270435.1"
## [89] "chrKI270438.1" "chrKI270442.1" "chrKI270448.1" "chrKI270465.1"
## [93] "chrKI270466.1" "chrKI270467.1" "chrKI270468.1" "chrKI270507.1"
## [97] "chrKI270508.1" "chrKI270509.1" "chrKI270510.1" "chrKI270511.1"
## [101] "chrKI270512.1" "chrKI270515.1" "chrKI270516.1" "chrKI270517.1"
## [105] "chrKI270518.1" "chrKI270519.1" "chrKI270521.1" "chrKI270522.1"
## [109] "chrKI270528.1" "chrKI270529.1" "chrKI270530.1" "chrKI270538.1"
## [113] "chrKI270539.1" "chrKI270544.1" "chrKI270579.1" "chrKI270580.1"
## [117] "chrKI270582.1" "chrKI270583.1" "chrKI270584.1" "chrKI270587.1"
## [121] "chrKI270588.1" "chrKI270589.1" "chrKI270590.1" "chrKI270591.1"
## [125] "chrKI270593.1" "chrKI270706.1" "chrKI270707.1" "chrKI270708.1"
## [129] "chrKI270709.1" "chrKI270710.1" "chrKI270711.1" "chrKI270712.1"
## [133] "chrKI270713.1" "chrKI270714.1" "chrKI270715.1" "chrKI270716.1"
## [137] "chrKI270717.1" "chrKI270718.1" "chrKI270719.1" "chrKI270720.1"
## [141] "chrKI270721.1" "chrKI270722.1" "chrKI270723.1" "chrKI270724.1"
## [145] "chrKI270725.1" "chrKI270726.1" "chrKI270727.1" "chrKI270728.1"
## [149] "chrKI270729.1" "chrKI270730.1" "chrKI270731.1" "chrKI270732.1"
## [153] "chrKI270733.1" "chrKI270734.1" "chrKI270735.1" "chrKI270736.1"
## [157] "chrKI270737.1" "chrKI270738.1" "chrKI270739.1" "chrKI270740.1"
## [161] "chrKI270741.1" "chrKI270742.1" "chrKI270743.1" "chrKI270744.1"
## [165] "chrKI270745.1" "chrKI270746.1" "chrKI270747.1" "chrKI270748.1"
## [169] "chrKI270749.1" "chrKI270750.1" "chrKI270751.1" "chrKI270752.1"
## [173] "chrKI270753.1" "chrKI270754.1" "chrKI270755.1" "chrKI270756.1"
## [177] "chrKI270757.1"- a (matrix) of M (Methylation) values, describing the number of read supporting methylation covering a single loci. Each row in this matrix is a methylation loci and each column is a sample.
Meth_matrix <- bsseq::getCoverage(bsseq_obj, type = 'M')
head(Meth_matrix)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1 UCLA1_D4hPGCLC_rep2
## [1,] 3 1 3 0
## [2,] 0 0 0 0
## [3,] 1 1 1 0
## [4,] 0 0 0 0
## [5,] 1 2 1 0
## [6,] 0 0 0 0
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 0 0 0 0
## [2,] 0 0 0 0
## [3,] 1 0 1 0
## [4,] 0 0 0 0
## [5,] 1 1 1 0
## [6,] 0 0 0 0
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 1 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 0 0
## [2,] 0 0
## [3,] 0 0
## [4,] 0 0
## [5,] 0 0
## [6,] 0 0- a (matrix) of Cov (Coverage) values, describing the total number of reads covering a single loci. Each row in this matrix is a methylation loci and each column is a sample.
Cov_matrix <- bsseq::getCoverage(bsseq_obj, type = 'Cov')
head(Cov_matrix)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1 UCLA1_D4hPGCLC_rep2
## [1,] 42 26 41 15
## [2,] 6 4 8 0
## [3,] 43 26 40 15
## [4,] 5 4 6 0
## [5,] 41 30 45 16
## [6,] 5 0 5 0
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 15 20 17 18
## [2,] 0 0 0 0
## [3,] 14 20 17 18
## [4,] 0 0 0 0
## [5,] 18 22 15 19
## [6,] 0 0 0 0
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 10 19 5
## [2,] 0 0 0
## [3,] 10 21 5
## [4,] 0 0 0
## [5,] 10 20 6
## [6,] 0 0 0
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 0 6
## [2,] 0 4
## [3,] 0 6
## [4,] 0 4
## [5,] 0 8
## [6,] 0 00/0 gives NaN, this is the case when 0 methylated reads mapping on that region over 0 total reads mapping on that region.
Methylation Proportions/Levels = Meth_matrix/Cov_matrix
MethPerc_matrix <- (Meth_matrix/Cov_matrix)*100
head(MethPerc_matrix)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1 UCLA1_D4hPGCLC_rep2
## [1,] 7.142857 3.846154 7.317073 0
## [2,] 0.000000 0.000000 0.000000 NaN
## [3,] 2.325581 3.846154 2.500000 0
## [4,] 0.000000 0.000000 0.000000 NaN
## [5,] 2.439024 6.666667 2.222222 0
## [6,] 0.000000 NaN 0.000000 NaN
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 0.000000 0.000000 0.000000 0
## [2,] NaN NaN NaN NaN
## [3,] 7.142857 0.000000 5.882353 0
## [4,] NaN NaN NaN NaN
## [5,] 5.555556 4.545455 6.666667 0
## [6,] NaN NaN NaN NaN
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 0 0.000000 0
## [2,] NaN NaN NaN
## [3,] 0 4.761905 0
## [4,] NaN NaN NaN
## [5,] 0 0.000000 0
## [6,] NaN NaN NaN
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] NaN 0
## [2,] NaN 0
## [3,] NaN 0
## [4,] NaN 0
## [5,] NaN 0
## [6,] NaN NaNThis matrix can also be obtained like this:
head(getMeth(BSseq = bsseq_obj, type = "raw"))## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1 UCLA1_D4hPGCLC_rep2
## [1,] 0.07142857 0.03846154 0.07317073 0
## [2,] 0.00000000 0.00000000 0.00000000 NaN
## [3,] 0.02325581 0.03846154 0.02500000 0
## [4,] 0.00000000 0.00000000 0.00000000 NaN
## [5,] 0.02439024 0.06666667 0.02222222 0
## [6,] 0.00000000 NaN 0.00000000 NaN
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 0.00000000 0.00000000 0.00000000 0
## [2,] NaN NaN NaN NaN
## [3,] 0.07142857 0.00000000 0.05882353 0
## [4,] NaN NaN NaN NaN
## [5,] 0.05555556 0.04545455 0.06666667 0
## [6,] NaN NaN NaN NaN
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 0 0.00000000 0
## [2,] NaN NaN NaN
## [3,] 0 0.04761905 0
## [4,] NaN NaN NaN
## [5,] 0 0.00000000 0
## [6,] NaN NaN NaN
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] NaN 0
## [2,] NaN 0
## [3,] NaN 0
## [4,] NaN 0
## [5,] NaN 0
## [6,] NaN NaNCheck whether there are regions with NaNs for all samples
Question: Is it true that there are no regions uncovered by any samples?
all(apply(X = is.na(MethPerc_matrix), MARGIN = 1, FUN = sum) < ncol(MethPerc_matrix))## [1] TRUENumber of CpGs
sample_anno$NumberCpGsPerSample <- colSums(Cov_matrix!=0)
sample_anno$NumberCpGsinBSseq <- rep(nrow(bsseq_obj), ncol(bsseq_obj))
sample_anno$PercCpGPerSampleinBSeq <- round((1-(colSums(Cov_matrix==0)/nrow(bsseq_obj)))*100, digits = 2)
datatable(sample_anno, extensions = "Buttons",
options = list(paging = TRUE,
scrollX=TRUE,
searching = TRUE,
ordering = TRUE,
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf'),
pageLength=10,
lengthMenu=c(3,5,10) ))3. Exploratory plots
The number of CpGs covered by at least one sample in the BSSeq object is 54853237.
3.3 Violin Plot showing mean
violin_plot(MethPerc_matrix, stat="mean", subset = 3000000) #+ ggplot2::scale_fill_manual(values = sample_colors)
## Subsetting taking random CpGs
## No id variables; using all as measure variables
3.4 Violin Plot showing median
violin_plot(MethPerc_matrix, subset = 3000000) #+ ggplot2::scale_fill_manual(values = sample_colors)
## Subsetting taking random CpGs
## No id variables; using all as measure variables
Mean values for Coverage
apply(Cov_matrix, 2, mean, na.rm = TRUE)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1
## 18.949867 14.656110 17.058471
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## 10.312723 11.888623 12.704628
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 10.051916 9.910634 7.491432
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 10.262596 6.193677 6.849366
## UCLA1_hESC_KO_CDKO1312
## 7.953816Median values for Coverage
apply(Cov_matrix, 2, median, na.rm = TRUE)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1
## 19 15 17
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## 10 12 13
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 10 10 7
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 10 6 7
## UCLA1_hESC_KO_CDKO1312
## 8Mean values for Methylation Proportions
apply(MethPerc_matrix, 2, mean, na.rm = TRUE)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1
## 2.395617 2.515441 4.120048
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## 3.698001 3.540895 2.695498
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 3.022686 3.901569 4.024746
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 3.919115 1.641083 1.529154
## UCLA1_hESC_KO_CDKO1312
## 1.433811Median values for Methylation Proportions
apply(MethPerc_matrix, 2, median, na.rm = TRUE)## UCLA1_hESC_rep1 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep1
## 0 0 0
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## 0 0 0
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 0 0 0
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 0 0 0
## UCLA1_hESC_KO_CDKO1312
## 03.5 PCA plots
pca_res <- do_PCA(MethPerc_matrix)plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = "Line", custom_colors = cellTypes_colors, point_size = 5)
5. Saving and Session Info
Saving BSSeq objects :
- Complete
saveRDS(bsseq_obj, file = paste0 (OutputFolder, "bsseq_obj_complete.rds"))- Filtered by sample and by CpGs covered by all samples kept
pData(bsseq_obj)$Seq <- "WGBS"
saveRDS(methylSig::filter_loci_by_group_coverage(bs = bsseq_obj,group_column = "Seq", min_samples_per_group = c("WGBS" = dim(bsseq_obj)[2])), file = paste0 (OutputFolder, "bsseq_obj_sharedbyall.rds")) #same of doing group_column = 'Type', min_samples_per_group = c('D4hPGCLC' = 6, 'hESC' = 4, 'hESC_KO' = 3))- With CpGs covered by 75% of all samples kept (independently on cell types)
saveRDS(methylSig::filter_loci_by_group_coverage(bs = bsseq_obj,group_column = "Seq", min_samples_per_group = c("WGBS" = round(0.75*dim(bsseq_obj)[2], digits=0))), file = paste0 (OutputFolder, "bsseq_obj_sharedby75ofall.rds"))- With CpGs covered by at least 60% of samples of the same cell type
saveRDS(methylSig::filter_loci_by_group_coverage(bs = bsseq_obj,group_column = "Type", min_samples_per_group = c('D4hPGCLC' = 4, 'hESC' = 3, 'hESC_KO' = 2)), file = paste0 (OutputFolder, "bsseq_obj_sharedby60CellType.rds"))Session Info
SessionInfo <- sessionInfo()
Date <- date()Date## [1] "Tue Apr 29 13:50:08 2025"SessionInfo## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] methylSig_1.10.0 dmrseq_1.18.1
## [3] DT_0.28 data.table_1.14.8
## [5] tidyr_1.3.0 dplyr_1.1.2
## [7] ggrepel_0.9.3 ggplot2_3.4.2
## [9] bsseq_1.34.0 SummarizedExperiment_1.28.0
## [11] Biobase_2.58.0 MatrixGenerics_1.10.0
## [13] matrixStats_1.0.0 GenomicRanges_1.50.2
## [15] GenomeInfoDb_1.34.9 IRanges_2.32.0
## [17] S4Vectors_0.36.2 BiocGenerics_0.44.0
##
## loaded via a namespace (and not attached):
## [1] AnnotationHub_3.6.0 BiocFileCache_2.6.1
## [3] plyr_1.8.8 splines_4.2.1
## [5] crosstalk_1.2.0 BiocParallel_1.32.6
## [7] digest_0.6.33 foreach_1.5.2
## [9] htmltools_0.5.5 fansi_1.0.4
## [11] magrittr_2.0.3 memoise_2.0.1
## [13] BSgenome_1.66.3 tzdb_0.4.0
## [15] limma_3.54.2 Biostrings_2.66.0
## [17] readr_2.1.4 R.utils_2.12.2
## [19] prettyunits_1.1.1 colorspace_2.1-0
## [21] blob_1.2.4 rappdirs_0.3.3
## [23] xfun_0.39 crayon_1.5.2
## [25] RCurl_1.98-1.12 jsonlite_1.8.7
## [27] annotatr_1.24.0 iterators_1.0.14
## [29] glue_1.6.2 gtable_0.3.3
## [31] zlibbioc_1.44.0 XVector_0.38.0
## [33] DelayedArray_0.24.0 Rhdf5lib_1.20.0
## [35] HDF5Array_1.26.0 scales_1.2.1
## [37] DBI_1.1.3 rngtools_1.5.2
## [39] Rcpp_1.0.11 DSS_2.46.0
## [41] xtable_1.8-4 progress_1.2.2
## [43] bumphunter_1.40.0 bit_4.0.5
## [45] htmlwidgets_1.6.2 httr_1.4.6
## [47] RColorBrewer_1.1-3 ellipsis_0.3.2
## [49] farver_2.1.1 pkgconfig_2.0.3
## [51] XML_3.99-0.14 R.methodsS3_1.8.2
## [53] sass_0.4.7 dbplyr_2.3.3
## [55] locfit_1.5-9.7 utf8_1.2.3
## [57] labeling_0.4.2 tidyselect_1.2.0
## [59] rlang_1.1.1 reshape2_1.4.4
## [61] later_1.3.1 AnnotationDbi_1.60.2
## [63] munsell_0.5.0 BiocVersion_3.16.0
## [65] tools_4.2.1 cachem_1.0.8
## [67] cli_3.6.1 generics_0.1.3
## [69] RSQLite_2.3.1 evaluate_0.21
## [71] stringr_1.5.0 fastmap_1.1.1
## [73] yaml_2.3.7 outliers_0.15
## [75] knitr_1.43 bit64_4.0.5
## [77] purrr_1.0.1 KEGGREST_1.38.0
## [79] nlme_3.1-162 doRNG_1.8.6
## [81] sparseMatrixStats_1.10.0 mime_0.12
## [83] R.oo_1.25.0 xml2_1.3.5
## [85] biomaRt_2.54.1 compiler_4.2.1
## [87] rstudioapi_0.15.0 filelock_1.0.2
## [89] curl_5.0.1 png_0.1-8
## [91] interactiveDisplayBase_1.36.0 tibble_3.2.1
## [93] bslib_0.5.0 stringi_1.7.12
## [95] highr_0.10 GenomicFeatures_1.50.4
## [97] lattice_0.21-8 Matrix_1.6-0
## [99] permute_0.9-7 vctrs_0.6.3
## [101] pillar_1.9.0 lifecycle_1.0.3
## [103] rhdf5filters_1.10.1 BiocManager_1.30.20
## [105] jquerylib_0.1.4 bitops_1.0-7
## [107] httpuv_1.6.11 rtracklayer_1.58.0
## [109] R6_2.5.1 BiocIO_1.8.0
## [111] promises_1.2.0.1 codetools_0.2-19
## [113] gtools_3.9.4 rhdf5_2.42.1
## [115] rjson_0.2.21 withr_2.5.0
## [117] regioneR_1.30.0 GenomicAlignments_1.34.1
## [119] Rsamtools_2.14.0 GenomeInfoDbData_1.2.9
## [121] parallel_4.2.1 hms_1.1.3
## [123] grid_4.2.1 rmarkdown_2.23
## [125] DelayedMatrixStats_1.20.0 shiny_1.7.4.1
## [127] restfulr_0.0.15