Some Covered CpG sites is extremely high for RRBS with BiSeq package

[reprogrammer]

I am using BiSeq package following bismark tools for RRBS analysis. I termed Covfiles as a BSraw file to include all the used 20 samples, and below is what it returns when I call covStatistics(Covfiles)

> Covfiles
class: BSraw
dim: 5542384 20
metadata(0):
assays(2): totalReads methReads
rownames(5542384): 1 2 ... 5542383 5542384
rowData names(0):
colnames(20): 5AD0 5AD1 ... PBD5 PBD7
colData names(1): group

> covStatistics(Covfiles)
$Covered_CpG_sites
5AD0 5AD1 5AD3 5AD5 5AD7 FAD0 FAD1 FAD3
3202085 3176893 3068230 3138142 3046918 1352823 3302600 2980850
FAD5 FAD7 NSD0 NSD1 NSD3 NSD5 NSD7 PBD0
2997963 3098611 3400952 3413895 2913465 3164717 3108568 3114307
PBD1 PBD3 PBD5 PBD7
2392880 3189085 2932880 3188629

$Median_coverage
5AD0 5AD1 5AD3 5AD5 5AD7 FAD0 FAD1 FAD3 FAD5 FAD7 NSD0 NSD1 NSD3 NSD5
5 5 5 5 4 3 4 4 4 5 4 3 2 2
NSD7 PBD0 PBD1 PBD3 PBD5 PBD7
5 6 3 7 5 6


But I found some site is covered extremely higher than others when I use

covBoxplots(Covfiles, col = "cornflowerblue", las = 2)

I don't quite sure whether this is right or not. Is anyone could tell me how is my data looks like, and why there are so many strange sites?


> sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.5 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel stats4 stats graphics grDevices utils
[7] datasets methods base

other attached packages:
[1] BiSeq_1.12.0 Formula_1.2-1
[3] SummarizedExperiment_1.2.2 Biobase_2.32.0
[5] GenomicRanges_1.24.2 GenomeInfoDb_1.8.1
[7] IRanges_2.6.0 S4Vectors_0.10.1
[9] BiocGenerics_0.18.0

loaded via a namespace (and not attached):
[1] XVector_0.12.0 bitops_1.0-6
[3] tools_3.3.0 zlibbioc_1.18.0
[5] annotate_1.50.0 RSQLite_1.0.0
[7] lattice_0.20-33 Matrix_1.2-6
[9] DBI_0.4-1 rtracklayer_1.32.0
[11] Biostrings_2.40.2 lmtest_0.9-34
[13] grid_3.3.0 nnet_7.3-12
[15] globaltest_5.26.0 flexmix_2.3-13
[17] AnnotationDbi_1.34.3 XML_3.98-1.4
[19] survival_2.39-2 BiocParallel_1.6.2
[21] lokern_1.1-6 Rsamtools_1.24.0
[23] modeltools_0.2-21 sfsmisc_1.1-0
[25] GenomicAlignments_1.8.3 splines_3.3.0
[27] xtable_1.8-2 betareg_3.0-5
[29] sandwich_2.3-4 RCurl_1.95-4.8
[31] zoo_1.7-13

[dpryan]

That happens quite a bit. My guess is that those are PCR duplicates. I've normally filtered out the top ~0.1% of covered sites from RRBS datasets.

[reprogrammer]

That happens quite a bit. My guess is that those are PCR duplicates. I've normally filtered out the top ~0.1% of covered sites from RRBS datasets.

Thank you Ryan, May I ask how do you filter the top ~0.1% of covered sites? Is it also done by BiSeq package?

[dpryan]

I usually do it in R, so presumably one can do that with the BSraw object.