Closed bschilder closed 2 years ago
bschilder
commented
2 years ago Ah, actually I see what's happening. reference= used to take a single GRanges object, but it now takes a named list only.
Might be good to add a check at the beginning to see if the supplied reference is actually a GRanges object, and if so turn it into a list (with a single element) and give it a name based on the object name.
bschilder
commented
2 years ago Ok, got much further than before by supplying reference as a list, but still hit another error. Is this because reference can only handle a single element?
EpiCompare::EpiCompare(peakfiles = peakfiles,
picard_files = picardfiles,
blacklist = hg19_blacklist,
genome_build = "hg38",
reference = ref_list,
stat_plot = TRUE,
chromHMM_plot = TRUE,
chipseeker_plot = TRUE,
enrichment_plot = TRUE,
save_output = TRUE,
chromHMM_annotation = "K562",
output_dir = save_dir <- here::here("reports","encode_ac.encode_me3"))Saving 10 x 14 in image
Saving 10 x 14 in image
Rows: 2 Columns: 9
── Column specification ─────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): X1, X4, X6
dbl (5): X2, X3, X5, X7, X8
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
/Users/schilder/Library/Caches/org.R-project.R/R/BiocFileCache
does not exist, create directory? (yes/no): yes
adding K562's chrHMM to local cache,future invocations will use local image
adding rname 'http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmK562HMM.bed.gz'
|===============================================================================| 100%
Rows: 2 Columns: 9
── Column specification ─────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): X1, X4, X6
dbl (5): X2, X3, X5, X7, X8
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
|===============================================================================| 100%
snapshotDate(): 2022-04-21
downloading 1 resources
retrieving 1 resource
|===============================================================================| 100%
loading from cache
require("rtracklayer")
Working on: phase_1_05_jan_2022.S_1_R1
Working on: phase_1_05_jan_2022.S_1_R2
Working on: phase_1_05_jan_2022.S_1_R1
Working on: phase_1_05_jan_2022.S_1_R2
Working on: phase_2_03_feb_2022.S_4_R1
Working on: phase_2_03_feb_2022.S_4_R1
Working on: phase_2_28_jan_2022.S_2_R1
Working on: phase_2_28_jan_2022.S_3_R1
Working on: phase_2_28_jan_2022.S_4_R1
Working on: phase_2_28_jan_2022.S_5_R1
Working on: phase_2_28_jan_2022.S_6_R1
Working on: phase_2_28_jan_2022.S_2_R1
Working on: phase_2_28_jan_2022.S_3_R1
Working on: phase_2_28_jan_2022.S_4_R1
Working on: phase_2_28_jan_2022.S_5_R1
Working on: phase_2_28_jan_2022.S_6_R1
Working on: phase_3_18_feb_2022.S_10_R1
Working on: phase_3_18_feb_2022.S_11_R1
Working on: phase_3_18_feb_2022.S_8_R1
Working on: phase_3_18_feb_2022.S_9_R1
Working on: phase_3_18_feb_2022.S_10_R1
Working on: phase_3_18_feb_2022.S_11_R1
Working on: phase_3_18_feb_2022.S_8_R1
Working on: phase_3_18_feb_2022.S_9_R1
Working on: phase_4_24_mar_2022.s1_R1
Working on: phase_4_24_mar_2022.s11_R1
Working on: phase_4_24_mar_2022.s4_R1
Working on: phase_4_24_mar_2022.s6_R1
Working on: phase_4_24_mar_2022.s8_R1
Working on: phase_4_24_mar_2022.s9_R1
Working on: phase_4_24_mar_2022.s1_R1
Working on: phase_4_24_mar_2022.s11_R1
Working on: phase_4_24_mar_2022.s4_R1
Working on: phase_4_24_mar_2022.s6_R1
Working on: phase_4_24_mar_2022.s8_R1
Working on: phase_4_24_mar_2022.s9_R1
Working on: encode_ac
Working on: encode_me3
Quitting from lines 333-358 (EpiCompare.Rmd)
Error in dimnames(x) <- dn :
length of 'dimnames' [1] not equal to array extentHm, same issue as before when reference is named list with only a single element. So not sure why I'm getting this dimnames error.
Same thing occurs when i remove picardfiles as an arg, so that doesn't seem to be it either.
Based on the latest error, the source of the bug is somewhere in here: https://github.com/neurogenomics/EpiCompare/blob/2b84d5eba23aec7ca6c5db1bfbd27e4d1606dc29/inst/markdown/EpiCompare.Rmd#L332
serachoi1230
commented
2 years ago I've seen that error before, it's very strange, it comes and goes. Sometimes it runs all checks fine without errors but when I try run it, it gives me that error. Will have a look!
bschilder
commented
2 years ago Here is the full trackback, which might be helpful.
Seems the source of the issue is here: https://github.com/neurogenomics/EpiCompare/blob/3f59ce8e736d75ab3c0768fa617237dfffb8e696/R/overlap_stat_plot.R#L53
encode_ac
Saving 10 x 13.7 in image
Quitting from lines 296-311 (EpiCompare.Rmd)
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'x' in selecting a method for function 'elementMetadata': GRanges objects don't support [[, as.list(), lapply(), or unlist() at the
moment
35.
h(simpleError(msg, call))
34.
.handleSimpleError(function (cond)
.Internal(C_tryCatchHelper(addr, 1L, cond)), "GRanges objects don't support [[, as.list(), lapply(), or unlist() at the\n moment",
base::quote(getListElement(x, i, ...)))
33.
stop(wmsg(class(x), " objects don't support [[, as.list(), ",
"lapply(), or unlist() at the moment"))
32.
getListElement(x, i, ...)
31.
getListElement(x, i, ...)
30.
reference[[1]]
29.
reference[[1]] at overlap_stat_plot.R#53
28.
S4Vectors::elementMetadata(reference[[1]]) at overlap_stat_plot.R#53
27.
ncol(S4Vectors::elementMetadata(reference[[1]])) at overlap_stat_plot.R#53
26.
overlap_stat_plot(reference = reference_tidy, peaklist = peaklist_tidy,
annotation = txdb) at <text>#4
25.
eval(expr, envir, enclos)
24.
eval(expr, envir, enclos)
23.
eval_with_user_handlers(expr, envir, enclos, user_handlers)
22.
withVisible(eval_with_user_handlers(expr, envir, enclos, user_handlers))
21.
withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler)
20.
handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler))
19.
timing_fn(handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler)))
18.
evaluate_call(expr, parsed$src[[i]], envir = envir, enclos = enclos,
debug = debug, last = i == length(out), use_try = stop_on_error !=
2L, keep_warning = keep_warning, keep_message = keep_message,
output_handler = output_handler, include_timing = include_timing)
17.
evaluate::evaluate(...)
16.
evaluate(code, envir = env, new_device = FALSE, keep_warning = !isFALSE(options$warning),
keep_message = !isFALSE(options$message), stop_on_error = if (is.numeric(options$error)) options$error else {
if (options$error && options$include)
0L ...
15.
in_dir(input_dir(), expr)
14.
in_input_dir(evaluate(code, envir = env, new_device = FALSE,
keep_warning = !isFALSE(options$warning), keep_message = !isFALSE(options$message),
stop_on_error = if (is.numeric(options$error)) options$error else {
if (options$error && options$include) ...
13.
eng_r(options)
12.
block_exec(params)
11.
call_block(x)
10.
process_group.block(group)
9.
process_group(group)
8.
withCallingHandlers(if (tangle) process_tangle(group) else process_group(group),
error = function(e) {
setwd(wd)
cat(res, sep = "\n", file = output %n% "") ...
7.
process_file(text, output)
6.
knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
5.
rmarkdown::render(input = markdown_path, output_dir = output_dir,
output_file = output_filename, quiet = TRUE, params = list(peakfile = peakfiles,
genome_build = genome_build, blacklist = blacklist, picard_files = picard_files,
reference = reference, upset_plot = upset_plot, stat_plot = stat_plot, ... at EpiCompare.R#130
4.
EpiCompare::EpiCompare(peakfiles = peakfiles, picard_files = picardfiles,
blacklist = hg19_blacklist, genome_build = "hg38", reference = ref_list[[nm]],
stat_plot = TRUE, chromHMM_plot = TRUE, chipseeker_plot = TRUE,
enrichment_plot = TRUE, save_output = TRUE, chromHMM_annotation = "K562", ...
3.
FUN(X[[i]], ...)
2.
lapply(names(ref_list), function(nm) {
message("\n", nm)
save_dir <- here::here("reports", nm)
dir.create(save_dir, showWarnings = FALSE, recursive = TRUE) ...
1.
lapply(names(ref_list), function(nm) {
message("\n", nm)
save_dir <- here::here("reports", nm)
dir.create(save_dir, showWarnings = FALSE, recursive = TRUE) ...Ok, so I figured this out. It had to do with the fact that reference[[1]] was being called regardless of whether reference was a list or an unlisted GRanges object. It looks like this was being done because overlap_stat_plot can only use a single reference at a time.
I've added a function that validates the reference format and returns the format required by that particular function: prepare_reference.R
However now I'm getting the original error again. I'm able to replicate the specific source of the error using just my peakfiles.
peaklist <- readRDS("~/Downloads/peakfiles.rds")
my_plot <- plot_chromHMM(peaklist = peaklist,
cell_line = "K562",
genome_build = "hg19") Here's the traceback:
**Working on: encode_ac
Quitting from lines 336-361 (EpiCompare.Rmd)
Error in dimnames(x) <- dn :
length of 'dimnames' [1] not equal to array extent
27.
`rownames<-`(`*tmp*`, value = names(peaklist)) at plot_chromHMM.R#81
26.
plot_chromHMM(peaklist = peaklist_tidy, chromHMM_annotation = chromHMM_list,
genome_build = params$genome_build, interact = params$interact) at <text>#6
25.
eval(expr, envir, enclos)
24.
eval(expr, envir, enclos)
23.
eval_with_user_handlers(expr, envir, enclos, user_handlers)
22.
withVisible(eval_with_user_handlers(expr, envir, enclos, user_handlers))
21.
withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler)
20.
handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler))
19.
timing_fn(handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr,
envir, enclos, user_handlers)), warning = wHandler, error = eHandler,
message = mHandler)))
18.
evaluate_call(expr, parsed$src[[i]], envir = envir, enclos = enclos,
debug = debug, last = i == length(out), use_try = stop_on_error !=
2L, keep_warning = keep_warning, keep_message = keep_message,
output_handler = output_handler, include_timing = include_timing)
17.
evaluate::evaluate(...)
16.
evaluate(code, envir = env, new_device = FALSE, keep_warning = !isFALSE(options$warning),
keep_message = !isFALSE(options$message), stop_on_error = if (is.numeric(options$error)) options$error else {
if (options$error && options$include)
0L ...
15.
in_dir(input_dir(), expr)
14.
in_input_dir(evaluate(code, envir = env, new_device = FALSE,
keep_warning = !isFALSE(options$warning), keep_message = !isFALSE(options$message),
stop_on_error = if (is.numeric(options$error)) options$error else {
if (options$error && options$include) ...
13.
eng_r(options)
12.
block_exec(params)
11.
call_block(x)
10.
process_group.block(group)
9.
process_group(group)
8.
withCallingHandlers(if (tangle) process_tangle(group) else process_group(group),
error = function(e) {
setwd(wd)
cat(res, sep = "\n", file = output %n% "") ...
7.
process_file(text, output)
6.
knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
5.
rmarkdown::render(input = markdown_path, output_dir = output_dir,
output_file = output_filename, quiet = TRUE, params = list(peakfile = peakfiles,
genome_build = genome_build, blacklist = blacklist, picard_files = picard_files,
reference = reference, upset_plot = upset_plot, stat_plot = stat_plot, ... at EpiCompare.R#130
4.
EpiCompare::EpiCompare(peakfiles = peakfiles, picard_files = picardfiles,
blacklist = hg19_blacklist, genome_build = "hg38", reference = ref_list[nm],
stat_plot = TRUE, chromHMM_plot = TRUE, chipseeker_plot = TRUE,
enrichment_plot = TRUE, save_output = TRUE, chromHMM_annotation = "K562", ...
3.
FUN(X[[i]], ...)
2.
lapply(names(ref_list), function(nm) {
message("\n", nm)
save_dir <- here::here("reports", nm)
dir.create(save_dir, showWarnings = FALSE, recursive = TRUE) ...
1.
lapply(names(ref_list), function(nm) {
message("\n", nm)
save_dir <- here::here("reports", nm)
dir.create(save_dir, showWarnings = FALSE, recursive = TRUE) ...**Ok, figured it out.
It was stemming from the fact that some of my peak GRanges in peaklist were not uniquely named.
In plot_chromHMM there's a step where you name the rows of the genomation::heatTargetAnnotation output matrix:
https://github.com/neurogenomics/EpiCompare/blob/dc31e81730b9c79adfc42ae743deb93399f8e09b/R/plot_chromHMM.R#L81
This causes an error because the matrix somehow aggregates items with the same name, producing a matrix with less rows than items in peaklist. But your manual naming step isn't really necessary since the matrix already uses the the peaklist names to name its rows, so I removed it.
In addition, I've added an extra check to make all peaklist names are unique, and if not, it gives the user a message saying they're being forced to be unique.
I've also gone through and simplified some of your code along the way, speeding it up wherever possible with more efficient equivalents.
bschilder
commented
2 years ago On a related note, the duplicate peak names were ultimately stemming from EpiCompare::gather_files.
nf-core/cutandrun creates duplicates of "peaks.stringent" files in different subfolders:
../04_reporting/ and ../03_peak_calling/
I didn't account for this when naming each file. I'll modify gather_files accordingly.
Can't seem to find anything wrong with my input files, but I keep getting this error.
Reprex
I've just uploaded the peak and picard list objects as .rds files to our lab folder.
Error
Session info