In cimbusch/tissueTreg: TWGBS and RNA-seq data from tissue T regulatory cells from mice

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

TWGBS data

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.

Per sample

We first load necessary libraries and data:

library(ExperimentHub)
library(bsseq)

eh <- ExperimentHub()
BS.obj.ex.fit <- eh[["EH1072"]]

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)

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.

Per tissue / cell type

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"]]

colnames(BS.obj.ex.fit.combined)

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)

RNA-seq data

Per sample

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"]]

EnsemblIDs are given as well as gene symbols:

head(assay(se_rpkms))
head(rowData(se_rpkms))

We can for example obtain the RPKM values for Foxp3 in this way:

assay(se_rpkms)[rowData(se_rpkms)$id_symbol=="Foxp3",]

colData() also contains information about the tissue and cell type the replicate belongs to:

head(colData(se_rpkms))

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"))

sessionInfo()

sessionInfo()

cimbusch/tissueTreg documentation built on Sept. 28, 2023, 5:20 p.m.