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

InputFolder <- params$InputFolder
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", 
                   "hiPSC_rep5" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_p5_rep1" = "#ff7f00", "hEGCLC_p5_rep2" = "#ff7f00", "hEGCLC_p5_rep3" = "#ff7f00",
                   "hEGCLC_p10_rep1" = "#ffc400", "hEGCLC_p10_rep2" = "#ffc400", "hEGCLC_p10_rep3" = "#ffc400")

cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs_p5 = "#ff7f00", hEGCLCs_p10 = "#ffc400")

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 <- data.frame("SampleID" = c("hiPSC_rep5", "hEGCLC_p10_rep1", "hEGCLC_p10_rep2", "hEGCLC_p10_rep3"), "OriginalID" = c("hiPSC", "hEGCLC_clone2", "hEGCLC_clone9", "hEGCLC_clone31"), "Type" = c("hiPSCs", "hEGCLCs_p10", "hEGCLCs_p10", "hEGCLCs_p10"), "Line" = rep("CTL08A", 4), "passage" = c("26", "10", "10", "10"))
sample_anno$Type <- factor(sample_anno$Type)
rownames(sample_anno) <- sample_anno$SampleID

Match order

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

Adding previous experiment samples

bdg_files2 <-  list.files(path = "~/DataDir/0.Preprocessing/Output", pattern = '*\\.bedGraph', full.names = TRUE) #output from our preprocessing

bdg_files_filtered2 <- bdg_files2[grep("filtered", bdg_files2)]

sample_anno2 <- read.table("~/DataDir/1.PreliminaryAnalysis/Input/SampleAnno.txt", sep='\t', header=TRUE)

sample_anno2[sample_anno2$Type %in% "hEGCLCs", "Type"] <- "hEGCLCs_p5"
sample_anno2$SampleID <- sub("hEGCLC_", "hEGCLC_p5_", sample_anno2$SampleID)
rownames(sample_anno2) <- sample_anno2$SampleID

Keeping only bedGraph files of selected samples.

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

Match order

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

sample_anno <- plyr::rbind.fill(sample_anno2, sample_anno)
rownames(sample_anno) <- sample_anno$SampleID
bdg_files_filtered <- c(bdg_files_filtered2, bdg_files_filtered)

sample_anno <- rbind(sample_anno[1:4,], sample_anno[15,], sample_anno[5:14,], sample_anno[16:18,])

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 18.
The number of CpGs in the BSSeq object is 3734728.

2.2 Exploration of BSseq object

bsseq::pData(bsseq_obj)

## DataFrame with 18 rows and 12 columns
##                        SampleID     OriginalID        Type       Day     Batch
##                     <character>    <character> <character> <integer> <integer>
## hiPSC_rep1           hiPSC_rep1             2S      hiPSCs        NA         2
## hiPSC_rep2           hiPSC_rep2             3S      hiPSCs        NA         3
## hiPSC_rep3           hiPSC_rep3             5S      hiPSCs        NA         5
## hiPSC_rep4           hiPSC_rep4             4S      hiPSCs        NA         4
## hiPSC_rep5           hiPSC_rep5          hiPSC      hiPSCs        NA        NA
## ...                         ...            ...         ...       ...       ...
## hEGCLC_p5_rep2   hEGCLC_p5_rep2       3_EG9_SL  hEGCLCs_p5        NA         3
## hEGCLC_p5_rep3   hEGCLC_p5_rep3      3_EG31_SL  hEGCLCs_p5        NA         3
## hEGCLC_p10_rep1 hEGCLC_p10_rep1  hEGCLC_clone2 hEGCLCs_p10        NA        NA
## hEGCLC_p10_rep2 hEGCLC_p10_rep2  hEGCLC_clone9 hEGCLCs_p10        NA        NA
## hEGCLC_p10_rep3 hEGCLC_p10_rep3 hEGCLC_clone31 hEGCLCs_p10        NA        NA
##                        Line     Clone      Medium  PercOxigen PercEfficiency
##                 <character> <integer> <character> <character>    <character>
## hiPSC_rep1           CTL08A        NA                                       
## hiPSC_rep2           CTL08A        NA                                       
## hiPSC_rep3           CTL08A        NA                                       
## hiPSC_rep4           CTL08A        NA                                       
## hiPSC_rep5           CTL08A        NA          NA          NA             NA
## ...                     ...       ...         ...         ...            ...
## hEGCLC_p5_rep2       CTL08A         9     StemFit      low_ox               
## hEGCLC_p5_rep3       CTL08A        31     StemFit      low_ox               
## hEGCLC_p10_rep1      CTL08A        NA          NA          NA             NA
## hEGCLC_p10_rep2      CTL08A        NA          NA          NA             NA
## hEGCLC_p10_rep3      CTL08A        NA          NA          NA             NA
##                       Notes     passage
##                 <character> <character>
## hiPSC_rep1                           NA
## hiPSC_rep2                           NA
## hiPSC_rep3                           NA
## hiPSC_rep4                           NA
## hiPSC_rep5               NA          26
## ...                     ...         ...
## hEGCLC_p5_rep2                       NA
## hEGCLC_p5_rep3                       NA
## hEGCLC_p10_rep1          NA          10
## hEGCLC_p10_rep2          NA          10
## hEGCLC_p10_rep3          NA          10

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 hiPSC_rep5 iMeLC_rep1
## [1,]         16         18         19         21         46         18
## [2,]         18         20         22         25         50         24
## [3,]         23         21         23         25         60         27
## [4,]         23         22         25         26         62         22
## [5,]         24         21         26         27         58         27
## [6,]         23         22         27         27         58         25
##      iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,]         28         27         16          13          11          23
## [2,]         32         28         17          17          12          28
## [3,]         39         30         19          19          16          28
## [4,]         40         30         21          18          16          25
## [5,]         40         30         23          20          15          25
## [6,]         42         30         23          19          17          26
##      hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,]             40             31             27              32
## [2,]             45             36             32              37
## [3,]             49             40             44              47
## [4,]             52             43             42              52
## [5,]             51             46             42              53
## [6,]             53             49             47              53
##      hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,]              19              26
## [2,]              22              27
## [3,]              29              31
## [4,]              28              32
## [5,]              28              33
## [6,]              28              32

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 hiPSC_rep5 iMeLC_rep1
## [1,]         19         21         23         25         50         24
## [2,]         18         20         22         25         51         25
## [3,]         23         22         23         26         62         28
## [4,]         23         23         26         27         62         27
## [5,]         24         22         27         29         58         28
## [6,]         23         23         28         30         58         26
##      iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,]         35         31         19          20          16          29
## [2,]         34         29         18          17          12          28
## [3,]         39         32         22          22          16          29
## [4,]         41         32         24          21          16          28
## [5,]         42         32         25          24          16          27
## [6,]         42         32         27          23          17          30
##      hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,]             45             43             37              34
## [2,]             46             37             34              38
## [3,]             54             42             46              49
## [4,]             55             46             47              53
## [5,]             55             48             48              53
## [6,]             56             51             49              54
##      hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,]              24              29
## [2,]              25              27
## [3,]              30              31
## [4,]              29              33
## [5,]              29              33
## [6,]              29              34

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 hiPSC_rep5 iMeLC_rep1
## [1,]   84.21053   85.71429   82.60870   84.00000   92.00000   75.00000
## [2,]  100.00000  100.00000  100.00000  100.00000   98.03922   96.00000
## [3,]  100.00000   95.45455  100.00000   96.15385   96.77419   96.42857
## [4,]  100.00000   95.65217   96.15385   96.29630  100.00000   81.48148
## [5,]  100.00000   95.45455   96.29630   93.10345  100.00000   96.42857
## [6,]  100.00000   95.65217   96.42857   90.00000  100.00000   96.15385
##      iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,]   80.00000   87.09677   84.21053    65.00000       68.75    79.31034
## [2,]   94.11765   96.55172   94.44444   100.00000      100.00   100.00000
## [3,]  100.00000   93.75000   86.36364    86.36364      100.00    96.55172
## [4,]   97.56098   93.75000   87.50000    85.71429      100.00    89.28571
## [5,]   95.23810   93.75000   92.00000    83.33333       93.75    92.59259
## [6,]  100.00000   93.75000   85.18519    82.60870      100.00    86.66667
##      hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,]       88.88889       72.09302       72.97297        94.11765
## [2,]       97.82609       97.29730       94.11765        97.36842
## [3,]       90.74074       95.23810       95.65217        95.91837
## [4,]       94.54545       93.47826       89.36170        98.11321
## [5,]       92.72727       95.83333       87.50000       100.00000
## [6,]       94.64286       96.07843       95.91837        98.14815
##      hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,]        79.16667        89.65517
## [2,]        88.00000       100.00000
## [3,]        96.66667       100.00000
## [4,]        96.55172        96.96970
## [5,]        96.55172       100.00000
## [6,]        96.55172        94.11765

This matrix can also be obtained like this:

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

##      hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 hiPSC_rep5 iMeLC_rep1
## [1,]  0.8421053  0.8571429  0.8260870  0.8400000  0.9200000  0.7500000
## [2,]  1.0000000  1.0000000  1.0000000  1.0000000  0.9803922  0.9600000
## [3,]  1.0000000  0.9545455  1.0000000  0.9615385  0.9677419  0.9642857
## [4,]  1.0000000  0.9565217  0.9615385  0.9629630  1.0000000  0.8148148
## [5,]  1.0000000  0.9545455  0.9629630  0.9310345  1.0000000  0.9642857
## [6,]  1.0000000  0.9565217  0.9642857  0.9000000  1.0000000  0.9615385
##      iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,]  0.8000000  0.8709677  0.8421053   0.6500000      0.6875   0.7931034
## [2,]  0.9411765  0.9655172  0.9444444   1.0000000      1.0000   1.0000000
## [3,]  1.0000000  0.9375000  0.8636364   0.8636364      1.0000   0.9655172
## [4,]  0.9756098  0.9375000  0.8750000   0.8571429      1.0000   0.8928571
## [5,]  0.9523810  0.9375000  0.9200000   0.8333333      0.9375   0.9259259
## [6,]  1.0000000  0.9375000  0.8518519   0.8260870      1.0000   0.8666667
##      hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,]      0.8888889      0.7209302      0.7297297       0.9411765
## [2,]      0.9782609      0.9729730      0.9411765       0.9736842
## [3,]      0.9074074      0.9523810      0.9565217       0.9591837
## [4,]      0.9454545      0.9347826      0.8936170       0.9811321
## [5,]      0.9272727      0.9583333      0.8750000       1.0000000
## [6,]      0.9464286      0.9607843      0.9591837       0.9814815
##      hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,]       0.7916667       0.8965517
## [2,]       0.8800000       1.0000000
## [3,]       0.9666667       1.0000000
## [4,]       0.9655172       0.9696970
## [5,]       0.9655172       1.0000000
## [6,]       0.9655172       0.9411765

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 3734728.

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 = FALSE,
                          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 = FALSE,
                          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      hiPSC_rep5 
##        42.57113        53.97726        56.50605        63.05298       227.83526 
##      iMeLC_rep1      iMeLC_rep2      iMeLC_rep3      iMeLC_rep4     hPGCLC_rep1 
##        44.68860        96.10764        53.65240        58.50788        42.69414 
##     hPGCLC_rep2     hPGCLC_rep3  hEGCLC_p5_rep1  hEGCLC_p5_rep2  hEGCLC_p5_rep3 
##        27.78454        48.42680        91.93421        87.54020        84.00855 
## hEGCLC_p10_rep1 hEGCLC_p10_rep2 hEGCLC_p10_rep3 
##       182.41638       163.61658       162.38895

Median values for Coverage

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

##      hiPSC_rep1      hiPSC_rep2      hiPSC_rep3      hiPSC_rep4      hiPSC_rep5 
##              37              47              49              54             203 
##      iMeLC_rep1      iMeLC_rep2      iMeLC_rep3      iMeLC_rep4     hPGCLC_rep1 
##              40              83              46              51              40 
##     hPGCLC_rep2     hPGCLC_rep3  hEGCLC_p5_rep1  hEGCLC_p5_rep2  hEGCLC_p5_rep3 
##              28              46              83              77              76 
## hEGCLC_p10_rep1 hEGCLC_p10_rep2 hEGCLC_p10_rep3 
##             162             148             149

Mean values for Methylation Proportions

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

##      hiPSC_rep1      hiPSC_rep2      hiPSC_rep3      hiPSC_rep4      hiPSC_rep5 
##        46.79386        47.79021        47.15759        47.42451        48.12409 
##      iMeLC_rep1      iMeLC_rep2      iMeLC_rep3      iMeLC_rep4     hPGCLC_rep1 
##        47.21074        46.07298        46.07659        44.61482        43.78058 
##     hPGCLC_rep2     hPGCLC_rep3  hEGCLC_p5_rep1  hEGCLC_p5_rep2  hEGCLC_p5_rep3 
##        44.20949        44.09465        44.67909        45.27622        45.21548 
## hEGCLC_p10_rep1 hEGCLC_p10_rep2 hEGCLC_p10_rep3 
##        47.52772        48.62960        47.20482

Median values for Methylation Proportions

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

##      hiPSC_rep1      hiPSC_rep2      hiPSC_rep3      hiPSC_rep4      hiPSC_rep5 
##        46.66667        52.94118        52.17391        52.34375        53.96825 
##      iMeLC_rep1      iMeLC_rep2      iMeLC_rep3      iMeLC_rep4     hPGCLC_rep1 
##        50.00000        50.00000        46.87500        43.63636        38.46154 
##     hPGCLC_rep2     hPGCLC_rep3  hEGCLC_p5_rep1  hEGCLC_p5_rep2  hEGCLC_p5_rep3 
##        40.74074        42.85714        45.94595        48.67257        47.91667 
## hEGCLC_p10_rep1 hEGCLC_p10_rep2 hEGCLC_p10_rep3 
##        50.00000        51.53218        50.00000

3.5 PCA plots

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, show_labels = TRUE)

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 18.
The number of CpGs in the BSSeq object is 3734728.
The number of CpGs covered by at least 50% of samples in the BSSeq object is 3706258.
The number of CpGs covered by at least 75% of samples in the BSSeq object is 3644178.
The number of CpGs covered by all samples in the BSSeq object is 3417412.

table(sample_anno$Type)

## 
## hEGCLCs_p10  hEGCLCs_p5      hiPSCs     hPGCLCs      iMeLCs 
##           3           3           5           3           4

5. Saving and Session Info

Saving BSSeq objects :

Complete

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 Jul 24 14:23:00 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: July 24, 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.6 Dendogram plot

4. Filtering BSSeq object by shared CpGs

5. Saving and Session Info

Exploratory Data Analysis

Gaja Matassa

Date: July 24, 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.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