Entering edit mode
Hi,
I have started using MEDIPS and have managed to run several of the functions: MEDIPS.saturation, MEDIPS.correlation, and MEDIPS.seqCoverage. However, when I try MEDIPS.CpGenrich I get an error.
This is my code:
PaleMerge1_CpGenrich = MEDIPS.CpGenrich(file="../BAM_files/PaleMerge1_R1_trimmoSE_bowtie241_aphid_v3_Q30_sorted.bam", BSgenome=BSgenome, extend=150, shift=shift, uniq=uniq, paired = FALSE)
This is the output with the error:
Reading bam alignment PaleMerge1_R1_trimmoSE_bowtie241_aphid_v3_Q30_sorted.bam
Total number of imported short reads: 80500164
Extending reads...
Creating GRange Object...
NaNs producedKeep at most NaN read(s) mapping to the same genomic location
Number of remaining reads: 80500164
Loading chromosome lengths for BSgenome.Apisum.NCBI.acyPis3...
Calculating CpG density for given regions...
Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'XStringSet': solving row 9086: 'allow.nonnarrowing' is FALSE and the supplied start (132645921) is > refwidth + 1
Here is the Traceback:
20.
h(simpleError(msg, call))
19.
.handleSimpleError(function (cond) .Internal(C_tryCatchHelper(addr, 1L, cond)), "solving row 9086: 'allow.nonnarrowing' is FALSE and the supplied start (132645921) is > refwidth + 1", base::quote(.Call2("C_solve_user_SEW", refwidths, start, end, width, translate.negative.coord, allow.nonnarrowing, ...
18.
.Call2("C_solve_user_SEW", refwidths, start, end, width, translate.negative.coord, allow.nonnarrowing, PACKAGE = "IRanges")
17.
solveUserSEW(refwidths, start = start, end = end, width = width)
16.
.extractFromBSgenomeSingleSequences(x, sseq_args$names, sseq_args$start, sseq_args$end, sseq_args$width, sseq_args$strand)
15.
.local(x, ...)
14.
getSeq(dataset, names = as.vector(seqnames(x)), start = start(x), end = end(x), as.character = TRUE)
13.
getSeq(dataset, names = as.vector(seqnames(x)), start = start(x), end = end(x), as.character = TRUE)
12.
XStringSet("DNA", x, start = start, end = end, width = width, use.names = use.names)
11.
DNAStringSet(getSeq(dataset, names = as.vector(seqnames(x)), start = start(x), end = end(x), as.character = TRUE))
10.
FUN(extractROWS(unlisted_X, IRanges(X_elt_start[i], X_elt_end[i])), ...)
9.
FUN(X[[i]], ...)
8.
lapply(non_empty_idx, function(i) FUN(extractROWS(unlisted_X, IRanges(X_elt_start[i], X_elt_end[i])), ...))
7.
lapply(non_empty_idx, function(i) FUN(extractROWS(unlisted_X, IRanges(X_elt_start[i], X_elt_end[i])), ...))
6.
lapply_CompressedList(X, FUN, ...)
5.
IRanges::lapply(GRangesList(GRange.Reads), function(x) { ranges(x) <- restrict(ranges(x), end = chr_lengths[which(chromosomes %in% as.vector(seqnames(x)))]) y = DNAStringSet(getSeq(dataset, names = as.vector(seqnames(x)), ...
4.
IRanges::lapply(GRangesList(GRange.Reads), function(x) { ranges(x) <- restrict(ranges(x), end = chr_lengths[which(chromosomes %in% as.vector(seqnames(x)))]) y = DNAStringSet(getSeq(dataset, names = as.vector(seqnames(x)), ...
3.
unlist(IRanges::lapply(GRangesList(GRange.Reads), function(x) { ranges(x) <- restrict(ranges(x), end = chr_lengths[which(chromosomes %in% as.vector(seqnames(x)))]) y = DNAStringSet(getSeq(dataset, names = as.vector(seqnames(x)), ...
2.
matrix(unlist(IRanges::lapply(GRangesList(GRange.Reads), function(x) { ranges(x) <- restrict(ranges(x), end = chr_lengths[which(chromosomes %in% as.vector(seqnames(x)))]) y = DNAStringSet(getSeq(dataset, names = as.vector(seqnames(x)), ...
1.
MEDIPS.CpGenrich(file = "../BAM_files/PaleMerge1_R1_trimmoSE_bowtie241_aphid_v3_Q30_sorted.bam", BSgenome = BSgenome, extend = 150, shift = shift, uniq = uniq, paired = FALSE)
The same BAM files are used with the successfully run functions. I made my own BSgenome package using the forge.
During the installation of MEDIPS there following message was shown several times - don't know if this is relevant, but looks odd:
No methods found in package ‘IRanges’ for request: ‘values’ when loading ‘MEDIPS’
Thank you for any help regarding this. Ian
sessionInfo( )
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
[6] LC_MESSAGES=en_GB.UTF-8 LC_PAPER=en_GB.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] BSgenome.Apisum.NCBI.acyPis3_3.0 MEDIPS_1.46.0 Rsamtools_2.10.0 BSgenome_1.62.0
[5] rtracklayer_1.54.0 Biostrings_2.62.0 XVector_0.34.0 GenomicRanges_1.46.0
[9] GenomeInfoDb_1.30.0 IRanges_2.28.0 S4Vectors_0.32.0 BiocGenerics_0.40.0
loaded via a namespace (and not attached):
[1] MatrixGenerics_1.6.0 Biobase_2.54.0 httr_1.4.2 bit64_4.0.5 gtools_3.9.2
[6] assertthat_0.2.1 BiocManager_1.30.16 BiocFileCache_2.2.0 blob_1.2.2 GenomeInfoDbData_1.2.7
[11] yaml_2.2.1 progress_1.2.2 pillar_1.6.4 RSQLite_2.2.8 lattice_0.20-45
[16] glue_1.4.2 digest_0.6.28 htmltools_0.5.2 preprocessCore_1.56.0 Matrix_1.3-4
[21] XML_3.99-0.8 pkgconfig_2.0.3 biomaRt_2.50.0 zlibbioc_1.40.0 purrr_0.3.4
[26] BiocParallel_1.28.0 tibble_3.1.5 KEGGREST_1.34.0 generics_0.1.1 ellipsis_0.3.2
[31] cachem_1.0.6 SummarizedExperiment_1.24.0 magrittr_2.0.1 crayon_1.4.2 memoise_2.0.0
[36] evaluate_0.14 fansi_0.5.0 xml2_1.3.2 tools_4.1.2 prettyunits_1.1.1
[41] hms_1.1.1 BiocIO_1.4.0 lifecycle_1.0.1 matrixStats_0.61.0 stringr_1.4.0
[46] DelayedArray_0.20.0 AnnotationDbi_1.56.1 compiler_4.1.2 rlang_0.4.12 grid_4.1.2
[51] RCurl_1.98-1.5 rstudioapi_0.13 rjson_0.2.20 rappdirs_0.3.3 bitops_1.0-7
[56] rmarkdown_2.11 DNAcopy_1.68.0 restfulr_0.0.13 DBI_1.1.1 curl_4.3.2
[61] R6_2.5.1 GenomicAlignments_1.30.0 knitr_1.36 dplyr_1.0.7 fastmap_1.1.0
[66] bit_4.0.4 utf8_1.2.2 filelock_1.0.2 stringi_1.7.5 parallel_4.1.2
[71] Rcpp_1.0.7 vctrs_0.3.8 png_0.1-7 xfun_0.27 dbplyr_2.1.1
[76] tidyselect_1.1.1
Hi,
I am sorry to hear about that problem with the MEDIPS.CpGenrich function.
The MEDIPS.CpGenrich() function used to depend on RangedData, which is not available any more. I had updated MEDIPS to function without RangedData and this should be available in MEDIPS >= 1.42.0. However, I do not know why it does not work for you anymore.
Based on the error/ warning massage you receive while installing MEDIPS, there might be a problem with the dependencies on IRanges. When I check the build report on https://www.bioconductor.org/checkResults/3.14/bioc-LATEST/MEDIPS/ I don't see such warnings/ errors, though. Have you tried to run the MEDIPS.CpGenrich function on the example data and does that work? Or is this maybe a problem with your specific BSgenome?
I also would like to stress that we have developed a new package ‘qsea’, which has a much improved method for estimating %methylation from MeDIP-seq data. It does not have the CG enrichment QC, but is much more flexible in terms of parameter settings etc.
Thank you and all the best, Lukas
P.S. Please note, if there is a broken dependency of MEDIPS on other packages, I will not have the time to fix this problem in the near future. Thank you for your understanding.
All the best, Lukas
Thank you for your reply. MEDIPS does work with the example sequences, so there is something wrong with my data. I'll start by using one of the biggest contigs, as there are many that are many small ones, the smallest being 200bp.