diff --git a/.Rbuildignore b/.Rbuildignore
index 01b8059..b5171a0 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -12,4 +12,5 @@
 ^\.RData$
 ^\.Ruserdata$
 ^\.Rproj\.user$
-^\.github
\ No newline at end of file
+^\.github
+^CITATION\.cff$
diff --git a/CITATION.cff b/CITATION.cff
new file mode 100644
index 0000000..2e50e64
--- /dev/null
+++ b/CITATION.cff
@@ -0,0 +1,33 @@
+cff-version: 1.2.0
+message: "If you use methylTFR, please cite it as below."
+type: software
+title: "methylTFR: Quantification of DNA methylation signatures in TFBS"
+abstract: >-
+  methylTFR is an R/Bioconductor package for analyzing DNA methylation
+  signatures within transcription factor binding sites (TFBS) in individual
+  cells or samples. It computes a GC-bias-corrected deviation score of
+  methylation levels at TF binding sites.
+authors:
+  - family-names: Gunduz
+    given-names: Irem B.
+    email: irembgunduz@gmail.com
+    orcid: "https://orcid.org/0000-0003-2641-0916"
+  - family-names: Murugan
+    given-names: Sarath Kumar
+    orcid: "https://orcid.org/0000-0001-7378-7869"
+  - family-names: Mueller
+    given-names: Fabian
+    orcid: "https://orcid.org/0000-0001-5809-2321"
+version: 0.99.8
+date-released: "2026-08-19"
+doi: 10.18129/B9.bioc.methylTFR
+license: MIT
+repository-code: "https://github.com/EpigenomeInformatics/methylTFR"
+url: "https://epigenomeinformatics.github.io/methylTFR/"
+keywords:
+  - R
+  - Bioconductor
+  - DNA methylation
+  - transcription factor binding sites
+  - single-cell
+  - epigenetics
diff --git a/DESCRIPTION b/DESCRIPTION
index 872da72..f95b550 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,18 @@
 Package: methylTFR
 Title: Quantification of DNA Methylation Signatures in TFBS
-Version: 0.99.7
+Version: 0.99.8
 Date: 2026-08-19
 Authors@R: c(
     person("Irem B.", "Gündüz", , "irembgunduz@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-2641-0916")),
     person("Sarath Kumar", "Murugan", , "sarath.murugan@outlook.com", role = "aut",
            comment = c(ORCID = "0000-0001-7378-7869")),
-    person("Fabian", "Muller", , "fabian.mueller@uni-saarland.de", role = "aut",
-           comment = c(ORCID = "0000-0001-5809-2321"))
+    person("Fabian", "Mueller", , "fabian.mueller@uni-saarland.de", role = "aut",
+           comment = c(ORCID = "0000-0001-5809-2321")),
+    person("ERA-NET Transcan-Neu III - EPILUNAR", role = "fnd",
+           comment = "Grant 01KT2409"),
+    person("Saarland University", role = "fnd",
+           comment = "NanoBioMed Young Investigator Grant")
   )
 Description: methylTFR is an extensive R package designed specifically for
     analyzing DNA methylation signatures within transcription factor
diff --git a/NEWS.md b/NEWS.md
index f9b6f94..ac24e64 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,27 @@
+# methylTFR 0.99.8
+
+BUG FIXES
+
+* `read_methylome(type = "encode")` returned methylation scores that were not
+  fractions. `parse_encode()` divided column 11 by column 10, but in a bedMethyl
+  file column 10 is read coverage and column 11 is the percentage of reads
+  methylated, so a fully methylated CpG covered by five reads scored 20. The
+  score is now the percentage over 100.
+* `parse_encode()` skipped the first line of every file unconditionally. Files
+  downloaded from ENCODE carry no header, so the first record was silently
+  discarded. The header is now detected.
+* `read_methylome()` stops if the parsed scores fall outside 0 to 1 rather than
+  passing them to the deviation calculation.
+* The `run_methylTFR_RnBeads()` example built an `RnBiseqSet` without setting
+  `identifiers.column`, so `samples()` returned row numbers and the example
+  stopped on the sample-identifier check. It now sets the option before
+  building the set and restores it afterwards.
+* The bundled `encode.tsv.gz` example has been rebuilt. Columns 10 and 11 are
+  now named `coverage` and `percentMeth` rather than `blockCount` and
+  `blockSizes`, and the percentages span 0 to 100, so the tests exercise the
+  case the old example passed by luck.
+
+
 # methylTFR 0.99.7
 
 NEW FEATURES
diff --git a/R/data_reader.R b/R/data_reader.R
index a57365f..3a37b76 100644
--- a/R/data_reader.R
+++ b/R/data_reader.R
@@ -61,6 +61,21 @@ read_methylome <- function(filename, type, cov_threshold = 1) {
 
     # Filter low coverage sites
     gr_obj <- gr_obj[gr_obj$coverage >= cov_threshold]
+
+    # A methylation score is a fraction. Anything else means the file did not
+    # have the layout the requested type expects, and returning it silently
+    # would put nonsense into every deviation computed downstream.
+    if (length(gr_obj) > 0) {
+        rng <- range(gr_obj$score, na.rm = TRUE)
+        if (is.finite(rng[1]) && (rng[1] < 0 || rng[2] > 1)) {
+            stop(
+                "Methylation scores parsed from ", filename, " span ",
+                signif(rng[1], 3), " to ", signif(rng[2], 3),
+                ", which is not a fraction. Check that '", type,
+                "' is the right format for this file."
+            )
+        }
+    }
     return(gr_obj)
 }
 
@@ -131,17 +146,30 @@ parse_bismarkcov <- function(filename) {
 
 #' @keywords internal
 parse_encode <- function(filename) {
+    # bedMethyl: column 10 is read coverage and column 11 is the percentage of
+    # reads methylated. Columns are taken by position, and the header is
+    # detected rather than assumed, because files downloaded from ENCODE carry
+    # no header line while the bundled example does.
     msites <- data.table::fread(
         filename,
-        header = FALSE, skip = 1, showProgress = FALSE
+        header = "auto", showProgress = FALSE
     )
     if (ncol(msites) < 11) {
         logger::log_warn(sprintf("%s is an invalid encode file!", filename))
         stop("encode file must contain 11 columns!")
     }
-    mscore <- round(msites$V11 / msites$V10, 6)
-    cov <- msites$V10
-    msites <- msites[, .(V1, V2, V3, V6)]
+    cov <- as.numeric(msites[[10]])
+    pct <- as.numeric(msites[[11]])
+    if (any(pct > 100, na.rm = TRUE) || any(pct < 0, na.rm = TRUE)) {
+        stop(
+            "Column 11 of ", filename, " is outside 0-100, so it is not a ",
+            "methylation percentage. Check that this is a bedMethyl file."
+        )
+    }
+    mscore <- round(pct / 100, 6)
+    msites <- data.table::data.table(
+        V1 = msites[[1]], V2 = msites[[2]], V3 = msites[[3]], V6 = msites[[6]]
+    )
     list(msites = msites, mscore = mscore, cov = cov)
 }
 
diff --git a/R/rnbeads_interface.R b/R/rnbeads_interface.R
index a8c4cfd..eefc707 100644
--- a/R/rnbeads_interface.R
+++ b/R/rnbeads_interface.R
@@ -79,6 +79,12 @@ resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) {
 #'
 #' if (requireNamespace("RnBeads", quietly = TRUE) &&
 #'     requireNamespace("RnBeads.hg38", quietly = TRUE)) {
+#'     # identifiers.column is a session option, not a property of the object,
+#'     # so it has to be set before the set is built or samples() falls back
+#'     # to row numbers and methylTFR cannot name the columns it returns
+#'     old_ids <- RnBeads::rnb.getOption("identifiers.column")
+#'     RnBeads::rnb.options(identifiers.column = "sampleName")
+#'
 #'     sites <- data.frame(
 #'         chr = as.character(GenomicRanges::seqnames(msites)),
 #'         start = GenomicRanges::start(msites),
@@ -88,8 +94,10 @@ resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) {
 #'     rnb_set <- RnBeads::RnBiseqSet(
 #'         pheno = data.frame(sampleName = "sample_1"),
 #'         sites = sites,
-#'         meth = matrix(msites$score, ncol = 1),
-#'         covg = matrix(msites$coverage, ncol = 1),
+#'         meth = matrix(msites$score, ncol = 1,
+#'                       dimnames = list(NULL, "sample_1")),
+#'         covg = matrix(msites$coverage, ncol = 1,
+#'                       dimnames = list(NULL, "sample_1")),
 #'         assembly = "hg38",
 #'         summarize.regions = FALSE
 #'     )
@@ -100,6 +108,7 @@ resolve_rnb_sample_ann <- function(rnb_set, sample_ann, sample_ids) {
 #'         gcfreqs = gcfreqs,
 #'         gc_dist = gcdist
 #'     )
+#'     RnBeads::rnb.options(identifiers.column = old_ids)
 #'     deviations(devs)
 #' }
 #' @seealso \code{\link{run_methyltfr}} for running methylTFR from
@@ -149,29 +158,57 @@ run_methylTFR_RnBeads <- function(
 #' @return A character vector of sample identifiers.
 #' @keywords internal
 rnb_sample_ids <- function(rnb_set) {
-    ids <- tryCatch(
-        colnames(RnBeads::meth(rnb_set, type = "sites", i = 1L)),
-        error = function(e) NULL
-    )
+    # Same as RnBeads' samples(): the column named by the
+    # "identifiers.column" option. samples() is not exported by all
+    # RnBeads versions, so it is not called directly.
+    ids <- tryCatch({
+        id_col <- RnBeads::rnb.getOption("identifiers.column")
+        ph <- RnBeads::pheno(rnb_set)
+        if (!is.null(id_col) && id_col %in% colnames(ph)) ph[[id_col]]
+    }, error = function(e) NULL)
     if (length(ids) == 0) {
         ids <- tryCatch(
-            rownames(RnBeads::pheno(rnb_set)),
+            colnames(RnBeads::meth(rnb_set, type = "sites", i = 1L)),
             error = function(e) NULL
         )
     }
     if (length(ids) == 0) {
-        nsamples <- tryCatch(
-            nrow(RnBeads::pheno(rnb_set)),
-            error = function(e) 0L
+        ids <- tryCatch(
+            rownames(RnBeads::pheno(rnb_set)),
+            error = function(e) NULL
         )
-        if (length(nsamples) == 1 && nsamples > 0) {
-            ids <- paste0("sample_", seq_len(nsamples))
-        }
     }
     if (length(ids) == 0) {
         stop("Could not determine sample identifiers from the RnBSet object.")
     }
-    return(as.character(ids))
+    ids <- as.character(ids)
+    if (any(grepl("\\.idat$|_Grn|_Red", ids, ignore.case = TRUE))) {
+        stop(
+            "Sample identifiers look like IDAT channel files, not samples: '",
+            ids[1], "'. The RnBSet was imported with the raw channel files ",
+            "as samples, so each array appears twice and the methylation ",
+            "values are not comparable. Re-import with a sample sheet and ",
+            "an identifiers.column before running methylTFR."
+        )
+    }
+    if (identical(ids, as.character(seq_along(ids)))) {
+        stop(
+            "Sample identifiers are the row numbers 1..", length(ids),
+            ", so the RnBSet carries no usable sample names in this ",
+            "session. identifiers.column is a session option rather than ",
+            "a property of the saved object, so set it before loading, ",
+            "for example rnb.options(identifiers.column = 'Sample.ID'), ",
+            "or assign colnames on the methylation matrix yourself."
+        )
+    }
+    if (anyDuplicated(ids)) {
+        stop(
+            "The RnBSet has duplicated sample identifiers, for example '",
+            ids[anyDuplicated(ids)], "'. Resolve or average the duplicates ",
+            "before running methylTFR."
+        )
+    }
+    return(ids)
 }
 
 
diff --git a/README.Rmd b/README.Rmd
index 1c89905..bfe86c2 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -57,7 +57,7 @@ library(methylTFRAnnotationHg38) # annotation package for hg38
 library(methylTFR)
 
 gcfreqs <- getGCfreq(motifSet = "jaspar2020")
-gc_dist <- getGenomeGC("hg38")
+gc_dist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 sample_dir <- file.path("samples_dir")
@@ -78,9 +78,19 @@ deviations <- run_methyltfr(sample_ann, # sample annotation file
 
 ## Citation
 
-Below is the citation output from using `citation('methylTFR')` in R. Please
-run this yourself to check for any updates on how to cite __methylTFR__.
+If you use `methylTFR` in your work, please cite:
 
-```{r 'citation', eval = requireNamespace('methylTFR')}
-print(citation("methylTFR"), bibtex = TRUE)
+> Gunduz IB, Murugan SK, Mueller F (2026). methylTFR: Quantification of DNA methylation signatures in TFBS. R package version 0.99.8. https://github.com/EpigenomeInformatics/methylTFR
+
+```bibtex
+@Manual{methylTFR,
+  title = {methylTFR: Quantification of DNA methylation signatures in TFBS},
+  author = {Irem B. Gunduz and Sarath Kumar Murugan and Fabian Mueller},
+  year = {2026},
+  note = {R package version 0.99.8},
+  url = {https://github.com/EpigenomeInformatics/methylTFR},
+  doi = {10.18129/B9.bioc.methylTFR},
+}
 ```
+
+Run `citation("methylTFR")` in R to get the entry for the version you have installed.
diff --git a/README.md b/README.md
index 7560092..9d94f79 100644
--- a/README.md
+++ b/README.md
@@ -52,7 +52,7 @@ library(methylTFRAnnotationHg38) # annotation package for hg38
 library(methylTFR)
 
 gcfreqs <- getGCfreq(motifSet = "jaspar2020")
-gc_dist <- getGenomeGC("hg38")
+gc_dist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 sample_dir <- file.path("samples_dir")
@@ -73,27 +73,19 @@ deviations <- run_methyltfr(sample_ann, # sample annotation file
 
 ## Citation
 
-Below is the citation output from using `citation('methylTFR')` in R. Please
-run this yourself to check for any updates on how to cite __methylTFR__.
+If you use `methylTFR` in your work, please cite:
 
+> Gunduz IB, Murugan SK, Mueller F (2026). methylTFR: Quantification of DNA methylation signatures in TFBS. R package version 0.99.8. https://github.com/EpigenomeInformatics/methylTFR
 
-```r
-print(citation("methylTFR"), bibtex = TRUE)
-#> 
-#> To cite package 'methylTFR' in publications use:
-#> 
-#>   Irem B. Gündüz, Sarath Kumar Murugan and Fabian Muller (2026).
-#>   methylTFR: Quantification of DNA Methylation Signatures in TFBS. R
-#>   package version 0.99.0.
-#>   https://github.com/EpigenomeInformatics/methylTFR
-#> 
-#> A BibTeX entry for LaTeX users is
-#> 
-#>   @Manual{,
-#>     title = {methylTFR: Quantification of DNA Methylation Signatures in TFBS},
-#>     author = {Irem B. Gündüz and Sarath Kumar Murugan and Fabian Muller},
-#>     year = {2026},
-#>     note = {R package version 0.99.0},
-#>     url = {https://github.com/EpigenomeInformatics/methylTFR},
-#>   }
+```bibtex
+@Manual{methylTFR,
+  title = {methylTFR: Quantification of DNA methylation signatures in TFBS},
+  author = {Irem B. Gunduz and Sarath Kumar Murugan and Fabian Mueller},
+  year = {2026},
+  note = {R package version 0.99.8},
+  url = {https://github.com/EpigenomeInformatics/methylTFR},
+  doi = {10.18129/B9.bioc.methylTFR},
+}
 ```
+
+Run `citation("methylTFR")` in R to get the entry for the version you have installed.
diff --git a/inst/CITATION b/inst/CITATION
index 7a6e064..1c40421 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -1,20 +1,33 @@
-pkgVer <- function(pkg) {
-    if (!exists("meta") || is.null(meta)) meta <- packageDescription(pkg)
-    ver <- meta$Version
-    paste0('https://github.com/EpigenomeInformatics/methylTFR', ' - R package version ', ver)
+if (!exists("meta") || is.null(meta)) {
+    meta <- packageDescription("methylTFR")
 }
 
-c(
-    bibentry(bibtype = "Manual",
-        title = "methylTFR: Quantification of DNA Methylation Signatures in TFBS",
-        author = personList(
-            as.person("Irem B. Gunduz"),
-            as.person("Sarath Kumar Murugan"),
-            as.person("Fabian Muller")
-        ),
-        year = format(Sys.Date(), "%Y"),
-        url = "https://github.com/EpigenomeInformatics/methylTFR",
-        note = pkgVer('methylTFR'),
-        doi = "10.18129/B9.bioc.methylTFR"
+pkg_year <- if (!is.null(meta$Date)) {
+    substr(meta$Date, 1, 4)
+} else if (!is.null(meta$`Date/Publication`)) {
+    substr(meta$`Date/Publication`, 1, 4)
+} else {
+    format(Sys.Date(), "%Y")
+}
+pkg_note <- sprintf("R package version %s", meta$Version)
+pkg_title <- "methylTFR: Quantification of DNA methylation signatures in TFBS"
+pkg_url <- "https://github.com/EpigenomeInformatics/methylTFR"
+
+bibentry(
+    bibtype = "Manual",
+    header = "To cite methylTFR in publications use:",
+    title = pkg_title,
+    author = c(
+        person("Irem B.", "Gunduz"),
+        person("Sarath Kumar", "Murugan"),
+        person("Fabian", "Mueller")
+    ),
+    year = pkg_year,
+    note = pkg_note,
+    url = pkg_url,
+    doi = "10.18129/B9.bioc.methylTFR",
+    textVersion = paste0(
+        "Gunduz IB, Murugan SK, Mueller F (", pkg_year, "). ",
+        pkg_title, ". ", pkg_note, ". ", pkg_url
     )
-)
\ No newline at end of file
+)
diff --git a/inst/extdata/BATF_gcfreqs.rda b/inst/extdata/BATF_gcfreqs.rda
index a88e5e2..a06a693 100644
Binary files a/inst/extdata/BATF_gcfreqs.rda and b/inst/extdata/BATF_gcfreqs.rda differ
diff --git a/inst/extdata/BATF_tf_bindsites.rda b/inst/extdata/BATF_tf_bindsites.rda
index 2a66b9f..6d96cce 100644
Binary files a/inst/extdata/BATF_tf_bindsites.rda and b/inst/extdata/BATF_tf_bindsites.rda differ
diff --git a/inst/extdata/encode.tsv.gz b/inst/extdata/encode.tsv.gz
index 95a0573..7e176cc 100644
Binary files a/inst/extdata/encode.tsv.gz and b/inst/extdata/encode.tsv.gz differ
diff --git a/inst/extdata/example_data.rda b/inst/extdata/example_data.rda
index e0c4ef2..68c96b7 100644
Binary files a/inst/extdata/example_data.rda and b/inst/extdata/example_data.rda differ
diff --git a/inst/extdata/gcdist_BATF.rda b/inst/extdata/gcdist_BATF.rda
index b980b75..3523550 100644
Binary files a/inst/extdata/gcdist_BATF.rda and b/inst/extdata/gcdist_BATF.rda differ
diff --git a/inst/extdata/gcdist_subset.rda b/inst/extdata/gcdist_subset.rda
index 4fce124..14f567c 100644
Binary files a/inst/extdata/gcdist_subset.rda and b/inst/extdata/gcdist_subset.rda differ
diff --git a/inst/extdata/msites_sub.rda b/inst/extdata/msites_sub.rda
index 49def96..852e575 100644
Binary files a/inst/extdata/msites_sub.rda and b/inst/extdata/msites_sub.rda differ
diff --git a/inst/extdata/tc_mem.rda b/inst/extdata/tc_mem.rda
index 1fcb510..d75d363 100644
Binary files a/inst/extdata/tc_mem.rda and b/inst/extdata/tc_mem.rda differ
diff --git a/inst/extdata/tc_naive.rda b/inst/extdata/tc_naive.rda
index bb4af32..24446ae 100644
Binary files a/inst/extdata/tc_naive.rda and b/inst/extdata/tc_naive.rda differ
diff --git a/inst/scripts/read_methylome_data.R b/inst/scripts/read_methylome_data.R
index 40954aa..1add6ab 100644
--- a/inst/scripts/read_methylome_data.R
+++ b/inst/scripts/read_methylome_data.R
@@ -1,9 +1,16 @@
 # Create simulated data formats for the methylTFR tests
 output_dir <- "/scratch/icbb/igunduz/methylTFR_manuscript/github/methylTFR/inst/extdata/"
 
-# Encode data
+# Encode data, bedMethyl (bed9+2). Columns 10 and 11 are read coverage and the
+# percentage of reads methylated; the earlier version of this file named them
+# blockCount and blockSizes after the generic BED12 fields.
+#
+# The percentages deliberately span 0 to 100 and are paired with low as well as
+# high coverage. Dividing column 11 by column 10, which is what parse_encode()
+# used to do, returns 20 for the fourth row, so this example fails against that
+# bug instead of passing by luck.
 Encode <- data.frame(
-  chr = c("chr1", "chr1", "chr1", "chr1", "chr1", "chr1"),
+  chrom = rep("chr1", 6),
   start = c(1000170, 1000190, 1000191, 1000198, 1000199, 1000206),
   end = c(1000171, 1000191, 1000192, 1000199, 1000200, 1000207),
   name = rep("HepG2_B1__GC_", 6),
@@ -12,8 +19,8 @@ Encode <- data.frame(
   thickStart = c(1000170, 1000190, 1000191, 1000198, 1000199, 1000206),
   thickEnd = c(1000171, 1000191, 1000192, 1000199, 1000200, 1000207),
   itemRgb = rep("255,255,0", 6),
-  blockCount = c(62, 62, 31, 62, 31, 31),
-  blockSizes = c(6, 3, 0, 10, 0, 10)
+  coverage = c(62, 62, 31, 5, 31, 10),
+  percentMeth = c(6, 3, 0, 100, 55, 100)
 )
 
 # bismarkCytosine data
diff --git a/inst/scripts/vignette_data.R b/inst/scripts/vignette_data.R
index 40043fc..cd36936 100644
--- a/inst/scripts/vignette_data.R
+++ b/inst/scripts/vignette_data.R
@@ -14,7 +14,7 @@ suppressPackageStartupMessages({
 
 # Load the annotation
 gcfreqs <-  getGCfreq(motifSet = "jaspar2020_distal")
-gcdist <- getGenomeGC("hg38")
+gcdist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 # Read the data as GRanges object
@@ -54,8 +54,8 @@ tools::resaveRdaFiles("/scratch/icbb/igunduz/methylTFR_manuscript/github/methylT
 ############################################################################
 # Load the full dataset
 msites <- read_methylome(filename,"bissnp",1)
-enhancer <- readRDS("/icbb/projects/share/annotations/methylTFRAnnotationHg38/inst/extdata/distal_regions.RDS")
-gcdist <- getGenomeGC("hg38")
+enhancer <- readRDS("/scratch/icbb/igunduz/methylTFR_manuscript/github/methylTFRAnnotationHg38_old/inst/extdata/distal_regions.RDS")
+gcdist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 # Check the overlaps with gc_freqs
diff --git a/man/run_methylTFR_RnBeads.Rd b/man/run_methylTFR_RnBeads.Rd
index 69e1f51..47dbd15 100644
--- a/man/run_methylTFR_RnBeads.Rd
+++ b/man/run_methylTFR_RnBeads.Rd
@@ -67,6 +67,12 @@ load(system.file("extdata", "gcdist_subset.rda", package = "methylTFR"))
 
 if (requireNamespace("RnBeads", quietly = TRUE) &&
     requireNamespace("RnBeads.hg38", quietly = TRUE)) {
+    # identifiers.column is a session option, not a property of the object,
+    # so it has to be set before the set is built or samples() falls back
+    # to row numbers and methylTFR cannot name the columns it returns
+    old_ids <- RnBeads::rnb.getOption("identifiers.column")
+    RnBeads::rnb.options(identifiers.column = "sampleName")
+
     sites <- data.frame(
         chr = as.character(GenomicRanges::seqnames(msites)),
         start = GenomicRanges::start(msites),
@@ -76,8 +82,10 @@ if (requireNamespace("RnBeads", quietly = TRUE) &&
     rnb_set <- RnBeads::RnBiseqSet(
         pheno = data.frame(sampleName = "sample_1"),
         sites = sites,
-        meth = matrix(msites$score, ncol = 1),
-        covg = matrix(msites$coverage, ncol = 1),
+        meth = matrix(msites$score, ncol = 1,
+                      dimnames = list(NULL, "sample_1")),
+        covg = matrix(msites$coverage, ncol = 1,
+                      dimnames = list(NULL, "sample_1")),
         assembly = "hg38",
         summarize.regions = FALSE
     )
@@ -88,6 +96,7 @@ if (requireNamespace("RnBeads", quietly = TRUE) &&
         gcfreqs = gcfreqs,
         gc_dist = gcdist
     )
+    RnBeads::rnb.options(identifiers.column = old_ids)
     deviations(devs)
 }
 }
diff --git a/tests/testthat/test_read_methylome.R b/tests/testthat/test_read_methylome.R
index 73ddb50..6384baf 100644
--- a/tests/testthat/test_read_methylome.R
+++ b/tests/testthat/test_read_methylome.R
@@ -44,4 +44,33 @@ test_that("read_methylome", {
 
     # Check the class of the encode object
     expect_s4_class(encode, "GRanges")
+
+    # Every row of the example is kept at the default coverage threshold
+    expect_equal(length(encode), 6)
+
+    # Column 11 is a percentage, so the score is that percentage over 100. Row
+    # four is 100 percent at coverage 5: dividing column 11 by column 10, as
+    # parse_encode() used to, returned 20 for it.
+    expect_equal(encode$score, c(0.06, 0.03, 0, 1, 0.55, 1))
+    expect_equal(encode$coverage, c(62, 62, 31, 5, 31, 10))
+    expect_true(all(encode$score >= 0 & encode$score <= 1))
+
+    # Coverage filtering reads column 10
+    encode_cov <- read_methylome(encode_path, "encode", cov_threshold = 20)
+    expect_equal(length(encode_cov), 4)
+})
+
+test_that("read_methylome rejects an encode file whose column 11 is not a percentage", {
+    bad <- tempfile(fileext = ".tsv")
+    on.exit(unlink(bad), add = TRUE)
+    write.table(
+        data.frame(
+            chrom = "chr1", start = 1000170, end = 1000171,
+            name = "x", score = 62, strand = "+",
+            thickStart = 1000170, thickEnd = 1000171, itemRgb = "255,255,0",
+            coverage = 5, percentMeth = 400
+        ),
+        bad, sep = "\t", row.names = FALSE, quote = FALSE
+    )
+    expect_error(read_methylome(bad, "encode"), "not a\\s+methylation percentage")
 })
diff --git a/vignettes/methylTFR.Rmd b/vignettes/methylTFR.Rmd
index d6c1e83..3075a3c 100644
--- a/vignettes/methylTFR.Rmd
+++ b/vignettes/methylTFR.Rmd
@@ -186,7 +186,7 @@ your working directory.
 library(methylTFRAnnotationHg38) # annotation package for hg38
 
 gcfreqs <- getGCfreq(motifSet = "jaspar2020")
-gc_dist <- getGenomeGC("hg38")
+gc_dist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 sample_dir <- file.path("samples_dir")
@@ -226,7 +226,7 @@ library(methylTFRAnnotationHg38)
 rnb_set <- load.rnb.set("reports/data_import_data/rnb.set_preprocessed")
 
 gcfreqs <- getGCfreq(motifSet = "jaspar2020")
-gc_dist <- getGenomeGC("hg38")
+gc_dist <- getGenomeGC()
 tf_bindsites <- getTFbindsites(motifSet = "jaspar2020")
 
 deviations <- run_methylTFR_RnBeads(