for (i in 1:length(params))
  print(paste('Parameter:', names(params)[i], ' - Value:', params[[i]], '- Class:', class(params[[i]])))
## [1] "Parameter: InputFolder  - Value: ~/DataDir/0.Preprocessing/Output/ - Class: character"
## [1] "Parameter: MetadataFolder  - Value: ~/DataDir/1.PreliminaryAnalysis/Input/ - Class: character"
## [1] "Parameter: OutputFolder  - Value: ~/DataDir/1.PreliminaryAnalysis/Output/WithSelectedEGCLCs/ - Class: character"

InputFolder <- params$InputFolder
MetadataFolder <- params$MetadataFolder
OutputFolder <- params$OutputFolder

if (dir.exists(OutputFolder) == FALSE) {
  dir.create(OutputFolder, recursive=TRUE)
}

Colors for figures

sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_rep1" = "#ff7f00", "hEGCLC_rep2" = "#ff7f00", "hEGCLC_rep3" = "#ff7f00")

cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs = "#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 = '*\\.bedGraph', full.names = TRUE) #output from our preprocessing

bdg_files_filtered <- bdg_files[grep("filtered", bdg_files)]

sample_anno <- read.table(paste0(MetadataFolder, "/SampleAnno.txt"), sep='\t', header=TRUE)
rownames(sample_anno) <- sample_anno$SampleID

Keeping only bedGraph files of selected samples.

bdg_files_filtered <- bdg_files_filtered[which(grepl(paste(sample_anno$OriginalID, collapse = "|"), bdg_files_filtered) )]

Match order

bdg_files_filtered <- bdg_files_filtered[order(match(sub(".filtered.bedGraph", "", basename(bdg_files_filtered)), sample_anno$OriginalID))]

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(sub(".filtered.bedGraph", "", basename(bdg_files_filtered)), sample_anno$OriginalID)){
  bsseq_obj = bsseq::read.bismark(
      files = bdg_files_filtered,
      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 14.
The number of CpGs in the BSSeq object is 3732776.

2.2 Exploration of BSseq object

bsseq::pData(bsseq_obj)

## DataFrame with 14 rows and 11 columns
##                SampleID  OriginalID        Type       Day     Batch        Line
##             <character> <character> <character> <integer> <integer> <character>
## hiPSC_rep1   hiPSC_rep1          2S      hiPSCs        NA         2      CTL08A
## hiPSC_rep2   hiPSC_rep2          3S      hiPSCs        NA         3      CTL08A
## hiPSC_rep3   hiPSC_rep3          5S      hiPSCs        NA         5      CTL08A
## hiPSC_rep4   hiPSC_rep4          4S      hiPSCs        NA         4      CTL08A
## iMeLC_rep1   iMeLC_rep1          2M      iMeLCs        NA         2      CTL08A
## ...                 ...         ...         ...       ...       ...         ...
## hPGCLC_rep2 hPGCLC_rep2          3P     hPGCLCs         6         3      CTL08A
## hPGCLC_rep3 hPGCLC_rep3          5P     hPGCLCs         6         5      CTL08A
## hEGCLC_rep1 hEGCLC_rep1    3_EG2_SL     hEGCLCs        NA         3      CTL08A
## hEGCLC_rep2 hEGCLC_rep2    3_EG9_SL     hEGCLCs        NA         3      CTL08A
## hEGCLC_rep3 hEGCLC_rep3   3_EG31_SL     hEGCLCs        NA         3      CTL08A
##                 Clone      Medium  PercOxigen PercEfficiency       Notes
##             <integer> <character> <character>    <character> <character>
## hiPSC_rep1         NA                                                   
## hiPSC_rep2         NA                                                   
## hiPSC_rep3         NA                                                   
## hiPSC_rep4         NA                                                   
## iMeLC_rep1         NA                                                   
## ...               ...         ...         ...            ...         ...
## hPGCLC_rep2        NA                                  0_029   low_input
## hPGCLC_rep3        NA                                 0_0152    ok_input
## hEGCLC_rep1         2     StemFit      low_ox                           
## hEGCLC_rep2         9     StemFit      low_ox                           
## hEGCLC_rep3        31     StemFit      low_ox

bsseq 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 GRanges object. GRanges are 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     10470      *
##   [3]     chr1     10483      *
##   [4]     chr1     10488      *
##   [5]     chr1     10492      *
##   [6]     chr1     10496      *
##   -------
##   seqinfo: 25 sequences from an unspecified genome; no seqlengths

Methylation calls are from these chromosomes:

levels(bsseq_obj@rowRanges@seqnames@values)

##  [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9" 
## [10] "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18"
## [19] "chr19" "chr20" "chr21" "chr22" "chrX"  "chrY"  "chrM"

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)

##      hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,]         16         18         19         21         18         28
## [2,]         18         20         22         25         24         32
## [3,]         23         21         23         25         27         39
## [4,]         23         22         25         26         22         40
## [5,]         24         21         26         27         27         40
## [6,]         23         22         27         27         25         42
##      iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,]         27         16          13          11          23          40
## [2,]         28         17          17          12          28          45
## [3,]         30         19          19          16          28          49
## [4,]         30         21          18          16          25          52
## [5,]         30         23          20          15          25          51
## [6,]         30         23          19          17          26          53
##      hEGCLC_rep2 hEGCLC_rep3
## [1,]          31          27
## [2,]          36          32
## [3,]          40          44
## [4,]          43          42
## [5,]          46          42
## [6,]          49          47

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)

##      hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,]         19         21         23         25         24         35
## [2,]         18         20         22         25         25         34
## [3,]         23         22         23         26         28         39
## [4,]         23         23         26         27         27         41
## [5,]         24         22         27         29         28         42
## [6,]         23         23         28         30         26         42
##      iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,]         31         19          20          16          29          45
## [2,]         29         18          17          12          28          46
## [3,]         32         22          22          16          29          54
## [4,]         32         24          21          16          28          55
## [5,]         32         25          24          16          27          55
## [6,]         32         27          23          17          30          56
##      hEGCLC_rep2 hEGCLC_rep3
## [1,]          43          37
## [2,]          37          34
## [3,]          42          46
## [4,]          46          47
## [5,]          48          48
## [6,]          51          49

0/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)

##      hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,]   84.21053   85.71429   82.60870   84.00000   75.00000   80.00000
## [2,]  100.00000  100.00000  100.00000  100.00000   96.00000   94.11765
## [3,]  100.00000   95.45455  100.00000   96.15385   96.42857  100.00000
## [4,]  100.00000   95.65217   96.15385   96.29630   81.48148   97.56098
## [5,]  100.00000   95.45455   96.29630   93.10345   96.42857   95.23810
## [6,]  100.00000   95.65217   96.42857   90.00000   96.15385  100.00000
##      iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,]   87.09677   84.21053    65.00000       68.75    79.31034    88.88889
## [2,]   96.55172   94.44444   100.00000      100.00   100.00000    97.82609
## [3,]   93.75000   86.36364    86.36364      100.00    96.55172    90.74074
## [4,]   93.75000   87.50000    85.71429      100.00    89.28571    94.54545
## [5,]   93.75000   92.00000    83.33333       93.75    92.59259    92.72727
## [6,]   93.75000   85.18519    82.60870      100.00    86.66667    94.64286
##      hEGCLC_rep2 hEGCLC_rep3
## [1,]    72.09302    72.97297
## [2,]    97.29730    94.11765
## [3,]    95.23810    95.65217
## [4,]    93.47826    89.36170
## [5,]    95.83333    87.50000
## [6,]    96.07843    95.91837

This matrix can also be obtained like this:

head(getMeth(BSseq = bsseq_obj, type = "raw"))

##      hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,]  0.8421053  0.8571429  0.8260870  0.8400000  0.7500000  0.8000000
## [2,]  1.0000000  1.0000000  1.0000000  1.0000000  0.9600000  0.9411765
## [3,]  1.0000000  0.9545455  1.0000000  0.9615385  0.9642857  1.0000000
## [4,]  1.0000000  0.9565217  0.9615385  0.9629630  0.8148148  0.9756098
## [5,]  1.0000000  0.9545455  0.9629630  0.9310345  0.9642857  0.9523810
## [6,]  1.0000000  0.9565217  0.9642857  0.9000000  0.9615385  1.0000000
##      iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,]  0.8709677  0.8421053   0.6500000      0.6875   0.7931034   0.8888889
## [2,]  0.9655172  0.9444444   1.0000000      1.0000   1.0000000   0.9782609
## [3,]  0.9375000  0.8636364   0.8636364      1.0000   0.9655172   0.9074074
## [4,]  0.9375000  0.8750000   0.8571429      1.0000   0.8928571   0.9454545
## [5,]  0.9375000  0.9200000   0.8333333      0.9375   0.9259259   0.9272727
## [6,]  0.9375000  0.8518519   0.8260870      1.0000   0.8666667   0.9464286
##      hEGCLC_rep2 hEGCLC_rep3
## [1,]   0.7209302   0.7297297
## [2,]   0.9729730   0.9411765
## [3,]   0.9523810   0.9565217
## [4,]   0.9347826   0.8936170
## [5,]   0.9583333   0.8750000
## [6,]   0.9607843   0.9591837

Check 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] TRUE

Number 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 3732776.

3.1 Coverage distribution for all samples

plotEmpiricalDistribution(bsseq_obj, 
                          bySample = TRUE,
                          testCovariate = NULL,
                          type = "Cov") + guides(linetype="none") + ggplot2::scale_color_manual(values = sample_colors)

plotEmpiricalDistribution(bsseq_obj, 
                          bySample = TRUE,
                          testCovariate = "Type",
                          type = "Cov") + guides(linetype="none") + ggplot2::scale_color_manual(values = cellTypes_colors)

3.2 Methylation levels distribution for all samples

plotEmpiricalDistribution(bsseq_obj, 
                          bySample = TRUE,
                          testCovariate = NULL,
                          adj = 3) + guides(linetype="none") + ggplot2::scale_color_manual(values = sample_colors)

plotEmpiricalDistribution(bsseq_obj, 
                          bySample = TRUE,
                          testCovariate = "Type",
                          adj = 3) + guides(linetype="none") + ggplot2::scale_color_manual(values = cellTypes_colors)

3.3 Violin Plot showing mean

violin_plot(MethPerc_matrix, stat="mean") + ggplot2::scale_fill_manual(values = sample_colors)
## No id variables; using all as measure variables

3.4 Violin Plot showing median

violin_plot(MethPerc_matrix) + ggplot2::scale_fill_manual(values = sample_colors)
## No id variables; using all as measure variables

Mean values for Coverage

apply(Cov_matrix, 2, mean, na.rm = TRUE)

##  hiPSC_rep1  hiPSC_rep2  hiPSC_rep3  hiPSC_rep4  iMeLC_rep1  iMeLC_rep2 
##    42.59339    54.00549    56.53560    63.08595    44.71197    96.15790 
##  iMeLC_rep3  iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1 
##    53.68046    58.53848    42.71647    27.79907    48.45212    91.98228 
## hEGCLC_rep2 hEGCLC_rep3 
##    87.58598    84.05248

Median values for Coverage

apply(Cov_matrix, 2, median, na.rm = TRUE)

##  hiPSC_rep1  hiPSC_rep2  hiPSC_rep3  hiPSC_rep4  iMeLC_rep1  iMeLC_rep2 
##          37          47          49          54          40          83 
##  iMeLC_rep3  iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1 
##          46          51          40          28          46          83 
## hEGCLC_rep2 hEGCLC_rep3 
##          77          76

Mean values for Methylation Proportions

apply(MethPerc_matrix, 2, mean, na.rm = TRUE)

##  hiPSC_rep1  hiPSC_rep2  hiPSC_rep3  hiPSC_rep4  iMeLC_rep1  iMeLC_rep2 
##    46.79386    47.79021    47.15759    47.42451    47.21074    46.07298 
##  iMeLC_rep3  iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1 
##    46.07659    44.61482    43.78058    44.20949    44.09465    44.67909 
## hEGCLC_rep2 hEGCLC_rep3 
##    45.27622    45.21548

Median values for Methylation Proportions

apply(MethPerc_matrix, 2, median, na.rm = TRUE)

##  hiPSC_rep1  hiPSC_rep2  hiPSC_rep3  hiPSC_rep4  iMeLC_rep1  iMeLC_rep2 
##    46.66667    52.94118    52.17391    52.34375    50.00000    50.00000 
##  iMeLC_rep3  iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1 
##    46.87500    43.63636    38.46154    40.74074    42.85714    45.94595 
## hEGCLC_rep2 hEGCLC_rep3 
##    48.67257    47.91667

3.5 PCA plots

3.5.1 All samples

pca_res <- do_PCA(MethPerc_matrix)
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5)

3.5.2 Without EGCLCs

pca_res <- do_PCA(MethPerc_matrix[, sample_anno[!sample_anno$Type%in%c("hEGCLCs"),]$SampleID])
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5)

3.5.3 Only EGCLCs and hiPSCs

pca_res <- do_PCA(MethPerc_matrix[, sample_anno[sample_anno$Type%in%c("hEGCLCs", "hiPSCs"),]$SampleID])
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5)

3.6 Dendogram plot

d <- hclust(dist(t(MethPerc_matrix),
                 method="euclidean"))
#d$labels <- rownames(pData(bsseq_obj))
d$labels <- rownames(sample_anno)
  #paste0(pData(bsseq_obj)$Type, "_", pData(bsseq_obj)$InputDNA)
dend <- as.dendrogram(d)
plot(dend)

4. Filtering BSSeq object by shared CpGs

The number of samples in the BSSeq object is 14.
The number of CpGs in the BSSeq object is 3732776.
The number of CpGs covered by at least 50% of samples in the BSSeq object is 3709122.
The number of CpGs covered by at least 75% of samples in the BSSeq object is 3680616.
The number of CpGs covered by all samples in the BSSeq object is 3481598.

table(sample_anno$Type)

## 
## hEGCLCs  hiPSCs hPGCLCs  iMeLCs 
##       3       4       3       4

The number of CpGs covered by at least 50% of samples of the same type in the BSSeq object is 3682480.
The number of CpGs covered by at least 75% of samples of the same type in the BSSeq object is 3630948.

5. Saving and Session Info

Saving BSSeq objects :

Complete

saveRDS(bsseq_obj, file = paste0 (OutputFolder, "bsseq_obj_complete.rds"))

Filtered by sample

saveRDS(bsseq_obj, file = paste0 (OutputFolder, "bsseq_obj.rds"))

Filtered by sample and by CpGs covered by all samples kept

saveRDS(methylSig::filter_loci_by_group_coverage(bs = bsseq_obj,group_column = "Line", min_samples_per_group = c("CTL08A" = dim(bsseq_obj)[2])), file = paste0 (OutputFolder, "bsseq_obj_sharedbyall.rds"))

With CpGs covered by 75% of all samples kept

saveRDS(methylSig::filter_loci_by_group_coverage(bs = bsseq_obj,group_column = "Line", min_samples_per_group = c("CTL08A" = round(0.75*dim(bsseq_obj)[2], digits=0))), file = paste0 (OutputFolder, "bsseq_obj_sharedby75ofall.rds"))

Session Info

SessionInfo <- sessionInfo()
Date <- date()

Date

## [1] "Thu Jan  9 17:42:23 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

Exploratory Data Analysis

Gaja Matassa

Date: January 09, 2025

1. Loading data

2. Reading bedgraph files with `read.bismark` from bsseq package

2.1 Generation of BSseq object

2.2 Exploration of BSseq object

3. Exploratory plots

3.1 Coverage distribution for all samples

3.2 Methylation levels distribution for all samples

3.3 Violin Plot showing mean

3.4 Violin Plot showing median

3.5 PCA plots

3.5.1 All samples

3.5.2 Without EGCLCs

3.5.3 Only EGCLCs and hiPSCs

3.6 Dendogram plot

4. Filtering BSSeq object by shared CpGs

5. Saving and Session Info

Exploratory Data Analysis

Gaja Matassa

Date: January 09, 2025

1. Loading data

2. Reading bedgraph files with read.bismark from bsseq package

2.1 Generation of BSseq object

2.2 Exploration of BSseq object

3. Exploratory plots

3.1 Coverage distribution for all samples

3.2 Methylation levels distribution for all samples

3.3 Violin Plot showing mean

3.4 Violin Plot showing median

3.5 PCA plots

3.5.1 All samples

3.5.2 Without EGCLCs

3.5.3 Only EGCLCs and hiPSCs

3.6 Dendogram plot

4. Filtering BSSeq object by shared CpGs

5. Saving and Session Info

2. Reading bedgraph files with `read.bismark` from bsseq package