Hi,
I'm importing fragments, peaks, and counts from chromap and trying to get all the metadata collected using signac functions.
For each sample, I'm doing the following
# Load genome annotations
annotations = GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) = 'UCSC'
genome(annotations) = "hg38"
# Load peaks as GRanges object
p = as.data.frame(read.table("peaks_merged.bed",header=F,sep="\t"))
colnames(p) = c("chr","start","stop")
peaks = makeGRangesFromDataFrame(p)
# Create Fragment objects, then Feature Matrix, then Seurat objects for each sample
p1_path <- "patient1_fragments.bed.gz"
p1_cells <- read_tsv(p1_path,col_names=c("chr","start","stop","cell","support"),col_types=c("-","-","-","-","c"),col_select="cell") %>% pull(cell) %>% unique()
names(x = p1_cells) <- paste0("Patient1-diag", p1_cells)
p1_frags <- CreateFragmentObject(path = p1_path, cells = p1_cells, max.lines=NULL)
p1_mat <- FeatureMatrix(
fragments = p1_frags,
features = peaks,
process_n = 20000,
sep = c("-", "-"),
verbose = TRUE
)
p1_assay <- CreateChromatinAssay(p1_mat, fragments = p1_frags, genome = 'hg38')
p1_seurat <- CreateSeuratObject(p1_assay, assay = "peaks")
p1_seurat$Sample <- "Patient1-diag"
p1_seurat$Patient <- "Patient1"
p1_seurat$Disease <- "Diag"
This all goes fine.
However if I try to compute Nucleosomal signal with NucleosomeSignal, I get the following error:
p1_seurat <- NucleosomeSignal(p1_seurat)
Found 3345708 cell barcodes
Done Processing 150 million linesError in ecdf(x = af$nucleosome_signal) :
'x' must have 1 or more non-missing values
I thought this maybe was related to #786 or #826, but my annotations/seqlevels all appear to be correct as is:
> seqlevelsStyle(annotations)
[1] "UCSC"
> seqlevels(annotations)
[1] "chrX" "chr20" "chr1" "chr6" "chr3" "chr7" "chr12" "chr11" "chr4" "chr17" "chr2"
[12] "chr16" "chr8" "chr19" "chr9" "chr13" "chr14" "chr5" "chr22" "chr10" "chrY" "chr18"
[23] "chr15" "chr21" "chrM"
> head(Fragments(p1_seurat)[[1]])
chrom start end barcode readCount
1 chr1 10084 10497 TACTGCCAGCTAACAA 1
2 chr1 10091 10301 CACCACTCAGAAAGCG 1
3 chr1 10091 10301 CACCACTCAGTGAGCG 1
4 chr1 10095 10507 GTGTCAAGTGATGCGA 1
5 chr1 10097 10512 AGTCAACCATAAAGTG 1
6 chr1 10151 10186 TAGCTTTGTTTGATCG 1
I also suspected it was because my seurat object didn't have CountFragments run on it, which I can do separately (successfully):
p1_fragmentInfo <- CountFragments(p1_path)
rownames(p1_fragmentInfo) <- p1_fragmentInfo$CB
p1_seurat$fragments <- p1_fragmentInfo[colnames(p1_seurat), "frequency_count"]
p1_seurat$mononucleosomal <- p1_fragmentInfo[colnames(p1_seurat), "mononucleosomal"]
p1_seurat$nucleosome_free <- p1_fragmentInfo[colnames(p1_seurat), "nucleosome_free"]
p1_seurat$reads_count <- p1_fragmentInfo[colnames(p1_seurat), "reads_count"]
p1_seurat <- FRiP(
object = p1_seurat,
assay = 'peaks',
total.fragments = "fragments"
)
p1_seurat$blacklist_fraction <- FractionCountsInRegion(
object = p1_seurat,
assay = 'peaks',
regions = blacklist_hg38
)
both FRiP and FractionCountsInRegion run successfully, however NucleosomeSignal and TSSEnrichment continue to fail with the same error.
I feel like there must be something these functions is expecting to find that isn't there, but I can't figure out what.
Any help would be appreciated!
sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.5 (Ootpa)
Matrix products: default
BLAS/LAPACK: /nas/longleaf/rhel8/apps/r/4.1.0/lib/libopenblas_haswellp-r0.3.5.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.18.3 AnnotationFilter_1.18.0
[4] GenomicFeatures_1.46.5 AnnotationDbi_1.56.2 Biobase_2.54.0
[7] Signac_1.7.0 SeuratObject_4.0.4 Seurat_4.1.0
[10] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 IRanges_2.28.0
[13] S4Vectors_0.32.4 BiocGenerics_0.40.0 forcats_0.5.1
[16] stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
[19] readr_2.1.2 tidyr_1.2.0 tibble_3.1.8
[22] ggplot2_3.3.6 tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] utf8_1.2.2 reticulate_1.25 tidyselect_1.1.2
[4] RSQLite_2.2.10 htmlwidgets_1.5.4 grid_4.1.0
[7] BiocParallel_1.28.3 Rtsne_0.15 munsell_0.5.0
[10] codetools_0.2-18 ica_1.0-2 future_1.24.0
[13] miniUI_0.1.1.1 withr_2.5.0 spatstat.random_2.1-0
[16] colorspace_2.0-3 filelock_1.0.2 rstudioapi_0.13
[19] ROCR_1.0-11 tensor_1.5 listenv_0.8.0
[22] MatrixGenerics_1.6.0 GenomeInfoDbData_1.2.7 polyclip_1.10-0
[25] bit64_4.0.5 parallelly_1.30.0 vctrs_0.4.1
[28] generics_0.1.2 BiocFileCache_2.2.1 R6_2.5.1
[31] DelayedArray_0.20.0 bitops_1.0-7 spatstat.utils_2.3-0
[34] cachem_1.0.6 assertthat_0.2.1 promises_1.2.0.1
[37] BiocIO_1.4.0 scales_1.2.0 gtable_0.3.0
[40] globals_0.14.0 goftest_1.2-3 rlang_1.0.4
[43] RcppRoll_0.3.0 splines_4.1.0 rtracklayer_1.54.0
[46] lazyeval_0.2.2 spatstat.geom_2.3-2 broom_1.0.0
[49] yaml_2.3.5 reshape2_1.4.4 abind_1.4-5
[52] modelr_0.1.8 backports_1.4.1 httpuv_1.6.5
[55] tools_4.1.0 ellipsis_0.3.2 spatstat.core_2.4-0
[58] RColorBrewer_1.1-3 ggridges_0.5.3 Rcpp_1.0.8.3
[61] plyr_1.8.7 progress_1.2.2 zlibbioc_1.40.0
[64] RCurl_1.98-1.6 prettyunits_1.1.1 rpart_4.1.16
[67] deldir_1.0-6 pbapply_1.5-0 cowplot_1.1.1
[70] zoo_1.8-9 SummarizedExperiment_1.24.0 haven_2.4.3
[73] ggrepel_0.9.1 cluster_2.1.2 fs_1.5.2
[76] magrittr_2.0.2 data.table_1.14.2 scattermore_0.8
[79] lmtest_0.9-40 reprex_2.0.1 RANN_2.6.1
[82] ProtGenerics_1.26.0 fitdistrplus_1.1-6 matrixStats_0.62.0
[85] hms_1.1.1 patchwork_1.1.1 mime_0.12
[88] xtable_1.8-4 XML_3.99-0.9 readxl_1.3.1
[91] gridExtra_2.3 compiler_4.1.0 biomaRt_2.50.3
[94] KernSmooth_2.23-20 crayon_1.5.1 htmltools_0.5.2
[97] mgcv_1.8-40 later_1.3.0 tzdb_0.2.0
[100] lubridate_1.8.0 DBI_1.1.2 dbplyr_2.1.1
[103] MASS_7.3-55 rappdirs_0.3.3 Matrix_1.4-0
[106] cli_3.3.0 parallel_4.1.0 igraph_1.3.3
[109] pkgconfig_2.0.3 GenomicAlignments_1.30.0 plotly_4.10.0
[112] spatstat.sparse_2.1-0 xml2_1.3.3 XVector_0.34.0
[115] rvest_1.0.2 digest_0.6.29 sctransform_0.3.3
[118] RcppAnnoy_0.0.19 spatstat.data_2.1-2 Biostrings_2.62.0
[121] cellranger_1.1.0 leiden_0.3.9 fastmatch_1.1-3
[124] uwot_0.1.11 restfulr_0.0.13 curl_4.3.2
[127] shiny_1.7.1 Rsamtools_2.10.0 rjson_0.2.21
[130] lifecycle_1.0.1 nlme_3.1-155 jsonlite_1.8.0
[133] viridisLite_0.4.0 fansi_1.0.3 pillar_1.7.0
[136] lattice_0.20-45 KEGGREST_1.34.0 fastmap_1.1.0
[139] httr_1.4.2 survival_3.2-13 glue_1.6.2
[142] png_0.1-7 bit_4.0.4 stringi_1.7.6
[145] blob_1.2.2 memoise_2.0.1 irlba_2.3.5
[148] future.apply_1.8.1
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4