Error in infer_grn()

[TerezaClarence]

Dear Pndo developers,

I was wondering if you could help me to debug my code - I am just following your tutorial and applying it on my data below, however I ran into an error irrespective of which setup I use:

pbmc10
An object of class Seurat 
601236 features across 116471 samples within 4 assays 
Active assay: peaks (349688 features, 0 variable features)
 3 other assays present: RNA, ATAC, SCT
 7 dimensional reductions calculated: pca, harmony.pca, umap.rnah, lsi, harmony.lsi, umap.atach, wnn.umap

#---initialize GRN object
pbmc10 <- initiate_grn(pbmc10)

library(BSgenome.Hsapiens.UCSC.hg38)
data(motifs)

pbmc10 <- find_motifs(
    pbmc10,
    pfm = motifs,
    genome = BSgenome.Hsapiens.UCSC.hg38
)

pbmc10 <- infer_grn(
    pbmc10#,
    #peak_to_gene_method = 'GREAT',
    #method = 'brms'
)

Selecting candidate regulatory regions near genes

Error in .subset2(x, i, exact = exact): attempt to select less than one element in integerOneIndex
Traceback:

1. infer_grn(pbmc10)
2. infer_grn.SeuratPlus(pbmc10)
3. fit_grn_models(object = object, genes = genes, network_name = network_name, 
 .     peak_to_gene_method = peak_to_gene_method, upstream = upstream, 
 .     downstream = downstream, extend = extend, only_tss = only_tss, 
 .     parallel = parallel, tf_cor = tf_cor, peak_cor = peak_cor, 
 .     aggregate_rna_col = aggregate_rna_col, aggregate_peaks_col = aggregate_peaks_col, 
 .     method = method, alpha = alpha, family = family, interaction_term = interaction_term, 
 .     adjust_method = adjust_method, scale = scale, verbose = verbose, 
 .     ...)
4. fit_grn_models.SeuratPlus(object = object, genes = genes, network_name = network_name, 
 .     peak_to_gene_method = peak_to_gene_method, upstream = upstream, 
 .     downstream = downstream, extend = extend, only_tss = only_tss, 
 .     parallel = parallel, tf_cor = tf_cor, peak_cor = peak_cor, 
 .     aggregate_rna_col = aggregate_rna_col, aggregate_peaks_col = aggregate_peaks_col, 
 .     method = method, alpha = alpha, family = family, interaction_term = interaction_term, 
 .     adjust_method = adjust_method, scale = scale, verbose = verbose, 
 .     ...)
5. aggregate_matrix(t(peaks_near_gene), groups = colnames(peaks_near_gene), 
 .     fun = "sum")
6. fast_aggregate(x = x, groupings = groups, fun = fun)
7. data.frame(interaction(groupings2, sep = "_"))
8. interaction(groupings2, sep = "_")
9. as.factor(args[[i]])
10. is.factor(x)
11. args[[i]]
12. `[[.data.frame`(args, i)
13. (function(x, i, exact) if (is.matrix(i)) as.matrix(x)[[i]] else .subset2(x, 
  .     i, exact = exact))(x, ..., exact = exact)

I am not sure what could be the problem, I tried to DietSeurat, so it only contains RNA and peaks or RNA and ATAC assay but I still got the same problem.

Thank you so much for your help.

Best regards,
Tereza

>> sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 21.04

Matrix products: default


locale:
[1] C

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

other attached packages:
 [1] BSgenome.Hsapiens.UCSC.hg38_1.4.4 BSgenome_1.62.0                  
 [3] rtracklayer_1.54.0                Biostrings_2.62.0                
 [5] XVector_0.34.0                    Pando_1.0.3                      
 [7] regioneR_1.26.1                   GenomicRanges_1.46.1             
 [9] GenomeInfoDb_1.30.1               IRanges_2.28.0                   
[11] S4Vectors_0.32.4                  BiocGenerics_0.40.0              
[13] dplyr_1.0.10                      Signac_1.9.0                     
[15] SeuratObject_4.1.3                Seurat_4.3.0                     
[17] dreamlet_0.0.39                   variancePartition_1.25.13        
[19] BiocParallel_1.28.3               limma_3.50.3                     
[21] ggplot2_3.4.0                    

loaded via a namespace (and not attached):
  [1] pbdZMQ_0.3-8                scattermore_0.8            
  [3] mashr_0.2.69                R.methodsS3_1.8.2          
  [5] tidyr_1.2.1                 bit64_4.0.5                
  [7] R.utils_2.12.2              irlba_2.3.5.1              
  [9] DelayedArray_0.20.0         data.table_1.14.6          
 [11] TFBSTools_1.32.0            KEGGREST_1.34.0            
 [13] RCurl_1.98-1.9              doParallel_1.0.17          
 [15] generics_0.1.3              RhpcBLASctl_0.21-247.1     
 [17] cowplot_1.1.1               RSQLite_2.2.20             
 [19] RANN_2.6.1                  future_1.30.0              
 [21] ggpointdensity_0.1.0        tzdb_0.3.0                 
 [23] bit_4.0.5                   spatstat.data_3.0-0        
 [25] httpuv_1.6.7                SummarizedExperiment_1.24.0
 [27] assertthat_0.2.1            DirichletMultinomial_1.36.0
 [29] viridis_0.6.2               hms_1.1.2                  
 [31] babelgene_22.9              evaluate_0.19              
 [33] promises_1.2.0.1            fansi_1.0.3                
 [35] restfulr_0.0.15             progress_1.2.2             
 [37] caTools_1.18.2              Rgraphviz_2.38.0           
 [39] igraph_1.3.5                DBI_1.1.3                  
 [41] htmlwidgets_1.6.0           spatstat.geom_3.0-3        
 [43] purrr_1.0.0                 ellipsis_0.3.2             
 [45] backports_1.4.1             annotate_1.72.0            
 [47] aod_1.3.2                   deldir_1.0-6               
 [49] sparseMatrixStats_1.6.0     MatrixGenerics_1.6.0       
 [51] vctrs_0.5.1                 SingleCellExperiment_1.16.0
 [53] Biobase_2.54.0              ROCR_1.0-11                
 [55] abind_1.4-5                 cachem_1.0.6               
 [57] withr_2.5.0                 ggforce_0.4.1              
 [59] progressr_0.12.0            sctransform_0.3.5          
 [61] GenomicAlignments_1.30.0    prettyunits_1.1.1          
 [63] goftest_1.2-3               cluster_2.1.4              
 [65] seqLogo_1.60.0              IRdisplay_1.1              
 [67] lazyeval_0.2.2              crayon_1.5.2               
 [69] spatstat.explore_3.0-5      edgeR_3.36.0               
 [71] pkgconfig_2.0.3             tweenr_2.0.2               
 [73] nlme_3.1-161                pals_1.7                   
 [75] rlang_1.0.6                 globals_0.16.2             
 [77] lifecycle_1.0.3             miniUI_0.1.1.1             
 [79] dichromat_2.0-0.1           invgamma_1.1               
 [81] polyclip_1.10-4             matrixStats_0.63.0         
 [83] lmtest_0.9-40               graph_1.72.0               
 [85] Matrix_1.5-3                IRkernel_1.3.1             
 [87] ashr_2.2-54                 Rhdf5lib_1.16.0            
 [89] boot_1.3-28.1               zoo_1.8-11                 
 [91] base64enc_0.1-3             ggridges_0.5.4             
 [93] GlobalOptions_0.1.2         png_0.1-8                  
 [95] viridisLite_0.4.1           rjson_0.2.21               
 [97] bitops_1.0-7                R.oo_1.25.0                
 [99] KernSmooth_2.23-20          rhdf5filters_1.6.0         
[101] EnrichmentBrowser_2.25.3    blob_1.2.3                 
[103] DelayedMatrixStats_1.16.0   shape_1.4.6                
[105] mixsqp_0.3-48               stringr_1.5.0              
[107] SQUAREM_2021.1              parallelly_1.33.0          
[109] spatstat.random_3.0-1       readr_2.1.3                
[111] CNEr_1.30.0                 rmeta_3.0                  
[113] scales_1.2.1                memoise_2.0.1              
[115] GSEABase_1.56.0             magrittr_2.0.3             
[117] plyr_1.8.8                  ica_1.0-3                  
[119] gplots_3.1.3                zlibbioc_1.40.0            
[121] compiler_4.1.1              BiocIO_1.4.0               
[123] RColorBrewer_1.1-3          clue_0.3-63                
[125] KEGGgraph_1.54.0            lme4_1.1-31                
[127] fitdistrplus_1.1-8          Rsamtools_2.10.0           
[129] cli_3.5.0                   listenv_0.9.0              
[131] patchwork_1.1.2             pbapply_1.6-0              
[133] MASS_7.3-58.1               tidyselect_1.2.0           
[135] stringi_1.7.8               zenith_0.99.1              
[137] yaml_2.3.6                  locfit_1.5-9.7             
[139] ggrepel_0.9.2               grid_4.1.1                 
[141] fastmatch_1.1-3             tools_4.1.1                
[143] future.apply_1.10.0         parallel_4.1.1             
[145] circlize_0.4.15             uuid_1.1-0                 
[147] TFMPvalue_0.0.9             foreach_1.5.2              
[149] gridExtra_2.3               farver_2.1.1               
[151] Rtsne_0.16                  ggraph_2.1.0               
[153] RcppZiggurat_0.1.6          digest_0.6.31              
[155] pracma_2.4.2                shiny_1.7.4                
[157] motifmatchr_1.16.0          Rcpp_1.0.9                 
[159] broom_1.0.2                 later_1.3.0                
[161] RcppAnnoy_0.0.20            httr_1.4.4                 
[163] AnnotationDbi_1.56.2        ComplexHeatmap_2.10.0      
[165] Rdpack_2.4                  colorspace_2.0-3           
[167] XML_3.99-0.13               tensor_1.5                 
[169] reticulate_1.26             truncnorm_1.0-8            
[171] splines_4.1.1               uwot_0.1.14                
[173] RcppRoll_0.3.0              spatstat.utils_3.0-1       
[175] graphlayouts_0.8.4          sp_1.5-1                   
[177] mapproj_1.2.11              plotly_4.10.1              
[179] xtable_1.8-4                poweRlaw_0.70.6            
[181] jsonlite_1.8.4              nloptr_2.0.3               
[183] tidygraph_1.2.2             Rfast_2.0.6                
[185] R6_2.5.1                    pillar_1.8.1               
[187] htmltools_0.5.4             mime_0.12                  
[189] glue_1.6.2                  fastmap_1.1.0              
[191] minqa_1.2.5                 codetools_0.2-18           
[193] maps_3.4.1                  mvtnorm_1.1-3              
[195] utf8_1.2.2                  lattice_0.20-45            
[197] spatstat.sparse_3.0-0       tibble_3.1.8               
[199] pbkrtest_0.5.1              leiden_0.4.3               
[201] gtools_3.9.4                GO.db_3.14.0               
[203] survival_3.4-0              repr_1.1.4                 
[205] munsell_0.5.0               GetoptLong_1.0.5           
[207] rhdf5_2.38.1                GenomeInfoDbData_1.2.7     
[209] iterators_1.0.14            HDF5Array_1.22.1           
[211] reshape2_1.4.4              gtable_0.3.1               
[213] msigdbr_7.5.1               rbibutils_2.2.11  

[joschif]

Hi @TerezaClarence , I'm not fully sure what causes the issue here, are you using a publically PBMC available dataset? If you point me to where you have it from then I can try to reproduce the error

[TerezaClarence]

Dear @joschif I am using my multiome dataset, data were loaded from cellranger 10x multiome and preprocessed with seurat and signac. I have no idea why it's giving me the error - I tried to repeatedly just diet the data so it contains RNA and ATC-assay, but it still gives me the same error. Changing the DefaultAssay also didn't help.

Best,
Tereza

[Sandman-1]

Hello @TerezaClarence! Just to make sure, the RNA assay does have variable features, correct?

[ahoffrichter]

Hello,
I used the scMultiome analysis vignette on our unpublished data in which in Step3 a gene regulatory network reconstruction is performed using pando.
When I run infer_grn I get the same error message as mentioned above:

  attempt to select less than one element in integerOneIndex

13: (function(x, i, exact) if (is.matrix(i)) as.matrix(x)[[i]] else .subset2(x, 
        i, exact = exact))(x, ..., exact = exact)
12: `[[.data.frame`(args, i)
11: args[[i]]
10: is.factor(x)
9: as.factor(args[[i]])
8: interaction(groupings2, sep = "_")
7: data.frame(interaction(groupings2, sep = "_"))
6: fast_aggregate(x = x, groupings = groups, fun = fun)
5: aggregate_matrix(t(peaks_near_gene), groups = colnames(peaks_near_gene), 
       fun = "sum")
4: fit_grn_models.SeuratPlus(object = object, genes = genes, network_name = network_name, 
       peak_to_gene_method = peak_to_gene_method, upstream = upstream, 
       downstream = downstream, extend = extend, only_tss = only_tss, 
       parallel = parallel, tf_cor = tf_cor, peak_cor = peak_cor, 
       aggregate_rna_col = aggregate_rna_col, aggregate_peaks_col = aggregate_peaks_col, 
       method = method, alpha = alpha, family = family, interaction_term = interaction_term, 
       adjust_method = adjust_method, scale = scale, verbose = verbose, 
       ...)
3: fit_grn_models(object = object, genes = genes, network_name = network_name, 
       peak_to_gene_method = peak_to_gene_method, upstream = upstream, 
       downstream = downstream, extend = extend, only_tss = only_tss, 
       parallel = parallel, tf_cor = tf_cor, peak_cor = peak_cor, 
       aggregate_rna_col = aggregate_rna_col, aggregate_peaks_col = aggregate_peaks_col, 
       method = method, alpha = alpha, family = family, interaction_term = interaction_term, 
       adjust_method = adjust_method, scale = scale, verbose = verbose, 
       ...)
2: infer_grn.SeuratPlus(seurat, parallel = F, tf_cor = 0.05, method = "glm", 
       family = "gaussian", scale = F, verbose = T)
1: infer_grn(seurat, parallel = F, tf_cor = 0.05, method = "glm", 
       family = "gaussian", scale = F, verbose = T)

> sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Ventura 13.3.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] doParallel_1.0.17                 iterators_1.0.14                  foreach_1.5.2                     Pando_1.0.3                       JASPAR2020_0.99.10               
 [6] TFBSTools_1.36.0                  presto_1.0.0                      data.table_1.14.8                 harmony_0.1.1                     Rcpp_1.0.10                      
[11] ensembldb_2.22.0                  AnnotationFilter_1.22.0           GenomicFeatures_1.50.4            AnnotationDbi_1.60.2              Biobase_2.58.0                   
[16] simspec_0.0.0.9000                BSgenome.Hsapiens.UCSC.hg38_1.4.5 BSgenome_1.66.3                   rtracklayer_1.58.0                Biostrings_2.66.0                
[21] XVector_0.38.0                    GenomicRanges_1.50.2              GenomeInfoDb_1.34.9               IRanges_2.32.0                    S4Vectors_0.36.2                 
[26] AnnotationHub_3.6.0               BiocFileCache_2.6.1               dbplyr_2.3.2                      BiocGenerics_0.44.0               Signac_1.9.0                     
[31] SeuratObject_4.1.3                Seurat_4.3.0                      patchwork_1.1.2                   lubridate_1.9.2                   forcats_1.0.0                    
[36] stringr_1.5.0                     dplyr_1.1.2                       purrr_1.0.1                       readr_2.1.4                       tidyr_1.3.0                      
[41] tibble_3.2.1                      ggplot2_3.4.2                     tidyverse_2.0.0                  

loaded via a namespace (and not attached):
  [1] Hmisc_5.1-0                   ica_1.0-3                     RcppRoll_0.3.0                ps_1.7.5                      Rsamtools_2.14.0             
  [6] lmtest_0.9-40                 rprojroot_2.0.3               crayon_1.5.2                  MASS_7.3-60                   nlme_3.1-162                 
 [11] backports_1.4.1               qlcMatrix_0.9.7               rlang_1.1.1                   ROCR_1.0-11                   irlba_2.3.5.1                
 [16] callr_3.7.3                   limma_3.54.2                  filelock_1.0.2                BiocParallel_1.32.6           rjson_0.2.21                 
 [21] CNEr_1.34.0                   bit64_4.0.5                   glue_1.6.2                    poweRlaw_0.70.6               sctransform_0.3.5            
 [26] processx_3.8.1                vipor_0.4.5                   chromVAR_1.20.2               spatstat.sparse_3.0-1         spatstat.geom_3.2-1          
 [31] tidyselect_1.2.0              SummarizedExperiment_1.28.0   usethis_2.1.6                 fitdistrplus_1.1-11           XML_3.99-0.14                
 [36] zoo_1.8-12                    GenomicAlignments_1.34.1      xtable_1.8-4                  magrittr_2.0.3                evaluate_0.21                
 [41] cli_3.6.1                     zlibbioc_1.44.0               rstudioapi_0.14               miniUI_0.1.1.1                sp_1.6-0                     
 [46] rpart_4.1.19                  fastmatch_1.1-3               maps_3.4.1                    shiny_1.7.4                   xfun_0.39                    
 [51] pkgbuild_1.4.0                cluster_2.1.4                 caTools_1.18.2                tidygraph_1.2.3               KEGGREST_1.38.0              
 [56] interactiveDisplayBase_1.36.0 ggrepel_0.9.3                 biovizBase_1.46.0             listenv_0.9.0                 TFMPvalue_0.0.9              
 [61] png_0.1-8                     future_1.32.0                 withr_2.5.0                   bitops_1.0-7                  slam_0.1-50                  
 [66] ggforce_0.4.1                 plyr_1.8.8                    sparsesvd_0.2-2               pracma_2.4.2                  pillar_1.9.0                 
 [71] cachem_1.0.8                  fs_1.6.2                      hdf5r_1.3.8                   vctrs_0.6.2                   ellipsis_0.3.2               
 [76] generics_0.1.3                devtools_2.4.5                tools_4.2.2                   foreign_0.8-84                beeswarm_0.4.0               
 [81] munsell_0.5.0                 tweenr_2.0.2                  DelayedArray_0.24.0           fastmap_1.1.1                 compiler_4.2.2               
 [86] pkgload_1.3.2                 abind_1.4-5                   httpuv_1.6.11                 sessioninfo_1.2.2             plotly_4.10.1                
 [91] GenomeInfoDbData_1.2.9        gridExtra_2.3                 lattice_0.21-8                ggpointdensity_0.1.0          deldir_1.0-9                 
 [96] utf8_1.2.3                    later_1.3.1                   jsonlite_1.8.4                scales_1.2.1                  docopt_0.7.1                 
[101] pbapply_1.7-0                 sparseMatrixStats_1.10.0      lazyeval_0.2.2                nabor_0.5.0                   promises_1.2.0.1             
[106] R.utils_2.12.2                goftest_1.2-3                 spatstat.utils_3.0-3          reticulate_1.28               checkmate_2.2.0              
[111] rmarkdown_2.21                cowplot_1.1.1                 Rtsne_0.16                    dichromat_2.0-0.1             uwot_0.1.14                  
[116] igraph_1.4.3                  survival_3.5-5                yaml_2.3.7                    htmltools_0.5.5               memoise_2.0.1                
[121] VariantAnnotation_1.44.1      profvis_0.3.8                 BiocIO_1.8.0                  graphlayouts_1.0.0            viridisLite_0.4.2            
[126] digest_0.6.31                 mime_0.12                     rappdirs_0.3.3                pals_1.7                      RSQLite_2.3.1                
[131] future.apply_1.11.0           mapproj_1.2.11                remotes_2.4.2                 urlchecker_1.0.1              blob_1.2.4                   
[136] R.oo_1.25.0                   splines_4.2.2                 Formula_1.2-5                 labeling_0.4.2                ProtGenerics_1.30.0          
[141] RCurl_1.98-1.12               hms_1.1.3                     colorspace_2.1-0              base64enc_0.1-3               BiocManager_1.30.20          
[146] ggbeeswarm_0.7.2              ggrastr_1.0.1                 nnet_7.3-19                   RANN_2.6.1                    ggseqlogo_0.1                
[151] fansi_1.0.4                   tzdb_0.4.0                    parallelly_1.35.0             R6_2.5.1                      grid_4.2.2                   
[156] ggridges_0.5.4                lifecycle_1.0.3               curl_5.0.0                    leiden_0.4.3                  motifmatchr_1.20.0           
[161] Matrix_1.5-4.1                desc_1.4.2                    RcppAnnoy_0.0.20              RColorBrewer_1.1-3            spatstat.explore_3.2-1       
[166] htmlwidgets_1.6.2             polyclip_1.10-4               biomaRt_2.54.1                timechange_0.2.0              seqLogo_1.64.0               
[171] globals_0.16.2                htmlTable_2.4.1               spatstat.random_3.1-5         progressr_0.13.0              codetools_0.2-19             
[176] matrixStats_0.63.0            GO.db_3.16.0                  gtools_3.9.4                  prettyunits_1.1.1             R.methodsS3_1.8.2            
[181] gtable_0.3.3                  DBI_1.1.3                     tensor_1.5                    httr_1.4.6                    KernSmooth_2.23-21           
[186] stringi_1.7.12                progress_1.2.2                reshape2_1.4.4                farver_2.1.1                  annotate_1.76.0              
[191] viridis_0.6.3                 DT_0.28                       xml2_1.3.4                    restfulr_0.0.15               scattermore_1.1              
[196] BiocVersion_3.16.0            bit_4.0.5                     MatrixGenerics_1.10.0         spatstat.data_3.0-1           ggraph_2.1.0                 
[201] pkgconfig_2.0.3               DirichletMultinomial_1.40.0   knitr_1.43    

[AmelZulji]

Not part of dev team but might be useful:

In my understanding infer_grn() will by default fit models for genes found as VariableFeatures. If VariableFeatures are not present in the object, it will throw the error that @TerezaClarence reported.

One can solve the problem by specifying genes using genes argument in the infer_grn().