search

satijalab / seurat

R toolkit for single cell genomics
http://www.satijalab.org/seurat
Other
2.29k stars 913 forks source link

Empty gene activity during integration vignette #5162

Closed rylanci closed 3 years ago

rylanci commented 3 years ago

Hello I am following along the RNA-ATAC integration vignette and interestingly the gene.activity matrix output from GeneActivity() is empty. It contains correct row names derived from the variable features in pbmc.rna, however no overlapping regions appear to be counted. The vignette functions as described until this step. Are you familiar with anything that may cause this result in GeneActivity()?

library(SeuratData)
# install the dataset and load requirements
InstallData("pbmcMultiome")

library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
library(cowplot)

# load both modalities
pbmc.rna <- LoadData("pbmcMultiome", "pbmc.rna")
pbmc.atac <- LoadData("pbmcMultiome", "pbmc.atac")

# repeat QC steps performed in the WNN vignette
pbmc.rna <- subset(pbmc.rna, seurat_annotations != "filtered")
pbmc.atac <- subset(pbmc.atac, seurat_annotations != "filtered")

# Perform standard analysis of each modality independently RNA analysis
pbmc.rna <- NormalizeData(pbmc.rna)
pbmc.rna <- FindVariableFeatures(pbmc.rna)
pbmc.rna <- ScaleData(pbmc.rna)
pbmc.rna <- RunPCA(pbmc.rna)
pbmc.rna <- RunUMAP(pbmc.rna, dims = 1:30)

# ATAC analysis add gene annotation information
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "hg38"
Annotation(pbmc.atac) <- annotations

# We exclude the first dimension as this is typically correlated with sequencing depth
pbmc.atac <- RunTFIDF(pbmc.atac)
pbmc.atac <- FindTopFeatures(pbmc.atac, min.cutoff = "q0")
pbmc.atac <- RunSVD(pbmc.atac)
pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

# quantify gene activity
gene.activities <- GeneActivity(pbmc.atac, features = VariableFeatures(pbmc.rna))
# add gene activities as a new assay
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = gene.activities)
# normalize gene activities
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac, features = rownames(pbmc.atac))

sum(pbmc.atac@meta.data$nCount_ACTIVITY)
sum(pbmc.atac@meta.data$nFeature_ACTIVITY)
0
0

Row names are gene names as anticipated

head(rownames(pbmc.atac@assays$ACTIVITY))
'CSF2RA''XG''MID1''NHS''PHEX''PTCHD1'

Here is the SessionInfo():

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/rlancione/miniconda3/envs/signac/lib/libopenblasp-r0.3.17.so

locale:
[1] 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.5                
 [3] EnsDb.Hsapiens.v86_2.99.0     ensembldb_2.16.4             
 [5] AnnotationFilter_1.16.0       GenomicFeatures_1.44.0       
 [7] AnnotationDbi_1.54.1          Biobase_2.52.0               
 [9] GenomicRanges_1.44.0          GenomeInfoDb_1.28.1          
[11] IRanges_2.26.0                S4Vectors_0.30.0             
[13] BiocGenerics_0.38.0           Signac_1.3.0                 
[15] SeuratObject_4.0.2            Seurat_4.0.3                 
[17] pbmcMultiome.SeuratData_0.1.0 SeuratData_0.2.1             

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                  reticulate_1.20            
  [3] tidyselect_1.1.1            RSQLite_2.2.7              
  [5] htmlwidgets_1.5.3           grid_4.1.0                 
  [7] docopt_0.7.1                BiocParallel_1.26.1        
  [9] Rtsne_0.15                  munsell_0.5.0              
 [11] codetools_0.2-18            ica_1.0-2                  
 [13] pbdZMQ_0.3-5                future_1.21.0              
 [15] miniUI_0.1.1.1              withr_2.4.2                
 [17] colorspace_2.0-2            filelock_1.0.2             
 [19] knitr_1.33                  uuid_0.1-4                 
 [21] rstudioapi_0.13             ROCR_1.0-11                
 [23] tensor_1.5                  listenv_0.8.0              
 [25] labeling_0.4.2              MatrixGenerics_1.4.2       
 [27] slam_0.1-48                 repr_1.1.3                 
 [29] GenomeInfoDbData_1.2.6      polyclip_1.10-0            
 [31] bit64_4.0.5                 farver_2.1.0               
 [33] parallelly_1.27.0           vctrs_0.3.8                
 [35] generics_0.1.0              xfun_0.25                  
 [37] biovizBase_1.40.0           BiocFileCache_2.0.0        
 [39] lsa_0.73.2                  ggseqlogo_0.1              
 [41] R6_2.5.0                    bitops_1.0-7               
 [43] spatstat.utils_2.2-0        cachem_1.0.5               
 [45] DelayedArray_0.18.0         assertthat_0.2.1           
 [47] promises_1.2.0.1            BiocIO_1.2.0               
 [49] scales_1.1.1                nnet_7.3-16                
 [51] gtable_0.3.0                globals_0.14.0             
 [53] goftest_1.2-2               rlang_0.4.11               
 [55] RcppRoll_0.3.0              splines_4.1.0              
 [57] rtracklayer_1.52.0          lazyeval_0.2.2             
 [59] dichromat_2.0-0             checkmate_2.0.0            
 [61] spatstat.geom_2.2-2         yaml_2.2.1                 
 [63] reshape2_1.4.4              abind_1.4-5                
 [65] backports_1.2.1             httpuv_1.6.1               
 [67] Hmisc_4.5-0                 tools_4.1.0                
 [69] ellipsis_0.3.2              spatstat.core_2.3-0        
 [71] RColorBrewer_1.1-2          ggridges_0.5.3             
 [73] Rcpp_1.0.7                  plyr_1.8.6                 
 [75] base64enc_0.1-3             progress_1.2.2             
 [77] zlibbioc_1.38.0             purrr_0.3.4                
 [79] RCurl_1.98-1.3              prettyunits_1.1.1          
 [81] rpart_4.1-15                deldir_0.2-10              
 [83] pbapply_1.4-3               zoo_1.8-9                  
 [85] SummarizedExperiment_1.22.0 ggrepel_0.9.1              
 [87] cluster_2.1.2               magrittr_2.0.1             
 [89] RSpectra_0.16-0             data.table_1.14.0          
 [91] scattermore_0.7             lmtest_0.9-38              
 [93] RANN_2.6.1                  SnowballC_0.7.0            
 [95] ProtGenerics_1.24.0         fitdistrplus_1.1-5         
 [97] matrixStats_0.60.0          hms_1.1.0                  
 [99] patchwork_1.1.1             mime_0.11                  
[101] evaluate_0.14               xtable_1.8-4               
[103] XML_3.99-0.6                jpeg_0.1-9                 
[105] sparsesvd_0.2               gridExtra_2.3              
[107] compiler_4.1.0              biomaRt_2.48.2             
[109] tibble_3.1.3                KernSmooth_2.23-20         
[111] crayon_1.4.1                htmltools_0.5.1.1          
[113] mgcv_1.8-36                 later_1.2.0                
[115] Formula_1.2-4               tidyr_1.1.3                
[117] DBI_1.1.1                   tweenr_1.0.2               
[119] dbplyr_2.1.1                MASS_7.3-54                
[121] rappdirs_0.3.3              Matrix_1.3-4               
[123] cli_3.0.1                   igraph_1.2.6               
[125] pkgconfig_2.0.3             GenomicAlignments_1.28.0   
[127] foreign_0.8-81              IRdisplay_1.0              
[129] plotly_4.9.4.1              spatstat.sparse_2.0-0      
[131] xml2_1.3.2                  XVector_0.32.0             
[133] VariantAnnotation_1.38.0    stringr_1.4.0              
[135] digest_0.6.27               sctransform_0.3.2          
[137] RcppAnnoy_0.0.19            spatstat.data_2.1-0        
[139] Biostrings_2.60.2           leiden_0.3.9               
[141] fastmatch_1.1-3             htmlTable_2.2.1            
[143] uwot_0.1.10                 restfulr_0.0.13            
[145] curl_4.3.2                  shiny_1.6.0                
[147] Rsamtools_2.8.0             rjson_0.2.20               
[149] lifecycle_1.0.0             nlme_3.1-152               
[151] jsonlite_1.7.2              BSgenome_1.60.0            
[153] viridisLite_0.4.0           fansi_0.5.0                
[155] pillar_1.6.2                lattice_0.20-44            
[157] KEGGREST_1.32.0             fastmap_1.1.0              
[159] httr_1.4.2                  survival_3.2-11            
[161] glue_1.4.2                  qlcMatrix_0.9.7            
[163] png_0.1-7                   bit_4.0.4                  
[165] ggforce_0.3.3               stringi_1.7.3              
[167] blob_1.2.2                  latticeExtra_0.6-29        
[169] memoise_2.0.0               IRkernel_1.2               
[171] dplyr_1.0.7                 irlba_2.3.3                
[173] future.apply_1.8.1
saketkc commented 3 years ago

Hi @rylanci, I am unable to reproduce this issue. What do the first few entries of your gene activity matrix look like?

> gene.activities[1:5, 5:10]
5 x 6 sparse Matrix of class "dgCMatrix"
       AAACAGCCAGGATAAC-1 AAACAGCCAGTAGGTG-1 AAACAGCCAGTTTACG-1 AAACAGCCATCCAGGT-1 AAACATGCAAGGTCCT-1 AAACATGCACCGGCTA-1
CSF2RA                  .                  3                  .                  1                  .                  .
XG                      .                  .                  .                  .                  .                  .
MID1                    .                  .                  .                  .                  .                  .
NHS                     .                  .                  .                  .                  .                  .
PHEX                    .                  .                  1                  .                  .                  .
no-response[bot] commented 3 years ago

This issue has been automatically closed because there has been no response to our request for more information from the original author. With only the information that is currently in the issue, we don't have enough information to take action. Please reach out if you have or find the answers we need so that we can investigate further.