search

yuchaojiang / TRIPOD

Nonparametric Interrogation of Transcriptional Regulation in Single-Cell RNA and Chromatin Accessibility Multiomic Data
GNU General Public License v3.0
7 stars 1 forks source link

getTrios returns empty #3

Closed DLGisch closed 2 years ago

DLGisch commented 2 years ago

Hello,

I did your last tip I fix the first problem. Now when the function getTrios returns empty and I don't receive any error. I am running the vignette data and the same two genes c("CCR7", "GNLY")

Whole cod:

library(Seurat)
library(Signac)
library(GenomicRanges)
library(dplyr)
library(ggplot2)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomeInfoDb)
library(patchwork)
library(BiocParallel)
library(TRIPOD)
library(stats)

pbmc <- readRDS(paste0(data,'pbmc.chromvar.annotated.rds'))
tripod.pbmc <- getObjectsForModelFit(object = pbmc, chr = paste0("chr", 1:22),
                                     convert = TRUE)

pbmc <- filterSeuratObject(object = pbmc, tripod.object = tripod.pbmc)
pbmc <- processSeuratObject(object = pbmc, dim.rna = 1:50, dim.atac = 2:50,
                            verbose = FALSE)
set.seed(123)

num.clusters <- optimizeResolution(
  object = pbmc,
  graph.name = "wsnn",
  assay.name = "WNN",
  resolutions = seq(10, 35, 5),
  min.num = 20
)
num.clusters

res <- 15
pbmc <- getClusters(object = pbmc, graph.name = "wsnn", algorithm = 3,
                    resolution = res, verbose = FALSE)

metacell.pbmc <- getMetacellMatrices(object = pbmc,
                                     cluster.name = "seurat_clusters"
)
color.pbmc <- getColors(object = pbmc)

pbmc.transcripts.gr=tripod.pbmc$transcripts.gr
pbmc.peaks.gr=tripod.pbmc$peaks.gr
pbmc.peakxmotif=tripod.pbmc$peakxmotif

pbmc.metacell.rna=metacell.pbmc$metacell.rna
pbmc.metacell.peak=metacell.pbmc$metacell.peak

pbmc.metacell.celltype=color.pbmc$metacell.celltype
pbmc.metacell.celltype.col=color.pbmc$metacell.celltype.col

genes <- c("CCR7", "GNLY")
xymats.list <- list()
xymats.tripod.Xt.list <- list()
xymats.tripod.Yj.list <- list()
for (g in genes){
xymats.list[[g]] <-   getXYMatrices(gene.name = g,
  ext.upstream =ext.upstream,
  ext.downstream = ext.downstream,
  transcripts.gr = pbmc.transcripts.gr,
  peaks.gr = pbmc.peaks.gr,
  metacell.rna = pbmc.metacell.rna,
  metacell.peak = pbmc.metacell.peak,
  peakxmotif = pbmc.peakxmotif,
  motifxTF = pbmc.motifxTF,
  metacell.celltype = pbmc.metacell.celltype,
  metacell.celltype.col =pbmc.metacell.celltype.col)

# run TRIPOD matching Xt
xymats.tripod.Xt.list[[g]] <-  fitModel(xymats.list[[g]],
  model.name = "TRIPOD",
  match.by = "Xt"
)

# run TRIPOD matching Yj
xymats.tripod.Yj.list[[g]] <- fitModel(
  xymats.list[[g]],
  model.name = "TRIPOD",
  match.by = "Yj"
)
}

xymats.tX1.pos.df <- getTrios(
  xymats.list = xymats.tripod.Xt.list,
  fdr.thresh = fdr.thresh,
  sign = sign,
  model.name = "TRIPOD",
  level = 1
)

Can you help me again? Thank you, Debora

yuchaojiang commented 2 years ago

I personally ran this code just a few days ago and had no issues. At which step did you encounter the problem? which object is empty?

DLGisch commented 2 years ago

xymats.tX1.pos.df

yuchaojiang commented 2 years ago

Can you send the entire script, output, and the potential error and part that doesn't make sense ? There is no way that I could tell what went wrong just based off of this. What is your fdr.thres? You did not specify it in the script above.

DLGisch commented 2 years ago

Hi,

I installed fresh in another machine, and got the exact same error.

Here is the complete output:

R version 4.2.1 (2022-06-23) -- "Funny-Looking Kid"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Seurat)
Attaching SeuratObject
Attaching sp
Warning messages:
1: R graphics engine version 15 is not supported by this version of RStudio. The Plots tab will be disabled until a newer version of RStudio is installed. 
2: In fun(libname, pkgname) :
  rgeos: versions of GEOS runtime 3.10.1-CAPI-1.16.0
and GEOS at installation 3.5.1-CAPI-1.9.1differ
> library(Signac)
> library(GenomicRanges)
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: ?BiocGenerics?

The following objects are masked from ?package:stats?:

    IQR, mad, sd, var, xtabs

The following objects are masked from ?package:base?:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname,
    do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect,
    is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int,
    pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ?S4Vectors?

The following objects are masked from ?package:base?:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ?IRanges?

The following object is masked from ?package:sp?:

    %over%

Loading required package: GenomeInfoDb
> library(dplyr)

Attaching package: ?dplyr?

The following objects are masked from ?package:GenomicRanges?:

    intersect, setdiff, union

The following object is masked from ?package:GenomeInfoDb?:

    intersect

The following objects are masked from ?package:IRanges?:

    collapse, desc, intersect, setdiff, slice, union

The following objects are masked from ?package:S4Vectors?:

    first, intersect, rename, setdiff, setequal, union

The following objects are masked from ?package:BiocGenerics?:

    combine, intersect, setdiff, union

The following objects are masked from ?package:stats?:

    filter, lag

The following objects are masked from ?package:base?:

    intersect, setdiff, setequal, union

> library(ggplot2)
Learn more about the underlying theory at https://ggplot2-book.org/
> library(EnsDb.Hsapiens.v86)
Loading required package: ensembldb
Loading required package: GenomicFeatures
Loading required package: AnnotationDbi
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite
    Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.

Attaching package: ?AnnotationDbi?

The following object is masked from ?package:dplyr?:

    select

Loading required package: AnnotationFilter

Attaching package: 'ensembldb'

The following object is masked from 'package:dplyr':

    filter

The following object is masked from 'package:stats':

    filter

> library(BSgenome.Hsapiens.UCSC.hg38)
Loading required package: BSgenome
Loading required package: Biostrings
Loading required package: XVector

Attaching package: 'Biostrings'

The following object is masked from 'package:base':

    strsplit

Loading required package: rtracklayer
> library(GenomeInfoDb)
> library(patchwork)
> library(BiocParallel)
> library(TRIPOD)
Warning messages:
1: replacing previous import 'GenomicRanges::union' by 'dplyr::union' when loading 'TRIPOD' 
2: replacing previous import 'GenomicRanges::intersect' by 'dplyr::intersect' when loading 'TRIPOD' 
3: replacing previous import 'GenomicRanges::setdiff' by 'dplyr::setdiff' when loading 'TRIPOD' 
> library(stats)
> date <-  format(Sys.time(), "%a_%b_%d")
> data <- './'
> pbmc <- readRDS(paste0(data,'pbmc.chromvar.annotated.rds'))
> tripod.pbmc <- getObjectsForModelFit(object = pbmc, chr = paste0("chr", 1:22),
+                                      convert = TRUE)
> transcripts.gr <- tripod.pbmc$transcripts.gr
> peaks.gr <- tripod.pbmc$peaks.gr
> motifxTF <- tripod.pbmc$motifxTF
> peakxmotif <- tripod.pbmc$peakxmotif
> pbmc <- filterSeuratObject(object = pbmc, tripod.object = tripod.pbmc)
> pbmc <- processSeuratObject(object = pbmc, dim.rna = 1:50, dim.atac = 2:50,
+                             verbose = FALSE)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
Computing nearest neighbor graph
Computing SNN
Computing nearest neighbor graph
Computing SNN
> set.seed(123)
> num.clusters <- optimizeResolution(
+   object = pbmc,
+   graph.name = "wsnn",
+   assay.name = "WNN",
+   resolutions = seq(10, 35, 5),
+   min.num = 20
+ )
> num.clusters
  resolution num_clusters num_below
1         10           61         2
2         15           84         1
3         20          111         5
4         25          138         8
5         30          165        13
6         35          201        26
> res <- 15
> pbmc <- getClusters(object = pbmc, graph.name = "wsnn", algorithm = 3,
+                     resolution = res, verbose = FALSE)
> metacell.pbmc <- getMetacellMatrices(object = pbmc,
+                                      cluster.name = "seurat_clusters"
+ )
> metacell.rna <- metacell.pbmc$rna
> metacell.peak <- metacell.pbmc$peak
> color.pbmc <- getColors(object = pbmc)
> metacell.celltype <- color.pbmc$metacell$celltype
> metacell.celltype.col <- color.pbmc$metacell$color
> # obtain highly variable genes based on SCT
> DefaultAssay(pbmc) <- "SCT"
> hvg.pbmc <- VariableFeatures(pbmc)
>  pbmc.transcripts.gr=transcripts.gr
> # saveRDS(pbmc.transcripts.gr, file =paste0(data,'pbmc.transcripts.gr.RDS'))
>  pbmc.peaks.gr=peaks.gr
> # saveRDS(pbmc.peaks.gr, file=paste0(data,'pbmc.peaks.gr.RDS'))
>  pbmc.motifxTF=motifxTF
> # saveRDS(pbmc.motifxTF, file=paste0(data,'pbmc.motifxTF.RDS'))
>  pbmc.peakxmotif=peakxmotif
> # saveRDS(pbmc.peakxmotif, file=paste0(data,'pbmc.peakxmotif.RDS'))
>  pbmc.metacell.rna=metacell.rna
> # saveRDS(pbmc.metacell.rna, file=paste0(data,'pbmc.metacell.rna.RDS'))
>  pbmc.metacell.peak=metacell.peak
> # saveRDS(pbmc.metacell.peak, file=paste0(data,'pbmc.metacell.peak.RDS'))
>  pbmc.metacell.celltype=metacell.celltype
> # saveRDS(pbmc.metacell.celltype, file=paste0(data,'pbmc.metacell.celltype.RDS'))
>  pbmc.metacell.celltype.col=metacell.celltype.col
> genes <- c("CCR7", "GNLY")
> ext.upstream <- ext.downstream <- 1e5
> xymats.list <- bplapply(
+   genes,
+   getXYMatrices,
+   ext.upstream = ext.upstream,
+   transcripts.gr = pbmc.transcripts.gr,
+   peaks.gr = pbmc.peaks.gr,
+   metacell.rna = pbmc.metacell.rna,
+   metacell.peak = pbmc.metacell.peak,
+   peakxmotif = pbmc.peakxmotif,
+   motifxTF = pbmc.motifxTF,
+   metacell.celltype = pbmc.metacell.celltype,
+   metacell.celltype.col = pbmc.metacell.celltype.col
+ )
> names(xymats.list) <- genes
> # run TRIPOD matching Xt
> xymats.tripod.Xt.list <- bplapply(
+   xymats.list,
+   fitModel,
+   model.name = "TRIPOD",
+   match.by = "Xt"
+ )
> names(xymats.tripod.Xt.list) <- genes
> # run TRIPOD matching Yj
> xymats.tripod.Yj.list <- bplapply(
+   xymats.list,
+   fitModel,
+   model.name = "TRIPOD",
+   match.by = "Yj"
+ )
> names(xymats.tripod.Yj.list) <- genes
> #xymats.tripod.Yj.list <- readRDS(paste0('pbmc/xymats.tripod.Yj.list.Wed_Oct_12.RDS'))
> # set FDR < 0.01
> fdr.thresh <- 0.01
> # focus on positive sign
> sign <- "positive"
> # TRIPOD level 1 matching Xt
> xymats.tX1.pos.df <- getTrios(
+   xymats.list = xymats.tripod.Xt.list,
+   fdr.thresh = fdr.thresh,
+   sign = sign,
+   model.name = "TRIPOD",
+   level = 1
+ )
> xymats.tX1.pos.df
[1] gene     peak_num TF_num   peak     TF       coef     pval     adj     
<0 rows> (or 0-length row.names)

and the sessionInfo is:

sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux

Matrix products: default
BLAS:   /geode2/soft/hps/rhel7/r/4.2.1/lib64/R/lib/libRblas.so
LAPACK: /geode2/soft/hps/rhel7/r/4.2.1/lib64/R/lib/libRlapack.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] TRIPOD_1.0.0                      BiocParallel_1.30.4              
 [3] patchwork_1.1.2                   BSgenome.Hsapiens.UCSC.hg38_1.4.4
 [5] BSgenome_1.64.0                   rtracklayer_1.56.1               
 [7] Biostrings_2.64.1                 XVector_0.36.0                   
 [9] EnsDb.Hsapiens.v86_2.99.0         ensembldb_2.20.2                 
[11] AnnotationFilter_1.20.0           GenomicFeatures_1.48.4           
[13] AnnotationDbi_1.58.0              Biobase_2.56.0                   
[15] ggplot2_3.3.6                     dplyr_1.0.10                     
[17] GenomicRanges_1.48.0              GenomeInfoDb_1.32.4              
[19] IRanges_2.30.1                    S4Vectors_0.34.0                 
[21] BiocGenerics_0.42.0               Signac_1.8.0                     
[23] sp_1.5-0                          SeuratObject_4.1.2               
[25] Seurat_4.2.0                     

loaded via a namespace (and not attached):
  [1] utf8_1.2.2                  reticulate_1.26             tidyselect_1.1.2           
  [4] RSQLite_2.2.18              htmlwidgets_1.5.4           grid_4.2.1                 
  [7] Rtsne_0.16                  munsell_0.5.0               codetools_0.2-18           
 [10] ica_1.0-3                   interp_1.1-3                future_1.28.0              
 [13] miniUI_0.1.1.1              withr_2.5.0                 spatstat.random_2.2-0      
 [16] colorspace_2.0-3            progressr_0.11.0            filelock_1.0.2             
 [19] knitr_1.39                  rstudioapi_0.13             ROCR_1.0-11                
 [22] tensor_1.5                  listenv_0.8.0               labeling_0.4.2             
 [25] MatrixGenerics_1.8.1        olsrr_0.5.3                 GenomeInfoDbData_1.2.8     
 [28] polyclip_1.10-0             farver_2.1.1                bit64_4.0.5                
 [31] parallelly_1.32.1           vctrs_0.4.1                 generics_0.1.3             
 [34] xfun_0.31                   BiocFileCache_2.4.0         R6_2.5.1                   
 [37] bitops_1.0-7                spatstat.utils_2.3-1        cachem_1.0.6               
 [40] DelayedArray_0.22.0         assertthat_0.2.1            promises_1.2.0.1           
 [43] BiocIO_1.6.0                scales_1.2.1                nnet_7.3-18                
 [46] rgeos_0.5-9                 gtable_0.3.1                globals_0.16.1             
 [49] goftest_1.2-3               rlang_1.0.4                 RcppRoll_0.3.0             
 [52] splines_4.2.1               lazyeval_0.2.2              checkmate_2.1.0            
 [55] spatstat.geom_2.4-0         yaml_2.3.5                  reshape2_1.4.4             
 [58] abind_1.4-5                 backports_1.4.1             httpuv_1.6.5               
 [61] Hmisc_4.7-1                 tools_4.2.1                 ellipsis_0.3.2             
 [64] spatstat.core_2.4-4         RColorBrewer_1.1-3          ggridges_0.5.4             
 [67] Rcpp_1.0.9                  plyr_1.8.7                  base64enc_0.1-3            
 [70] progress_1.2.2              zlibbioc_1.42.0             purrr_0.3.4                
 [73] RCurl_1.98-1.9              prettyunits_1.1.1           rpart_4.1.16               
 [76] deldir_1.0-6                pbapply_1.5-0               viridis_0.6.2              
 [79] cowplot_1.1.1               zoo_1.8-11                  SummarizedExperiment_1.26.1
 [82] ggrepel_0.9.1               cluster_2.1.4               magrittr_2.0.3             
 [85] data.table_1.14.2           scattermore_0.8             lmtest_0.9-40              
 [88] RANN_2.6.1                  ProtGenerics_1.28.0         fitdistrplus_1.1-8         
 [91] matrixStats_0.62.0          hms_1.1.2                   mime_0.12                  
 [94] xtable_1.8-4                XML_3.99-0.11               jpeg_0.1-9                 
 [97] shape_1.4.6                 gridExtra_2.3               compiler_4.2.1             
[100] biomaRt_2.52.0              tibble_3.1.8                KernSmooth_2.23-20         
[103] crayon_1.5.1                htmltools_0.5.3             mgcv_1.8-40                
[106] later_1.3.0                 Formula_1.2-4               tidyr_1.2.1                
[109] nbpMatching_1.5.1           DBI_1.1.3                   dbplyr_2.2.1               
[112] MASS_7.3-58.1               rappdirs_0.3.3              car_3.1-0                  
[115] Matrix_1.5-1                cli_3.4.1                   parallel_4.2.1             
[118] igraph_1.3.5                pkgconfig_2.0.3             GenomicAlignments_1.32.1   
[121] foreign_0.8-83              plotly_4.10.0               spatstat.sparse_2.1-1      
[124] foreach_1.5.2               xml2_1.3.3                  stringr_1.4.0              
[127] digest_0.6.29               sctransform_0.3.5           RcppAnnoy_0.0.19           
[130] spatstat.data_2.2-0         leiden_0.4.3                fastmatch_1.1-3            
[133] htmlTable_2.4.1             nortest_1.0-4               dendextend_1.16.0          
[136] uwot_0.1.14                 restfulr_0.0.15             curl_4.3.2                 
[139] shiny_1.7.2                 Rsamtools_2.12.0            rjson_0.2.21               
[142] lifecycle_1.0.1             nlme_3.1-159                jsonlite_1.8.0             
[145] carData_3.0-5               viridisLite_0.4.1           fansi_1.0.3                
[148] pillar_1.8.0                lattice_0.20-45             KEGGREST_1.36.3            
[151] fastmap_1.1.0               httr_1.4.3                  survival_3.4-0             
[154] glue_1.6.2                  FNN_1.1.3.1                 iterators_1.0.14           
[157] png_0.1-7                   glmnet_4.1-4                bit_4.0.4                  
[160] stringi_1.7.8               blob_1.2.3                  latticeExtra_0.6-30        
[163] memoise_2.0.1               irlba_2.3.5.1               future.apply_1.9.1

Thanks

yuchaojiang commented 2 years ago

Please follow the scripts first on here http://htmlpreview.github.io/?https://github.com/yuchaojiang/TRIPOD/blob/main/vignettes/preprocessing_pbmc.html and then on here http://htmlpreview.github.io/?https://github.com/yuchaojiang/TRIPOD/blob/main/vignettes/TRIPOD_pbmc.html Yours is not.

yuchaojiang commented 2 years ago

Specifically, you should use convert = FALSE but I stopped right here. Again, better to follow as is.

DLGisch commented 2 years ago

I missed that. Now it works. Thank you!