Exploratory Data Analysis of Hsu et al WGBS
Gaja Matassa
Date: April 23, 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/WGBS_Hsu - Class: character"
## [1] "Parameter: OutputFolder - Value: ~/DataDir/7.EscapeesExploration/WGBS_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[13]), 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[15]),
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[17]),
sapply(strsplit(basename(bdg_files), split = "-|_")[3:4], function(x) x[16]),
sapply(strsplit(basename(bdg_files), split = "-|_")[5:7], function(x) x[17]),
sapply(strsplit(basename(bdg_files), split = "-|_")[8:10], function(x) x[16]),
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 54594391.
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_rep2 GSM6731108 UCLA1 FALSE hESC
## UCLA1_hESC_rep3 GSM6731109 UCLA1 FALSE hESC
## UCLA1_D4hPGCLC_rep2 GSM6731110 UCLA1 FALSE D4hPGCLC
## UCLA1_D4hPGCLC_rep3 GSM6731111 UCLA1 FALSE D4hPGCLC
## UCLA2_hESC_rep1 GSM6731112 UCLA2 FALSE hESC
## ... ... ... ... ...
## UCLA2_D4hPGCLC_rep2 GSM6731116 UCLA2 FALSE D4hPGCLC
## UCLA2_D4hPGCLC_rep3 GSM6731117 UCLA2 FALSE D4hPGCLC
## UCLA1_hESC_KO_CDKO1102 GSM6731118 UCLA1 TRUE hESC_KO
## UCLA1_hESC_KO_CDKO1105 GSM6731119 UCLA1 TRUE hESC_KO
## UCLA1_hESC_KO_CDKO1312 GSM6731120 UCLA1 TRUE hESC_KO
## Replicate FileName
## <character> <character>
## UCLA1_hESC_rep2 rep2 GSM6731108_FMH-011-E..
## UCLA1_hESC_rep3 rep3 GSM6731109_FMH-011-E..
## UCLA1_D4hPGCLC_rep2 rep2 GSM6731110_FMH-011-E..
## UCLA1_D4hPGCLC_rep3 rep3 GSM6731111_FMH-011-E..
## UCLA2_hESC_rep1 rep1 GSM6731112_FMH-011-E..
## ... ... ...
## UCLA2_D4hPGCLC_rep2 rep2 GSM6731116_FMH-011-E..
## UCLA2_D4hPGCLC_rep3 rep3 GSM6731117_FMH-011-E..
## UCLA1_hESC_KO_CDKO1102 CDKO1102 GSM6731118_UCLA1-hTE..
## UCLA1_hESC_KO_CDKO1105 CDKO1105 GSM6731119_UCLA1-hTE..
## UCLA1_hESC_KO_CDKO1312 CDKO1312 GSM6731120_UCLA1-hTE..
## SampleID
## <character>
## UCLA1_hESC_rep2 UCLA1_hESC_rep2
## UCLA1_hESC_rep3 UCLA1_hESC_rep3
## UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep2
## UCLA1_D4hPGCLC_rep3 UCLA1_D4hPGCLC_rep3
## UCLA2_hESC_rep1 UCLA2_hESC_rep1
## ... ...
## 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: 176 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" "chrKI270315.1"
## [49] "chrKI270316.1" "chrKI270317.1" "chrKI270320.1" "chrKI270322.1"
## [53] "chrKI270330.1" "chrKI270333.1" "chrKI270336.1" "chrKI270337.1"
## [57] "chrKI270362.1" "chrKI270363.1" "chrKI270366.1" "chrKI270378.1"
## [61] "chrKI270382.1" "chrKI270383.1" "chrKI270384.1" "chrKI270385.1"
## [65] "chrKI270386.1" "chrKI270387.1" "chrKI270388.1" "chrKI270389.1"
## [69] "chrKI270390.1" "chrKI270391.1" "chrKI270392.1" "chrKI270393.1"
## [73] "chrKI270394.1" "chrKI270395.1" "chrKI270411.1" "chrKI270414.1"
## [77] "chrKI270417.1" "chrKI270418.1" "chrKI270419.1" "chrKI270420.1"
## [81] "chrKI270422.1" "chrKI270423.1" "chrKI270424.1" "chrKI270425.1"
## [85] "chrKI270429.1" "chrKI270435.1" "chrKI270438.1" "chrKI270442.1"
## [89] "chrKI270448.1" "chrKI270465.1" "chrKI270466.1" "chrKI270467.1"
## [93] "chrKI270468.1" "chrKI270507.1" "chrKI270508.1" "chrKI270509.1"
## [97] "chrKI270510.1" "chrKI270511.1" "chrKI270512.1" "chrKI270515.1"
## [101] "chrKI270516.1" "chrKI270517.1" "chrKI270518.1" "chrKI270519.1"
## [105] "chrKI270521.1" "chrKI270522.1" "chrKI270528.1" "chrKI270529.1"
## [109] "chrKI270530.1" "chrKI270538.1" "chrKI270539.1" "chrKI270544.1"
## [113] "chrKI270579.1" "chrKI270580.1" "chrKI270581.1" "chrKI270582.1"
## [117] "chrKI270583.1" "chrKI270584.1" "chrKI270587.1" "chrKI270588.1"
## [121] "chrKI270589.1" "chrKI270590.1" "chrKI270591.1" "chrKI270593.1"
## [125] "chrKI270706.1" "chrKI270707.1" "chrKI270708.1" "chrKI270709.1"
## [129] "chrKI270710.1" "chrKI270711.1" "chrKI270712.1" "chrKI270713.1"
## [133] "chrKI270714.1" "chrKI270715.1" "chrKI270716.1" "chrKI270717.1"
## [137] "chrKI270718.1" "chrKI270719.1" "chrKI270720.1" "chrKI270721.1"
## [141] "chrKI270722.1" "chrKI270723.1" "chrKI270724.1" "chrKI270725.1"
## [145] "chrKI270726.1" "chrKI270727.1" "chrKI270728.1" "chrKI270729.1"
## [149] "chrKI270730.1" "chrKI270731.1" "chrKI270732.1" "chrKI270733.1"
## [153] "chrKI270734.1" "chrKI270735.1" "chrKI270736.1" "chrKI270737.1"
## [157] "chrKI270738.1" "chrKI270739.1" "chrKI270740.1" "chrKI270741.1"
## [161] "chrKI270742.1" "chrKI270743.1" "chrKI270744.1" "chrKI270745.1"
## [165] "chrKI270746.1" "chrKI270747.1" "chrKI270748.1" "chrKI270749.1"
## [169] "chrKI270750.1" "chrKI270751.1" "chrKI270752.1" "chrKI270753.1"
## [173] "chrKI270754.1" "chrKI270755.1" "chrKI270756.1" "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_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3
## [1,] 14 8 5 11
## [2,] 0 0 0 0
## [3,] 11 8 5 12
## [4,] 0 0 0 0
## [5,] 13 7 5 14
## [6,] 0 0 0 0
## UCLA2_hESC_rep1 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 7 13 12 37
## [2,] 0 0 0 0
## [3,] 7 13 13 35
## [4,] 0 0 0 0
## [5,] 7 13 11 40
## [6,] 0 0 0 0
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 16 27 7
## [2,] 0 0 0
## [3,] 17 28 7
## [4,] 0 0 0
## [5,] 18 28 7
## [6,] 0 0 4
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 5 5
## [2,] 0 5
## [3,] 5 4
## [4,] 0 5
## [5,] 5 0
## [6,] 0 5- 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_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3
## [1,] 16 8 6 15
## [2,] 0 0 0 0
## [3,] 12 8 5 13
## [4,] 0 0 0 0
## [5,] 16 9 6 14
## [6,] 0 0 0 0
## UCLA2_hESC_rep1 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 7 13 13 39
## [2,] 0 0 0 0
## [3,] 7 13 13 38
## [4,] 0 0 0 0
## [5,] 7 13 13 41
## [6,] 0 0 0 0
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 21 27 8
## [2,] 0 0 0
## [3,] 21 28 8
## [4,] 0 0 0
## [5,] 20 28 7
## [6,] 0 0 4
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 5 5
## [2,] 0 5
## [3,] 5 5
## [4,] 0 5
## [5,] 5 0
## [6,] 0 50/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_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3
## [1,] 87.50000 100.00000 83.33333 73.33333
## [2,] NaN NaN NaN NaN
## [3,] 91.66667 100.00000 100.00000 92.30769
## [4,] NaN NaN NaN NaN
## [5,] 81.25000 77.77778 83.33333 100.00000
## [6,] NaN NaN NaN NaN
## UCLA2_hESC_rep1 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 100 100 92.30769 94.87179
## [2,] NaN NaN NaN NaN
## [3,] 100 100 100.00000 92.10526
## [4,] NaN NaN NaN NaN
## [5,] 100 100 84.61538 97.56098
## [6,] NaN NaN NaN NaN
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 76.19048 100 87.5
## [2,] NaN NaN NaN
## [3,] 80.95238 100 87.5
## [4,] NaN NaN NaN
## [5,] 90.00000 100 100.0
## [6,] NaN NaN 100.0
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 100 100
## [2,] NaN 100
## [3,] 100 80
## [4,] NaN 100
## [5,] 100 NaN
## [6,] NaN 100This matrix can also be obtained like this:
head(getMeth(BSseq = bsseq_obj, type = "raw"))## UCLA1_hESC_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3
## [1,] 0.8750000 1.0000000 0.8333333 0.7333333
## [2,] NaN NaN NaN NaN
## [3,] 0.9166667 1.0000000 1.0000000 0.9230769
## [4,] NaN NaN NaN NaN
## [5,] 0.8125000 0.7777778 0.8333333 1.0000000
## [6,] NaN NaN NaN NaN
## UCLA2_hESC_rep1 UCLA2_hESC_rep2 UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1
## [1,] 1 1 0.9230769 0.9487179
## [2,] NaN NaN NaN NaN
## [3,] 1 1 1.0000000 0.9210526
## [4,] NaN NaN NaN NaN
## [5,] 1 1 0.8461538 0.9756098
## [6,] NaN NaN NaN NaN
## UCLA2_D4hPGCLC_rep2 UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102
## [1,] 0.7619048 1 0.875
## [2,] NaN NaN NaN
## [3,] 0.8095238 1 0.875
## [4,] NaN NaN NaN
## [5,] 0.9000000 1 1.000
## [6,] NaN NaN 1.000
## UCLA1_hESC_KO_CDKO1105 UCLA1_hESC_KO_CDKO1312
## [1,] 1 1.0
## [2,] NaN 1.0
## [3,] 1 0.8
## [4,] NaN 1.0
## [5,] 1 NaN
## [6,] NaN 1.0Check 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 54594391.
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_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2
## 12.056125 8.232457 4.854130
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep1 UCLA2_hESC_rep2
## 8.670085 14.733661 7.227999
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 10.782898 12.312763 8.124827
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 8.428990 8.455229 7.465220
## UCLA1_hESC_KO_CDKO1312
## 8.424104Median values for Coverage
apply(Cov_matrix, 2, median, na.rm = TRUE)## UCLA1_hESC_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2
## 12 8 5
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep1 UCLA2_hESC_rep2
## 9 14 7
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 11 12 8
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 8 8 7
## UCLA1_hESC_KO_CDKO1312
## 8Mean values for Methylation Proportions
apply(MethPerc_matrix, 2, mean, na.rm = TRUE)## UCLA1_hESC_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2
## 76.98089 76.56637 73.39489
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep1 UCLA2_hESC_rep2
## 71.21324 78.60342 77.70786
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 78.03795 71.87689 73.10728
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 73.11973 77.41649 76.38198
## UCLA1_hESC_KO_CDKO1312
## 78.66226Median values for Methylation Proportions
apply(MethPerc_matrix, 2, median, na.rm = TRUE)## UCLA1_hESC_rep2 UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2
## 85.71429 85.71429 81.81818
## UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep1 UCLA2_hESC_rep2
## 80.00000 88.23529 87.50000
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## 87.50000 80.00000 81.81818
## UCLA2_D4hPGCLC_rep3 UCLA1_hESC_KO_CDKO1102 UCLA1_hESC_KO_CDKO1105
## 81.81818 85.71429 85.71429
## UCLA1_hESC_KO_CDKO1312
## 87.500003.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' = 5, 'hESC' = 5, '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' = 3, 'hESC' = 3, 'hESC_KO' = 2)), file = paste0 (OutputFolder, "bsseq_obj_sharedby60CellType.rds"))Session Info
SessionInfo <- sessionInfo()
Date <- date()Date## [1] "Wed Apr 23 17:28:59 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