...
|
...
|
@@ -26,15 +26,25 @@ knitr::opts_chunk$set(echo = TRUE)
|
26
|
26
|
options(DelayedArray.block.size=4500e6)
|
27
|
27
|
bs <- params$samples
|
28
|
28
|
nSamples <- length(Biobase::sampleNames(bs))
|
|
29
|
+
|
|
30
|
+
|
|
31
|
+
|
29
|
32
|
```
|
30
|
33
|
|
31
|
34
|
## Read information
|
32
|
35
|
```{r read_metrics,echo=FALSE,warning=FALSE,message=FALSE,fig.width=8, fig.height=7,fig.align='center'}
|
33
|
36
|
# Methylation distribution for all the cells
|
34
|
37
|
dat <- scmeth::readmetrics(bs)
|
|
38
|
+SummaryTable <- dat
|
|
39
|
+SummaryTable$sample <- sub('\\..*$', '', SummaryTable$sample)
|
35
|
40
|
o <- order(dat$total,dat$sample)
|
36
|
41
|
sampleOrder <- dat$sample[o]
|
37
|
|
-
|
|
42
|
+Summary_read <- summary(dat[,c('total','mapped','unmapped')])
|
|
43
|
+
|
|
44
|
+
|
|
45
|
+
|
|
46
|
+
|
|
47
|
+
|
38
|
48
|
m <- reshape2::melt(dat[,c("sample","mapped","unmapped")], id.vars="sample",
|
39
|
49
|
variable.name="Mapping_status")
|
40
|
50
|
m$sample <- factor(m$sample,levels=sampleOrder)
|
...
|
...
|
@@ -57,8 +67,10 @@ g
|
57
|
67
|
## Methylation bias plot
|
58
|
68
|
# Methylation distribution for all the cells
|
59
|
69
|
mbiasTableList <- scmeth::mbiasplot(dir=params$mbias)
|
60
|
|
-
|
61
|
70
|
mbiasDf <- do.call(rbind.data.frame, mbiasTableList)
|
|
71
|
+mbiasDf$methylation <- round(mbiasDf$methylation,4)
|
|
72
|
+write.table(mbiasDf, paste0(params$outdir,'/mbias_table.txt'),sep = "\t",quote=FALSE, row.names=FALSE)
|
|
73
|
+
|
62
|
74
|
meanTable <- stats::aggregate(methylation ~ position+read, data=mbiasDf, FUN=mean)
|
63
|
75
|
sdTable <- stats::aggregate(methylation ~ position+read, data=mbiasDf, FUN=sd)
|
64
|
76
|
seTable <- stats::aggregate(methylation ~ position+read, data=mbiasDf, FUN=function(x){sd(x)/sqrt(length(x))})
|
...
|
...
|
@@ -97,6 +109,8 @@ o <- order(covDf$coveredCpgs)
|
97
|
109
|
sampleOrder <- covDf$sample[o]
|
98
|
110
|
covDf$sample <- factor(covDf$sample, levels=sampleOrder)
|
99
|
111
|
|
|
112
|
+SummaryTable$CpG.coverage <- covVec
|
|
113
|
+
|
100
|
114
|
g <- ggplot2::ggplot(covDf, ggplot2::aes_string(x="''", y='coveredCpgs'))
|
101
|
115
|
g <- g+ggplot2::geom_boxplot()
|
102
|
116
|
g <- g+ggplot2::geom_jitter()
|
...
|
...
|
@@ -110,6 +124,7 @@ g
|
110
|
124
|
## Bisulfite conversion rate
|
111
|
125
|
```{r bs_conversion,echo=FALSE,warning=FALSE,message=FALSE,fig.height=5,fig.width=5,fig.align='center'}
|
112
|
126
|
bscDf <- scmeth::bsConversionPlot(bs)
|
|
127
|
+SummaryTable$Bisulfite.converion.Rate <- round(bscDf$bsc,4)
|
113
|
128
|
|
114
|
129
|
g <- ggplot2::ggplot(bscDf, ggplot2::aes_string(x="''", y='bsc'))
|
115
|
130
|
g <- g+ggplot2::geom_boxplot()
|
...
|
...
|
@@ -167,6 +182,10 @@ g
|
167
|
182
|
```{r cpg_discretization,echo=FALSE,warning=FALSE,message=FALSE,fig.height=5, fig.width=5, eval=TRUE,fig.align='center'}
|
168
|
183
|
# CpG Discretization
|
169
|
184
|
discretizedCpG <- cpgDiscretization(bs,subSample=params$nCpGs,offset=params$offset,coverageVec=covVec)
|
|
185
|
+
|
|
186
|
+SummaryTable$CpG.nonBinary.percentage <- round(discretizedCpG$discardPerc,4)
|
|
187
|
+write.table(SummaryTable, paste0(params$outdir,'/QC_Summary.txt'),sep = "\t",quote=FALSE, row.names=FALSE)
|
|
188
|
+
|
170
|
189
|
#percentDiscarded <- discretizedCpG[3]
|
171
|
190
|
percentDiscarded <- discretizedCpG[2]
|
172
|
191
|
sampleNames <- bsseq::sampleNames(bs)
|
...
|
...
|
@@ -212,6 +231,9 @@ g
|
212
|
231
|
##############################
|
213
|
232
|
# Downsample Plots
|
214
|
233
|
dsMatrix <- scmeth::downsample(bs,subSample=params$nCpGs,offset=params$offset)
|
|
234
|
+colnames(dsMatrix) <- SummaryTable$sample
|
|
235
|
+write.table(t(dsMatrix), paste0(params$outdir,'/Downsample_Table.txt'),sep = "\t",quote=FALSE)
|
|
236
|
+
|
215
|
237
|
meltedDsMatrix <- reshape2::melt(dsMatrix)
|
216
|
238
|
g <- ggplot2::ggplot(meltedDsMatrix, ggplot2::aes_string(x='Var1', y='value', group='Var2'))
|
217
|
239
|
g <- g+ggplot2::geom_line(alpha=0.2)
|