/
globalPlots.R
173 lines (160 loc) · 6.9 KB
/
globalPlots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
#' windows
#' @title Extract methylation values from tiled genomic windows
#' @description Obtain windows of individual smoothed methylation values
#' @param bs.filtered.bsseq Smoothed \code{bsseq} object with a testCovariate in \code{pData}
#' @param size The number of bases in the window (default is 2e4, which is 20 Kb)
#' @param goi A \code{BSgenome} object of the genome of interest (i.e. "BSgenome.Hsapiens.UCSC.hg38")
#' @return A matrix of smoothed individual methylation values
#' @importFrom magrittr %>%
#' @importFrom GenomeInfoDb seqlengths keepStandardChromosomes
#' @importFrom GenomicRanges tileGenome
#' @importFrom bsseq getMeth
#' @importClassesFrom bsseq BSseq
#' @importMethodsFrom bsseq pData
#' @export windows
#'
windows <- function(bs.filtered.bsseq = bs.filtered.bsseq,
size = 2e4,
goi = goi){
print(glue::glue("Obtaining {size/1000} Kb window individual smoothed methylation values from the {BSgenome::commonName(goi)} genome"))
goi %>%
GenomeInfoDb::seqlengths() %>%
GenomicRanges::tileGenome(tilewidth = size,
cut.last.tile.in.chrom = TRUE) %>%
GenomeInfoDb::keepStandardChromosomes(pruning.mode = "coarse") %>%
bsseq::getMeth(BSseq = bs.filtered.bsseq,
regions = .,
type = "smooth",
what = "perRegion") %>%
na.omit() %>%
return()
}
#' CGi
#' @title Extract methylation values from CpG islands
#' @description Obtain individual smoothed methylation values for CpG islands
#' @param bs.filtered.bsseq Smoothed \code{bsseq} object with a testCovariate in \code{pData}
#' @param genome A character vector of the genome of interest (i.e. "hg38")
#' @return A matrix of smoothed individual methylation values
#' @importFrom plyranges filter
#' @importFrom magrittr %>%
#' @importFrom bsseq getMeth
#' @importClassesFrom bsseq BSseq
#' @importMethodsFrom bsseq pData
#' @export CGi
#'
CGi <- function(bs.filtered.bsseq = bs.filtered.bsseq,
genome = genome){
print(glue::glue("Obtaining individual smoothed methylation values of CpG islands from {genome}"))
genome %>%
DMRichR::getCpGs() %>%
plyranges::filter(type == "islands") %>%
bsseq::getMeth(BSseq = bs.filtered.bsseq,
regions = .,
type = "smooth",
what = "perRegion") %>%
na.omit() %>%
return()
}
#' CpGs
#' @title Extract single CpG methylation values
#' @description Extracts single CpG individual smoothed methylation values
#' @param bs.filtered.bsseq Smoothed \code{bsseq} object with a testCovariate in \code{pData}
#' @return A matrix of smoothed individual methylation values
#' @importFrom bsseq getMeth
#' @importClassesFrom bsseq BSseq
#' @importMethodsFrom bsseq pData
#' @export CpGs
#'
CpGs <- function(bs.filtered.bsseq = bs.filtered.bsseq){
print(glue::glue("Obtaining smoothed methylation values for all covered CpGs"))
bs.filtered.bsseq %>%
bsseq::getMeth(BSseq = .,
type = "smooth",
what = "perBase") %>%
na.omit() %>%
return()
}
#' PCA
#' @title PCA plot of extracted methylation values
#' @description Performs and plots a PCA from individual smoothed methylation values
#' @param matrix A matrix of smoothed individual methylation values
#' @param testCovariate Factor of interest
#' @param bs.filtered.bsseq Smoothed \code{bsseq} object with a testCovariate in \code{pData}
#' @return A \code{ggplot} object that can be viewed by calling it,
#' saved with \code{ggplot2::ggsave()}, or further modified by adding \code{ggplot2} syntax
#' @importFrom PCAtools pca biplot
#' @importFrom glue glue
#' @importFrom Glimma glMDSPlot
#' @importFrom purrr pluck
#' @importClassesFrom bsseq BSseq
#' @importMethodsFrom bsseq pData
#' @references \url{https://stackoverflow.com/questions/40315227/how-to-solve-prcomp-default-cannot-rescale-a-constant-zero-column-to-unit-var/40317343}
#' @export PCA
#'
PCA <- function(matrix = matrix,
testCovariate = testCovariate,
bs.filtered.bsseq = bs.filtered.bsseq){
print(glue::glue("PCA of {length(matrix)} sites"))
matrix %>%
PCAtools::pca(scale = TRUE,
metadata = pData(bs.filtered.bsseq),
removeVar = 0.1) %>%
PCAtools::biplot(colby = testCovariate,
colkey = DMRichR::gg_color_hue(2) %>%
setNames(bs.filtered.bsseq %>%
pData() %>%
as.data.frame() %>%
purrr::pluck(testCovariate) %>%
unique() %>%
sort() %>%
rev()),
pointSize = 6,
labSize = 4,
legendPosition = 'top',
legendLabSize = 16,
legendIconSize = 8.0)
}
#' densityPlot
#' @title Density plot of extracted methylation values
#' @description Creates a density plot of the mean of individual smoothed methylation values
#' @param matrix A matrix of smoothed individual methylation values
#' @param group Ordered factor vector of sample groupings
#' @return A \code{ggplot} object that can be viewed by calling it,
#' saved with \code{ggplot2::ggsave()}, or further modified by adding \code{ggplot2} syntax.
#' @importFrom magrittr %>% set_colnames
#' @importFrom dplyr as_tibble select transmute contains mutate
#' @importFrom forcats fct_rev
#' @import ggplot2
#' @importFrom tidyr gather
#' @export densityPlot
#'
densityPlot <- function(matrix = matrix,
group = NA){
print(glue::glue("Density plot of {length(matrix)} sites"))
matrix %>%
dplyr::as_tibble() %>%
magrittr::set_colnames(paste(group, seq_along(1:length(group)))) %>%
dplyr::transmute(Group1 = dplyr::select(., dplyr::contains(levels(group)[1])) %>% rowMeans()*100,
Group2 = dplyr::select(., dplyr::contains(levels(group)[2])) %>% rowMeans()*100) %>%
magrittr::set_colnames(c(levels(group)[1], levels(group)[2])) %>%
tidyr::gather(key = "variable",
value = "value") %>%
#dplyr::mutate(variable = factor(.$variable)) %>%
dplyr::mutate(variable = factor(.$variable, levels = levels(group))) %>%
ggplot(aes(value, color = variable)) +
geom_density(size = 1.2) +
labs(x = "Percent Methylation",
y = "Density",
color = "Group") +
theme_classic() +
scale_x_continuous(expand = c(0.05,0.05),
breaks = c(0,25,50,75,100)) +
scale_y_continuous(expand = c(0.00,0.001)) +
theme(axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
strip.text = element_text(size = 20),
legend.text = element_text(size = 20),
legend.position = "bottom",
legend.title = element_text(size = 20)) %>%
return()
}