DiffBind consensus peak extraction
0
0
Entering edit mode
6 months ago
wei.zhang • 0

Hi all,

I am trying to extract peaks for each treatment, but met an error I don't know the reason, hope someone can help. The code is:

peaks.consensus <- dba.peakset(Epigenome2,bRetrieve = T) peaks.consensus GRanges object with 59140 ranges and 12 metadata columns: seqnames ranges strand | GNP1_1 GNP1_2 GNP1_3 GNP1_4 GNP2_1 GNP2_2 <Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> 1 chr1 2019459-2020459 | 63.8198 37.5433 45.4025 43.6208 57.3750 79.4296 2 chr1 2443568-2444568 | 44.6738 35.5139 24.6753 28.7501 30.1462 50.6363 3 chr1 2602018-2603018 | 53.1831 17.2496 36.5194 26.7673 28.2013 30.7790 4 chr1 2625942-2626942 | 25.5279 37.5433 43.4285 41.6380 34.0360 48.6506 5 chr1 2628332-2629332 | 87.2203 64.9397 48.3635 63.4484 78.7691 79.4296 ... ... ... ... . ... ... ... ... ... ... 59136 chrX 155828836-155829836 | 47.8648 42.6167 42.4415 44.6122 39.8708 60.5650 59137 chrX 155874227-155876157 | 104.2389 125.8206 89.8180 95.1726 111.8326 121.1301 59138 chrY 6910864-6911864 | 62.7561 36.5286 27.6363 17.8449 37.9259 36.7362 59139 chrY 11106918-11107918 | 67.0107 33.4845 35.5324 29.7414 61.2648 35.7433 59140 chrY 11110491-11111491 | 51.0558 29.4258 40.4674 45.6035 57.3750 53.6150 GNP2_3 GNP2_4 untreated_1 untreated_2 untreated_3 untreated_4 <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> 1 40.3346 45.4483 46.3299 52.4496 40.4039 43.3917 2 31.4807 26.6762 41.2941 27.2334 26.2626 28.5991 3 28.5294 31.6162 46.3299 27.2334 28.2827 24.6544 4 36.3996 32.6042 31.2223 40.3458 28.2827 40.4332 5 79.6855 51.3763 73.5236 73.6311 59.5958 63.1152 ... ... ... ... ... ... ... 59136 33.4482 42.4842 49.3515 57.4928 54.5453 34.5161 59137 110.1825 88.9205 118.8464 118.0115 91.9189 115.3824 59138 31.4807 25.6881 32.2295 28.2421 25.2525 24.6544 59139 28.5294 19.7601 37.2654 42.3631 32.3231 36.4885

59140 50.1724 29.6402 36.2582 52.4496 50.5049 33.5299

seqinfo: 35 sequences from an unspecified genome; no seqlengths

peaks.GNP1<-peaks.consensus[Epigenome2$called[,1]==1 & Epigenome2$called[,2]==1 & Epigenome2$called[,3]==1 & Epigenome2$called[,4]==1] Error: subscript is a logical vector with out-of-bounds TRUE values

DiffBind ChIP-seq • 359 views
ADD COMMENT
0
Entering edit mode

Can you tell us which version(s) you are using? sessionInfo()

ADD REPLY
0
Entering edit mode

Hi Rory, here is the information

sessionInfo() R version 4.0.5 (2021-03-31) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 20.04.2 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale: [1] LC_CTYPE=en_IE.UTF-8 LC_NUMERIC=C LC_TIME=en_IE.UTF-8 LC_COLLATE=en_IE.UTF-8 LC_MONETARY=en_IE.UTF-8
[6] LC_MESSAGES=en_IE.UTF-8 LC_PAPER=en_IE.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_IE.UTF-8 LC_IDENTIFICATION=C

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

other attached packages: [1] TxDb.Hsapiens.UCSC.hg38.knownGene_3.10.0 GenomicFeatures_1.42.3 AnnotationDbi_1.52.0
[4] biomaRt_2.46.3 clusterProfiler_3.18.1 ReactomePA_1.34.0
[7] ChIPseeker_1.26.2 DiffBind_3.0.15 SummarizedExperiment_1.20.0
[10] Biobase_2.50.0 MatrixGenerics_1.2.1 matrixStats_0.58.0
[13] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1
[16] S4Vectors_0.28.1 BiocGenerics_0.36.0

loaded via a namespace (and not attached): [1] utf8_1.2.1 tidyselect_1.1.0 RSQLite_2.2.6
[4] grid_4.0.5 BiocParallel_1.24.1 scatterpie_0.1.5
[7] munsell_0.5.0 base64url_1.4 systemPipeR_1.24.3
[10] withr_2.4.1 colorspace_2.0-0 GOSemSim_2.16.1
[13] Category_2.56.0 rstudioapi_0.13 DOSE_3.16.0
[16] bbmle_1.0.23.1 GenomeInfoDbData_1.2.4 mixsqp_0.3-43
[19] hwriter_1.3.2 polyclip_1.10-0 bit64_4.0.5
[22] farver_2.1.0 pheatmap_1.0.12 downloader_0.4
[25] coda_0.19-4 batchtools_0.9.15 vctrs_0.3.7
[28] generics_0.1.0 BiocFileCache_1.14.0 R6_2.5.0
[31] apeglm_1.12.0 graphlayouts_0.7.1 invgamma_1.1
[34] locfit_1.5-9.4 rsvg_2.1 bitops_1.0-6
[37] cachem_1.0.4 fgsea_1.16.0 DelayedArray_0.16.3
[40] assertthat_0.2.1 scales_1.1.1 ggraph_2.0.5
[43] enrichplot_1.10.2 gtable_0.3.0 tidygraph_1.2.0
[46] rlang_0.4.10 genefilter_1.72.1 splines_4.0.5
[49] rtracklayer_1.50.0 brew_1.0-6 checkmate_2.0.0
[52] BiocManager_1.30.12 yaml_2.2.1 reshape2_1.4.4
[55] backports_1.2.1 qvalue_2.22.0 RBGL_1.66.0
[58] tools_4.0.5 ggplot2_3.3.3 ellipsis_0.3.1
[61] gplots_3.1.1 RColorBrewer_1.1-2 Rcpp_1.0.6
[64] plyr_1.8.6 progress_1.2.2 zlibbioc_1.36.0
[67] purrr_0.3.4 RCurl_1.98-1.3 prettyunits_1.1.1
[70] openssl_1.4.3 viridis_0.6.0 ashr_2.2-47
[73] cowplot_1.1.1 ggrepel_0.9.1 magrittr_2.0.1
[76] data.table_1.14.0 DO.db_2.9 reactome.db_1.74.0
[79] truncnorm_1.0-8 mvtnorm_1.1-1 SQUAREM_2021.1
[82] amap_0.8-18 hms_1.0.0 xtable_1.8-4
[85] XML_3.99-0.6 emdbook_1.3.12 jpeg_0.1-8.1
[88] gridExtra_2.3 compiler_4.0.5 bdsmatrix_1.3-4
[91] tibble_3.1.0 KernSmooth_2.23-18 V8_3.4.0
[94] crayon_1.4.1 shadowtext_0.0.7 GOstats_2.56.0
[97] tidyr_1.1.3 geneplotter_1.68.0 DBI_1.1.1
[100] tweenr_1.0.2 dbplyr_2.1.1 MASS_7.3-53.1
[103] rappdirs_0.3.3 boot_1.3-27 ShortRead_1.48.0
[106] Matrix_1.3-2 igraph_1.2.6 pkgconfig_2.0.3
[109] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 rvcheck_0.1.8 GenomicAlignments_1.26.0
[112] numDeriv_2016.8-1.1 xml2_1.3.2 annotate_1.68.0
[115] XVector_0.30.0 AnnotationForge_1.32.0 stringr_1.4.0
[118] VariantAnnotation_1.36.0 digest_0.6.27 graph_1.68.0
[121] Biostrings_2.58.0 fastmatch_1.1-0 edgeR_3.32.1
[124] GSEABase_1.52.1 GreyListChIP_1.22.0 curl_4.3
[127] graphite_1.36.0 Rsamtools_2.6.0 gtools_3.8.2
[130] rjson_0.2.20 lifecycle_1.0.0 jsonlite_1.7.2
[133] viridisLite_0.4.0 askpass_1.1 limma_3.46.0
[136] BSgenome_1.58.0 fansi_0.4.2 pillar_1.6.0
[139] lattice_0.20-41 fastmap_1.1.0 httr_1.4.2
[142] plotrix_3.8-1 survival_3.2-10 GO.db_3.12.1
[145] glue_1.4.2 png_0.1-7 bit_4.0.4
[148] Rgraphviz_2.34.0 ggforce_0.3.3 stringi_1.5.3
[151] blob_1.2.1 DESeq2_1.30.1 latticeExtra_0.6-29
[154] caTools_1.18.2 memoise_2.0.0 DOT_0.1
[157] dplyr_1.0.5 irlba_2.3.3

ADD REPLY
0
Entering edit mode

I'd need to see the sequence of calls that led to this: the call to dba(), dba.count(), and any intervening calls.

If you could send me a link to the DBA object before and after counting I can have a look at what is going on.

ADD REPLY

Login before adding your answer.

Traffic: 1343 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6