The DNA methylation data have been stored as rds files and are bsseq objects. They contain the raw counts for each CpG position in the mm10 reference genome as well the preprocessed smoothed values for each sample. In addition there is a second file containing DNA methylation profiles per tissue, e.g. the replicates from the same tissue have been collapsed. Also here the smoothed values have been precomputed and are directly available.
We first load necessary libraries and data:
library(ExperimentHub)
library(bsseq)
eh <- ExperimentHub()
BS.obj.ex.fit <- eh[["EH1072"]]
#> see ?tissueTreg and browseVignettes('tissueTreg') for documentation
#> loading from cache
For example to visually inspect the DNA methylation state of the FoxP3 gene in the CNS2 region we first specify a single region:
regions <- GRanges(
seqnames = c("X"),
ranges = IRanges(start = c(7579676),
end = c(7595243)
)
)
We can then directly use the plotRegion function from the bsseq package for visualization based on smoothed values:
plotRegion(BS.obj.ex.fit, region=regions[1,], extend = 2000)
We assume that the three sample with higher DNA methylation in CNS2 are the T conventional samples. To check that visually we can color the samples by tissue and cell type.
We check the order in the object first:
colnames(BS.obj.ex.fit)
#> [1] "Fat-Treg-R1" "Fat-Treg-R2" "Fat-Treg-R3" "Liver-Treg-R1"
#> [5] "Liver-Treg-R2" "Liver-Treg-R3" "Skin-Treg-R1" "Skin-Treg-R2"
#> [9] "Skin-Treg-R3" "Lymph-N-Tcon-R1" "Lymph-N-Tcon-R2" "Lymph-N-Tcon-R3"
#> [13] "Lymph-N-Treg-R1" "Lymph-N-Treg-R2" "Lymph-N-Treg-R3"
And assign the colors:
pData <- pData(BS.obj.ex.fit)
pData$col <- rep(c("red", "blue", "green", "yellow", "orange"), rep(3,5))
pData(BS.obj.ex.fit) <- pData
plotRegion(BS.obj.ex.fit, region=regions[1,], extend = 2000)
This plot is confirming our assumption but we don’t have to plot all samples since they seem to be already consistent. We would now like to do the same region using smoothed values for each group.
The smoothed data has already been precomputed and stored in a another rds file which we first need to load:
BS.obj.ex.fit.combined <- eh[["EH1073"]]
#> see ?tissueTreg and browseVignettes('tissueTreg') for documentation
#> loading from cache
colnames(BS.obj.ex.fit.combined)
#> [1] "Fat-Treg" "Liver-Treg" "Lymph-N-Tcon" "Lymph-N-Treg" "Skin-Treg"
We now only have the tissue and cell type instead of the replicates. We assign the same colors:
pData <- pData(BS.obj.ex.fit.combined)
pData$col <- c("red", "blue", "yellow", "orange", "green")
pData(BS.obj.ex.fit.combined) <- pData
plotRegion(BS.obj.ex.fit.combined, region=regions[1,], extend = 2000)
The RNA-seq experiments are available and as a SummarizedExperiment object. Two objects are available for usage: RPKM values and htseq counts. We will use the RPKM values and would like to look up the expression of FoxP3:
We load the SummarizedExperiment object:
se_rpkms <- eh[["EH1074"]]
#> see ?tissueTreg and browseVignettes('tissueTreg') for documentation
#> loading from cache
EnsemblIDs are given as well as gene symbols:
head(assay(se_rpkms))
#> Fat-Treg-R1 Fat-Treg-R2 Fat-Treg-R3 Liver-Treg-R1
#> ENSMUSG00000030105 29 11 19 6
#> ENSMUSG00000042428 0 0 0 0
#> ENSMUSG00000096054 2 1 3 1
#> ENSMUSG00000046532 5 5 3 0
#> ENSMUSG00000020495 2 0 2 1
#> ENSMUSG00000058979 2 4 4 7
#> Liver-Treg-R2 Liver-Treg-R3 Lymph-N-Tcon-R1 Lymph-N-Tcon-R2
#> ENSMUSG00000030105 11 16 12 19
#> ENSMUSG00000042428 0 0 0 0
#> ENSMUSG00000096054 1 2 1 1
#> ENSMUSG00000046532 0 0 2 1
#> ENSMUSG00000020495 2 2 2 3
#> ENSMUSG00000058979 14 7 8 8
#> Lymph-N-Tcon-R3 Lymph-N-Treg-R1 Lymph-N-Treg-R2
#> ENSMUSG00000030105 19 18 22
#> ENSMUSG00000042428 0 0 0
#> ENSMUSG00000096054 1 1 2
#> ENSMUSG00000046532 1 0 0
#> ENSMUSG00000020495 2 3 2
#> ENSMUSG00000058979 14 8 8
#> Lymph-N-Treg-R3 Skin-Treg-R1 Skin-Treg-R2 Skin-Treg-R3
#> ENSMUSG00000030105 21 16 19 24
#> ENSMUSG00000042428 0 0 0 0
#> ENSMUSG00000096054 1 1 2 3
#> ENSMUSG00000046532 0 0 1 1
#> ENSMUSG00000020495 2 1 1 1
#> ENSMUSG00000058979 10 6 3 8
head(rowData(se_rpkms))
#> DataFrame with 6 rows and 1 column
#> id_symbol
#> <character>
#> ENSMUSG00000030105 Arl8b
#> ENSMUSG00000042428 Mgat3
#> ENSMUSG00000096054 Syne1
#> ENSMUSG00000046532 Ar
#> ENSMUSG00000020495 Smg8
#> ENSMUSG00000058979 Cecr5
We can for example obtain the RPKM values for Foxp3 in this way:
assay(se_rpkms)[rowData(se_rpkms)$id_symbol=="Foxp3",]
#> Fat-Treg-R1 Fat-Treg-R2 Fat-Treg-R3 Liver-Treg-R1 Liver-Treg-R2
#> 14 16 20 13 14
#> Liver-Treg-R3 Lymph-N-Tcon-R1 Lymph-N-Tcon-R2 Lymph-N-Tcon-R3 Lymph-N-Treg-R1
#> 18 1 1 1 15
#> Lymph-N-Treg-R2 Lymph-N-Treg-R3 Skin-Treg-R1 Skin-Treg-R2 Skin-Treg-R3
#> 13 15 25 14 21
colData() also contains information about the tissue and cell type the replicate belongs to:
head(colData(se_rpkms))
#> DataFrame with 6 rows and 1 column
#> tissue_cell
#> <character>
#> Fat-Treg-R1 Fat-Treg
#> Fat-Treg-R2 Fat-Treg
#> Fat-Treg-R3 Fat-Treg
#> Liver-Treg-R1 Liver-Treg
#> Liver-Treg-R2 Liver-Treg
#> Liver-Treg-R3 Liver-Treg
We can put this information together and visualize a simple barplot:
library(ggplot2)
library(reshape2)
foxp3_rpkm <- assay(se_rpkms)[rowData(se_rpkms)$id_symbol=="Foxp3",]
foxp3_rpkm_molten <- melt(foxp3_rpkm)
ggplot(data=foxp3_rpkm_molten, aes(x=rownames(foxp3_rpkm_molten), y=value, fill=colData(se_rpkms)$tissue_cell)) +
geom_bar(stat="identity") +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
xlab("Sample") +
ylab("RPKM") +
ggtitle("FoxP3 expression") +
guides(fill=guide_legend(title="tissue / cell type"))
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] reshape2_1.4.4 ggplot2_3.4.4
#> [3] DelayedArray_0.28.0 tissueTreg_1.22.0
#> [5] bsseq_1.38.0 SummarizedExperiment_1.32.0
#> [7] Biobase_2.62.0 MatrixGenerics_1.14.0
#> [9] matrixStats_1.0.0 GenomicRanges_1.54.0
#> [11] GenomeInfoDb_1.38.0 IRanges_2.36.0
#> [13] S4Vectors_0.40.0 ExperimentHub_2.10.0
#> [15] AnnotationHub_3.10.0 BiocFileCache_2.10.0
#> [17] dbplyr_2.3.4 BiocGenerics_0.48.0
#> [19] BiocStyle_2.30.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.1.3 bitops_1.0-7
#> [3] permute_0.9-7 rlang_1.1.1
#> [5] magrittr_2.0.3 compiler_4.3.1
#> [7] RSQLite_2.3.1 DelayedMatrixStats_1.24.0
#> [9] png_0.1-8 vctrs_0.6.4
#> [11] stringr_1.5.0 pkgconfig_2.0.3
#> [13] crayon_1.5.2 fastmap_1.1.1
#> [15] magick_2.8.1 XVector_0.42.0
#> [17] ellipsis_0.3.2 labeling_0.4.3
#> [19] utf8_1.2.4 Rsamtools_2.18.0
#> [21] promises_1.2.1 rmarkdown_2.25
#> [23] purrr_1.0.2 bit_4.0.5
#> [25] xfun_0.40 zlibbioc_1.48.0
#> [27] cachem_1.0.8 jsonlite_1.8.7
#> [29] blob_1.2.4 later_1.3.1
#> [31] rhdf5filters_1.14.0 Rhdf5lib_1.24.0
#> [33] BiocParallel_1.36.0 interactiveDisplayBase_1.40.0
#> [35] parallel_4.3.1 R6_2.5.1
#> [37] stringi_1.7.12 bslib_0.5.1
#> [39] rtracklayer_1.62.0 limma_3.58.0
#> [41] jquerylib_0.1.4 Rcpp_1.0.11
#> [43] bookdown_0.36 knitr_1.44
#> [45] R.utils_2.12.2 httpuv_1.6.12
#> [47] Matrix_1.6-1.1 tidyselect_1.2.0
#> [49] abind_1.4-5 yaml_2.3.7
#> [51] codetools_0.2-19 curl_5.1.0
#> [53] plyr_1.8.9 lattice_0.22-5
#> [55] tibble_3.2.1 withr_2.5.1
#> [57] shiny_1.7.5.1 KEGGREST_1.42.0
#> [59] evaluate_0.22 Biostrings_2.70.1
#> [61] pillar_1.9.0 BiocManager_1.30.22
#> [63] filelock_1.0.2 generics_0.1.3
#> [65] RCurl_1.98-1.12 BiocVersion_3.18.0
#> [67] sparseMatrixStats_1.14.0 munsell_0.5.0
#> [69] scales_1.2.1 gtools_3.9.4
#> [71] xtable_1.8-4 glue_1.6.2
#> [73] tools_4.3.1 BiocIO_1.12.0
#> [75] data.table_1.14.8 GenomicAlignments_1.38.0
#> [77] BSgenome_1.70.0 locfit_1.5-9.8
#> [79] XML_3.99-0.14 rhdf5_2.46.0
#> [81] grid_4.3.1 AnnotationDbi_1.64.0
#> [83] colorspace_2.1-0 GenomeInfoDbData_1.2.11
#> [85] HDF5Array_1.30.0 restfulr_0.0.15
#> [87] cli_3.6.1 rappdirs_0.3.3
#> [89] fansi_1.0.5 S4Arrays_1.2.0
#> [91] dplyr_1.1.3 gtable_0.3.4
#> [93] R.methodsS3_1.8.2 sass_0.4.7
#> [95] digest_0.6.33 SparseArray_1.2.0
#> [97] farver_2.1.1 rjson_0.2.21
#> [99] R.oo_1.25.0 memoise_2.0.1
#> [101] htmltools_0.5.6.1 lifecycle_1.0.3
#> [103] httr_1.4.7 statmod_1.5.0
#> [105] mime_0.12 bit64_4.0.5