stuart-lab / signac

R toolkit for the analysis of single-cell chromatin data
https://stuartlab.org/signac/
Other
327 stars 88 forks source link

`CoveragePlot` Error in `colnames<-`(`*tmp*`, value = start(x = region):end(x = region)) : attempt to set 'colnames' on an object with less than two dimensions #1770

Closed danli349 closed 2 months ago

danli349 commented 2 months ago

Hello:

pwm_set <- getMatrixSet(x = JASPAR2020, opts = list(species = 9606, all_versions = FALSE))
motif.matrix <- CreateMotifMatrix(features = granges(pbmc), pwm = pwm_set, genome = BSgenome.Mmusculus.UCSC.mm10,  
                                  use.counts = FALSE)
motif.object <- CreateMotifObject(data = motif.matrix, pwm = pwm_set)
pbmc <- SetAssayData(pbmc, assay = 'ATAC', slot = 'motifs', new.data = motif.object)

# Note that this step can take 30-60 minutes 
pbmc <- RunChromVAR(
  object = pbmc,
  genome = BSgenome.Mmusculus.UCSC.mm10
)
R[write to console]: Computing GC bias per region

R[write to console]: Selecting background regions

R[write to console]: Computing deviations from background

R[write to console]: Constructing chromVAR assay

R[write to console]: Warning:
R[write to console]:  Layer counts isn't present in the assay object; returning NULL
Fragments(pbmc)
[[1]]
A Fragment object for 1028 cells
CoveragePlot(pbmc, region = 'Actb', assay = 'ATAC', 
             expression.assay = 'SCT', peaks = FALSE)
Error in `colnames<-`(`*tmp*`, value = start(x = region):end(x = region)) : 
  attempt to set 'colnames' on an object with less than two dimensions

Thanks

timoast commented 2 months ago

Please include the full code you're running and the output of sessionInfo()

danli349 commented 2 months ago

@timoast

library(Seurat)
library(ggplot2)
library(patchwork)
library(Signac)
library(EnsDb.Mmusculus.v79)
library(dplyr)
library(BSgenome.Mmusculus.UCSC.mm10)

inputdata.10x <- Read10X_h5("./data/filtered_feature_bc_matrix.h5")

# extract RNA and ATAC data
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# Create Seurat object
pbmc <- CreateSeuratObject(counts = rna_counts)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^mt-")

# Now add in the ATAC-seq data
# we'll only use peaks in standard chromosomes
grange.counts <- StringToGRanges(rownames(atac_counts), sep = c(":", "-"))
grange.use <- seqnames(grange.counts) %in% standardChromosomes(grange.counts)
atac_counts <- atac_counts[as.vector(grange.use), ]
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "mm10"

frag.file <- "./data/atac_fragments.tsv.gz"
chrom_assay <- CreateChromatinAssay(
   counts = atac_counts,
   sep = c(":", "-"),
   genome = 'mm10',
   fragments = frag.file,
   min.cells = 10,
   annotation = annotations
 )
pbmc[["ATAC"]] <- chrom_assay

# RNA analysis
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc, verbose = FALSE) %>% RunPCA() %>% RunUMAP(dims = 1:50, 
                                                                    reduction.name = 'umap.rna', reduction.key = 'rnaUMAP_')

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

pbmc <- FindMultiModalNeighbors(pbmc, reduction.list = list("pca", "lsi"), dims.list = list(1:50, 2:50))
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3, verbose = FALSE)

pbmc$celltype <- Idents(pbmc)

CoveragePlot(pbmc, region = 'Erg', features = 'sct_Erg', 
             assay = 'ATAC', expression.assay = 'SCT', peaks = FALSE)
Error in `colnames<-`(`*tmp*`, value = start(x = region):end(x = region)) : 
  attempt to set 'colnames' on an object with less than two dimensions
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 20.04.1 LTS

Matrix products: default
BLAS/LAPACK: /home/dan/miniconda3/lib/libmkl_rt.so.1;  LAPACK version 3.9.0

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   

time zone: America/New_York
tzcode source: system (glibc)

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

other attached packages:
 [1] scales_1.3.0                       viridis_0.6.5                     
 [3] viridisLite_0.4.2                  RColorBrewer_1.1-3                
 [5] BSgenome.Mmusculus.UCSC.mm10_1.4.3 BSgenome_1.72.0                   
 [7] rtracklayer_1.64.0                 BiocIO_1.14.0                     
 [9] Biostrings_2.72.1                  XVector_0.44.0                    
[11] dplyr_1.1.4                        EnsDb.Mmusculus.v79_2.99.0        
[13] ensembldb_2.28.0                   AnnotationFilter_1.28.0           
[15] GenomicFeatures_1.56.0             AnnotationDbi_1.66.0              
[17] Biobase_2.64.0                     GenomicRanges_1.56.1              
[19] GenomeInfoDb_1.40.1                IRanges_2.38.1                    
[21] S4Vectors_0.42.1                   BiocGenerics_0.50.0               
[23] Signac_1.13.0                      patchwork_1.2.0                   
[25] ggplot2_3.5.1                      Seurat_5.1.0                      
[27] SeuratObject_5.0.2                 sp_2.1-4                          

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22            splines_4.4.1              
  [3] later_1.3.2                 bitops_1.0-8               
  [5] tibble_3.2.1                polyclip_1.10-7            
  [7] rpart_4.1.23                XML_3.99-0.17              
  [9] fastDummies_1.7.3           lifecycle_1.0.4            
 [11] hdf5r_1.3.11                globals_0.16.3             
 [13] lattice_0.22-6              MASS_7.3-61                
 [15] backports_1.5.0             magrittr_2.0.3             
 [17] rmarkdown_2.27              Hmisc_5.1-3                
 [19] plotly_4.10.4               yaml_2.3.10                
 [21] httpuv_1.6.15               sctransform_0.4.1          
 [23] spam_2.10-0                 spatstat.sparse_3.1-0      
 [25] reticulate_1.38.0           cowplot_1.1.3              
 [27] pbapply_1.7-2               DBI_1.2.3                  
 [29] abind_1.4-5                 zlibbioc_1.50.0            
 [31] Rtsne_0.17                  purrr_1.0.2                
 [33] biovizBase_1.52.0           RCurl_1.98-1.16            
 [35] nnet_7.3-19                 VariantAnnotation_1.50.0   
 [37] GenomeInfoDbData_1.2.12     ggrepel_0.9.5              
 [39] irlba_2.3.5.1               listenv_0.9.1              
 [41] spatstat.utils_3.0-5        goftest_1.2-3              
 [43] RSpectra_0.16-2             spatstat.random_3.3-1      
 [45] fitdistrplus_1.2-1          parallelly_1.38.0          
 [47] leiden_0.4.3.1              codetools_0.2-20           
 [49] DelayedArray_0.30.1         RcppRoll_0.3.1             
 [51] tidyselect_1.2.1            farver_2.1.2               
 [53] UCSC.utils_1.0.0            base64enc_0.1-3            
 [55] matrixStats_1.3.0           spatstat.explore_3.3-1     
 [57] GenomicAlignments_1.40.0    jsonlite_1.8.8             
 [59] Formula_1.2-5               progressr_0.14.0           
 [61] ggridges_0.5.6              survival_3.7-0             
 [63] systemfonts_1.1.0           ragg_1.3.2                 
 [65] ica_1.0-3                   Rcpp_1.0.13                
 [67] glue_1.7.0                  gridExtra_2.3              
 [69] SparseArray_1.4.8           xfun_0.46                  
 [71] MatrixGenerics_1.16.0       withr_3.0.1                
 [73] fastmap_1.2.0               fansi_1.0.6                
 [75] digest_0.6.36               R6_2.5.1                   
 [77] mime_0.12                   textshaping_0.3.6          
 [79] colorspace_2.1-1            scattermore_1.2            
 [81] tensor_1.5                  dichromat_2.0-0.1          
 [83] spatstat.data_3.1-2         RSQLite_2.3.7              
 [85] utf8_1.2.4                  tidyr_1.3.1                
 [87] generics_0.1.3              data.table_1.15.4          
 [89] httr_1.4.7                  htmlwidgets_1.6.4          
 [91] S4Arrays_1.4.1              uwot_0.2.2                 
 [93] pkgconfig_2.0.3             gtable_0.3.5               
 [95] blob_1.2.4                  lmtest_0.9-40              
 [97] htmltools_0.5.8.1           dotCall64_1.1-1            
 [99] ProtGenerics_1.36.0         png_0.1-8                  
[101] spatstat.univar_3.0-0       rstudioapi_0.16.0          
[103] knitr_1.48                  reshape2_1.4.4             
[105] rjson_0.2.21                checkmate_2.3.2            
[107] nlme_3.1-166                curl_5.2.1                 
[109] zoo_1.8-12                  cachem_1.1.0               
[111] stringr_1.5.1               KernSmooth_2.23-24         
[113] vipor_0.4.7                 parallel_4.4.1             
[115] miniUI_0.1.1.1              foreign_0.8-87             
[117] ggrastr_1.0.2               restfulr_0.0.15            
[119] pillar_1.9.0                grid_4.4.1                 
[121] vctrs_0.6.5                 RANN_2.6.1                 
[123] promises_1.3.0              xtable_1.8-4               
[125] cluster_2.1.6               beeswarm_0.4.0             
[127] htmlTable_2.4.3             evaluate_0.24.0            
[129] cli_3.6.3                   compiler_4.4.1             
[131] Rsamtools_2.20.0            rlang_1.1.4                
[133] crayon_1.5.3                future.apply_1.11.2        
[135] labeling_0.4.3              ggbeeswarm_0.7.2           
[137] plyr_1.8.9                  stringi_1.8.4              
[139] deldir_2.0-4                BiocParallel_1.38.0        
[141] munsell_0.5.1               lazyeval_0.2.2             
[143] spatstat.geom_3.3-2         Matrix_1.7-0               
[145] RcppHNSW_0.6.0              bit64_4.0.5                
[147] future_1.34.0               KEGGREST_1.44.1            
[149] shiny_1.9.1                 SummarizedExperiment_1.34.0
[151] ROCR_1.0-11                 igraph_2.0.3               
[153] memoise_2.0.1               fastmatch_1.1-4            
[155] bit_4.0.5                  

Thanks a lot

timoast commented 2 months ago

Here you're changing the gene annotation chromosome names to UCSC-style:

seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "mm10"

From your previous issue, your data is mapped to the NCBI genome. You need to make sure the chromosome names in the gene annotations match those used in your data.