GeneActivity() return all zero matrix on merged object.

[Ceeyouthere]

Hi,
When I run the GeneActivity() function on a merged object, the returned matrix contains all zero. After running the function step by step, I found that it might have something to do with the cell barcodes.

The original code from GeneActivity function return all zero matrix.

assay <- "atac"
counts <- FeatureMatrix(
  fragments = atac@assays[["atac"]]@fragments[1],
  features = transcripts,
  cells = colnames(x = atac[["atac"]]),
  sep = c(":", "-")) ## return all zero matrix

After commenting out the cells argument, the returned matrix contains data but no sample index in the dimname.

counts <- FeatureMatrix(
  fragments = atac@assays[["atac"]]@fragments[1],
  features = transcripts,
  #cells = colnames(x = atac[["atac"]]),
  sep = c(":", "-")) ## return counts matrix with data. 

I followed the merge object vigette here. , and the barcodes in my fragment object look like this,

## cell barcode in the fragment objects
> atac@assays[["atac"]]@fragments[[1]]@cells %>% head
atac1_AAACGAAAGAGCAGCT-1 atac1_AAACGAAAGGAATGGA-1 atac1_AAACGAAAGGTTTCTG-1 
    "AAACGAAAGAGCAGCT-1"     "AAACGAAAGGAATGGA-1"     "AAACGAAAGGTTTCTG-1" 
atac1_AAACGAAAGTCGACCC-1 atac1_AAACGAAAGTCGATAA-1 atac1_AAACGAACAAACCCTA-1 
    "AAACGAAAGTCGACCC-1"     "AAACGAAAGTCGATAA-1"     "AAACGAACAAACCCTA-1" 


## cell barcodes in the atac assay
>colnames(x = atac[["atac"]]) %>% head
[1] "atac1_AAACGAAAGGAATGGA-1" "atac1_AAACGAAAGGTTTCTG-1" "atac1_AAACGAACAAACCCTA-1"
[4] "atac1_AAACGAACAGGTGTCC-1" "atac1_AAACTCGGTCACTCTC-1" "atac1_AAACTCGGTTTGTAGC-1"

I am using Signac 1.0.0 by the way.

[timoast]

I don't understand, did you run GeneActivity() or the code posted above? I'd suggest running GeneActivity() as shown in the vignette, making sure that you're using the correct gene annotations for the genome

[Ceeyouthere]

gene.activities <- GeneActivity(
  object = atac,
  assay = "atac")

This returns a all-zero matrix.
atac is a merged objects with 4 fragment objects, created by following the merging vignette.

I ran the GeneActivity() first and the returned matrix was all zero. Then I tried to troubleshoot and run the GeneActivity() function manually. I basically copied your code line by line and fill in the argument. Everything was fine until the FeatureMatrix() part, which was shown above. The annotation should not be the problem , since it works on unmerged atac object.

I ended up making gene activity matrix before merging and everything worked.
Thank you for your time and quick reply.

Here is the full code of manually running GeneActivity()

annotation <- Annotation(atac)
transcripts <- CollapseToLongestTranscript(ranges = annotation)
transcripts <- transcripts[transcripts$gene_biotype == "protein_coding"]
transcripts <- Extend(x = transcripts, upstream = 2000, 
                        downstream = 0)
assay <- "atac"
counts <- FeatureMatrix(
  fragments = atac@assays[["atac"]]@fragments,
  features = transcripts,
  cells = colnames(x = atac[["atac"]]),

Here is how the annotation was added. I used a self-made EnsDb package. It should be fine since I used on other analysis before.

DefaultAssay(atac) <- "atac"
devtools::load_all(path = "~/Desktop/500_Data/503_DataAnalysis/2020-06-16-annotation-inspection/EnsDb.Mmusculus.v100.v2/EnsDb.Mmusculus.v100/")
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v100)
# change to UCSC style since the data 
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "mm10"

# add the gene information to the object
Annotation(atac) <- annotations

[timoast]

Thanks for the clarification. One potential cause could be if the cell names are incorrect in the fragment object, it won't be matching any of the fragments in the file correctly. Are you able to successfully run any other functions that use the fragment files, for example make a coverage plot of a region?

[timoast]

I was able to reproduce this by merging datasets and running GeneActivity(), so I'll look into this now and hopefully have a fix available soon

[timoast]

Should be fixed now on the develop branch

[decarlin-d4c]

Hi, I'm getting this bug even after installing the development branch. I have Signac_1.1.1.9014 installed, but I've tried it with the main branch (1.1.1) as well. Specifically, going through the ATAC integration here:

https://satijalab.org/seurat/articles/atacseq_integration_vignette.html

I get an all zeros matrix from the GeneActivity() call. Here's my sessionInfo():

R version 4.0.3 (2020-10-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 18.04.4 LTS

Matrix products: default
BLAS/LAPACK: /encrypted_1/decarlin/anaconda3/envs/seurat_atac2/lib/libopenblasp-r0.3.12.so

locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] cowplot_1.1.1 ggplot2_3.3.3
[3] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.14.1
[5] AnnotationFilter_1.14.0 GenomicFeatures_1.42.3
[7] AnnotationDbi_1.52.0 Biobase_2.50.0
[9] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
[11] IRanges_2.24.1 S4Vectors_0.28.1
[13] BiocGenerics_0.36.1 Signac_1.1.1.9014
[15] SeuratObject_4.0.0 Seurat_4.0.1
[17] pbmcMultiome.SeuratData_0.1.0 SeuratData_0.2.1

loaded via a namespace (and not attached):
[1] utf8_1.2.1 reticulate_1.19
[3] tidyselect_1.1.0 RSQLite_2.2.5
[5] htmlwidgets_1.5.3 docopt_0.7.1
[7] grid_4.0.3 BiocParallel_1.24.1
[9] Rtsne_0.15 munsell_0.5.0
[11] codetools_0.2-18 ica_1.0-2
[13] future_1.21.0 miniUI_0.1.1.1
[15] withr_2.4.2 colorspace_2.0-0
[17] OrganismDbi_1.32.0 knitr_1.32
[19] rstudioapi_0.13 ROCR_1.0-11
[21] tensor_1.5 listenv_0.8.0
[23] MatrixGenerics_1.2.1 slam_0.1-48
[25] GenomeInfoDbData_1.2.4 polyclip_1.10-0
[27] farver_2.1.0 bit64_4.0.5
[29] parallelly_1.24.0 vctrs_0.3.7
[31] generics_0.1.0 xfun_0.22
[33] biovizBase_1.38.0 BiocFileCache_1.14.0
[35] lsa_0.73.2 ggseqlogo_0.1
[37] R6_2.5.0 reshape_0.8.8
[39] bitops_1.0-6 spatstat.utils_2.1-0
[41] cachem_1.0.4 DelayedArray_0.16.3
[43] assertthat_0.2.1 promises_1.2.0.1
[45] scales_1.1.1 nnet_7.3-15
[47] gtable_0.3.0 globals_0.14.0
[49] goftest_1.2-2 ggbio_1.38.0
[51] rlang_0.4.10 RcppRoll_0.3.0
[53] splines_4.0.3 rtracklayer_1.50.0
[55] lazyeval_0.2.2 dichromat_2.0-0
[57] spatstat.geom_2.1-0 checkmate_2.0.0
[59] BiocManager_1.30.12 reshape2_1.4.4
[61] abind_1.4-5 backports_1.2.1
[63] httpuv_1.5.5 Hmisc_4.5-0
[65] RBGL_1.66.0 tools_4.0.3
[67] ellipsis_0.3.1 spatstat.core_2.1-2
[69] RColorBrewer_1.1-2 ggridges_0.5.3
[71] Rcpp_1.0.6 plyr_1.8.6
[73] base64enc_0.1-3 progress_1.2.2
[75] zlibbioc_1.36.0 purrr_0.3.4
[77] RCurl_1.98-1.3 ps_1.6.0
[79] prettyunits_1.1.1 rpart_4.1-15
[81] openssl_1.4.3 deldir_0.2-10
[83] pbapply_1.4-3 zoo_1.8-9
[85] SummarizedExperiment_1.20.0 ggrepel_0.9.1
[87] cluster_2.1.1 magrittr_2.0.1
[89] data.table_1.14.0 scattermore_0.7
[91] lmtest_0.9-38 RANN_2.6.1
[93] SnowballC_0.7.0 ProtGenerics_1.22.0
[95] fitdistrplus_1.1-3 matrixStats_0.58.0
[97] hms_1.0.0 patchwork_1.1.1
[99] mime_0.10 xtable_1.8-4
[101] XML_3.99-0.6 jpeg_0.1-8.1
[103] sparsesvd_0.2 gridExtra_2.3
[105] compiler_4.0.3 biomaRt_2.46.3
[107] tibble_3.1.1 KernSmooth_2.23-18
[109] crayon_1.4.1 htmltools_0.5.1.1
[111] mgcv_1.8-34 later_1.2.0
[113] Formula_1.2-4 tidyr_1.1.3
[115] DBI_1.1.1 tweenr_1.0.2
[117] dbplyr_2.1.1 MASS_7.3-53.1
[119] rappdirs_0.3.3 Matrix_1.3-2
[121] cli_2.4.0 igraph_1.2.6
[123] pkgconfig_2.0.3 GenomicAlignments_1.26.0
[125] foreign_0.8-81 plotly_4.9.3
[127] spatstat.sparse_2.0-0 xml2_1.3.2
[129] XVector_0.30.0 stringr_1.4.0
[131] VariantAnnotation_1.36.0 digest_0.6.27
[133] sctransform_0.3.2 RcppAnnoy_0.0.18
[135] graph_1.68.0 spatstat.data_2.1-0
[137] Biostrings_2.58.0 fastmatch_1.1-0
[139] leiden_0.3.7 htmlTable_2.1.0
[141] uwot_0.1.10 curl_4.3
[143] shiny_1.6.0 Rsamtools_2.6.0
[145] lifecycle_1.0.0 nlme_3.1-152
[147] jsonlite_1.7.2 viridisLite_0.4.0
[149] askpass_1.1 BSgenome_1.58.0
[151] fansi_0.4.2 pillar_1.6.0
[153] GGally_2.1.1 lattice_0.20-41
[155] fastmap_1.1.0 httr_1.4.2
[157] survival_3.2-10 glue_1.4.2
[159] qlcMatrix_0.9.7 png_0.1-7
[161] bit_4.0.4 ggforce_0.3.3
[163] stringi_1.5.3 blob_1.2.1
[165] latticeExtra_0.6-29 memoise_2.0.0
[167] dplyr_1.0.5 irlba_2.3.3
[169] future.apply_1.7.0

Any help would be much appreciated, thanks for the package.