error using tx_file in plotGeometry

[heziqing]

Hello, I used following code to filter molecules by cell type. And use plotGeometry to plot the filtered molecules, but it encountered an error.

error
Error in .get_tx_df(sfe, data_dir = NULL, tech = NULL, file = NULL, sample_id = sample_id, : GDAL Parquet driver is required to selectively read genes from GeoParquet file.

how to filter molecules

dir_use <- dir_list[[tmpsample]]
mols <- fread(paste0(dir_use,"/detected_transcripts.csv"))
  
tmp_cells <- read.csv(Bigcluster_ID[[tmpsample]],colClasses = "character")
tmpmols <-  mols[mols$cell_id %in% tmp_cells$cell_ID, ]
fwrite(tmpmols, file = paste0(outdir, tmpsample,"_filtered_transcripts.csv"), row.names = FALSE)

plot

tmp_mer <- readObject(paste0("/dath/heziqing/Vessels/Merfish/SFE/saveObject/",tmpsample,"_filter5_volume20")) 
tmptx <- formatTxSpots(paste0("/dath/heziqing/Vessels/Merfish/Gene_plot/filter_molecules/",tmpsample,"_filtered_transcripts.csv"))
p <- plotGeometry(tmp_mer,gene=tmpgene,colGeometryName = "cellSeg",rowGeometryName="txSpots",
                fill=FALSE,dark=FALSE,
                tx_size=tmpsize,
                tx_file = tmptx,
                border.size=0.03,bbox=box_list[[tmpsample]],tx_col = "darkred")

Here is the sessionInfo. Thanks in advances!

R version 4.4.2 (2024-10-31)
Platform: x86_64-conda-linux-gnu
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /share/home/heziqing/miniconda3/envs/SFE_devel/lib/libopenblasp-r0.3.28.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
 [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
[11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

time zone: Asia/Shanghai
tzcode source: system (glibc)

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

other attached packages:
 [1] readr_2.1.5                    gridExtra_2.3                 
 [3] data.table_1.16.4              zeallot_0.1.0                 
 [5] rlang_1.1.5                    terra_1.8-21                  
 [7] SpatialExperiment_1.17.0       alabaster.sfe_0.99.2001       
 [9] alabaster.base_1.6.1           zellkonverter_1.16.0          
[11] BiocNeighbors_2.0.1            BiocParallel_1.40.0           
[13] stringr_1.5.1                  patchwork_1.3.0               
[15] scuttle_1.16.0                 ggplot2_3.5.1                 
[17] lifecycle_1.0.4                RBioFormats_1.6.0             
[19] Matrix_1.6-5                   rjson_0.2.23                  
[21] spdep_1.3-10                   spData_2.3.3                  
[23] fs_1.6.5                       sf_1.0-19                     
[25] SingleCellExperiment_1.29.1    SummarizedExperiment_1.36.0   
[27] Biobase_2.66.0                 GenomicRanges_1.58.0          
[29] GenomeInfoDb_1.42.3            IRanges_2.40.1                
[31] S4Vectors_0.44.0               BiocGenerics_0.53.6           
[33] generics_0.1.3                 MatrixGenerics_1.18.1         
[35] matrixStats_1.5.0              SFEData_1.8.1                 
[37] Voyager_1.8.1                  SpatialFeatureExperiment_1.9.8

loaded via a namespace (and not attached):
  [1] splines_4.4.2             bitops_1.0-9             
  [3] filelock_1.0.3            tibble_3.2.1             
  [5] R.oo_1.27.0               basilisk.utils_1.18.0    
  [7] edgeR_4.4.2               lattice_0.22-6           
  [9] MASS_7.3-64               magrittr_2.0.3           
 [11] limma_3.62.2              yaml_2.3.10              
 [13] sp_2.2-0                  reticulate_1.40.0        
 [15] DBI_1.2.3                 multcomp_1.4-28          
 [17] abind_1.4-8               spatialreg_1.3-6         
 [19] zlibbioc_1.52.0           purrr_1.0.4              
 [21] R.utils_2.12.3            RCurl_1.98-1.16          
 [23] TH.data_1.1-3             rappdirs_0.3.3           
 [25] sandwich_3.1-1            GenomeInfoDbData_1.2.13  
 [27] alabaster.sce_1.7.0       units_0.8-5              
 [29] RSpectra_0.16-2           dqrng_0.4.1              
 [31] DelayedMatrixStats_1.28.1 codetools_0.2-20         
 [33] DropletUtils_1.26.0       DelayedArray_0.33.5      
 [35] xml2_1.3.6                tidyselect_1.2.1         
 [37] UCSC.utils_1.2.0          memuse_4.2-3             
 [39] farver_2.1.2              BiocFileCache_2.14.0     
 [41] jsonlite_1.8.9            e1071_1.7-16             
 [43] survival_3.8-3            tools_4.4.2              
 [45] ggnewscale_0.5.0          sfarrow_0.4.1            
 [47] Rcpp_1.0.14               glue_1.8.0               
 [49] SparseArray_1.7.5         EBImage_4.48.0           
 [51] dplyr_1.1.4               HDF5Array_1.34.0         
 [53] alabaster.spatial_1.7.1   withr_3.0.2              
 [55] BiocManager_1.30.25       fastmap_1.2.0            
 [57] basilisk_1.18.0           boot_1.3-31              
 [59] rhdf5filters_1.18.0       bluster_1.16.0           
 [61] digest_0.6.37             R6_2.5.1                 
 [63] colorspace_2.1-1          wk_0.9.4                 
 [65] LearnBayes_2.15.1         jpeg_0.1-10              
 [67] RSQLite_2.3.9             R.methodsS3_1.8.2        
 [69] class_7.3-23              httr_1.4.7               
 [71] htmlwidgets_1.6.4         S4Arrays_1.6.0           
 [73] pkgconfig_2.0.3           scico_1.5.0              
 [75] rJava_1.0-11              gtable_0.3.6             
 [77] blob_1.2.4                XVector_0.46.0           
 [79] htmltools_0.5.8.1         fftwtools_0.9-11         
 [81] alabaster.matrix_1.7.6    scales_1.3.0             
 [83] png_0.1-8                 tzdb_0.4.0               
 [85] coda_0.19-4.1             nlme_3.1-165             
 [87] curl_6.2.0                proxy_0.4-27             
 [89] cachem_1.1.0              zoo_1.8-12               
 [91] rhdf5_2.50.2              BiocVersion_3.20.0       
 [93] KernSmooth_2.23-26        parallel_4.4.2           
 [95] arrow_18.1.0.1            AnnotationDbi_1.68.0     
 [97] s2_1.1.7                  pillar_1.10.1            
 [99] alabaster.schemas_1.2.0   vctrs_0.6.5              
[101] dbplyr_2.5.0              beachmat_2.22.0          
[103] sfheaders_0.4.4           cluster_2.1.8            
[105] magick_2.8.5              mvtnorm_1.3-2            
[107] cli_3.6.3                 locfit_1.5-9.10          
[109] compiler_4.4.2            crayon_1.5.3             
[111] classInt_0.4-10           stringi_1.8.4            
[113] alabaster.se_1.7.0        deldir_2.0-4             
[115] assertthat_0.2.1          munsell_0.5.1            
[117] Biostrings_2.74.1         tiff_0.1-12              
[119] dir.expiry_1.14.0         ExperimentHub_2.14.0     
[121] hms_1.1.3                 sparseMatrixStats_1.18.0 
[123] bit64_4.6.0-1             Rhdf5lib_1.28.0          
[125] KEGGREST_1.46.0           statmod_1.5.0            
[127] alabaster.ranges_1.2.0    AnnotationHub_3.14.0     
[129] igraph_2.1.4              memoise_2.0.1            
[131] bit_4.5.0.1              

[lambdamoses]

This is because the Parquet GDAL driver is required to selectively read transcript spots for one gene. The default GDAL binary is quite old and does not have the Parquet driver. To make sure that Parquet driver is available on Mac, you can brew install arrow before brew install gdal. On Linux, it's more complicated since you need to compile GDAL from source. See https://github.com/rocker-org/rocker-versioned2/blob/master/scripts/experimental/install_dev_osgeo.sh Looks like you're on a server, so compiling may be challenging. You can use this Docker container which has GDAL with the Parquet driver: docker pull rocker/geospatial:dev-osgeo. On a server, you can use Singularity of Apptainer to use Docker images. If you haven't done so before, here's a guide: https://apptainer.org/docs/user/main/docker_and_oci.html https://bioconductor.org/help/docker/

Or for now you can use arrow::open_dataset() to open the parquet file without loading it into memory. Something like this

library(arrow)
library(tidyverse)
ds <- open_dataset("tx.parquet")
df <- ds |> filter(gene == "EPCAM") |> collect()

Then df has a geometry column which is a WKB binary. Convert it to sfc with sf::st_as_sfc().

Looks like that's kind of annoying, so long term solution would be partitioning the GeoParquet file, which solves the capacity issue as well and do away withe the GDAL Parquet driver. I'll try reimplementing parquet-related functionalities with open_dataset() and partitioning this week.

[heziqing]

Thanks for your reply. I used addTxSpots function to solve my problem.

dir_use <- dir_list[[tmpsample]]
tmp_mer <- readVizgen(dir_use, z = layer_list[[tmpsample]], image = "DAPI", add_molecules = FALSE,flip = "none",
                        min_area = 1,file_out=NULL,
                     BPPARAM=BiocParallel::MulticoreParam(20, tasks = 20L,force.GC = FALSE, progressbar = TRUE)) 

 tmp_mer <- addTxSpots(tmp_mer,
                        file= paste0("/dath/heziqing/Vessels/Merfish/Gene_plot/filter_molecules/",
                                     tmpsample,"_filtered_transcripts_EC.csv"))

But I want to plot a ficture just like this:

Here is the code I run, but it didnot work.

Allcols <- c("#FB0007","#23FF27","#F9FF0D","#14A1F3","#F500FF","#7F01FF","#FC4B0B","#22FFFF")
Allmarkers <- c("MFSD2A","ABCG1","ABCG2","SLC1A4","SLC7A1","SLC7A5","SLC38A3","SLCO2A1")
names(Allcols) <- Allmarkers
allgenes <- rowGeometries(tmp_mer)[[1]]$gene
tmpmarkers <- intersect(Allmarkers,allgenes)
tmpcols <- Allcols[tmpmarkers]
tmpsize <- 0.5

pdf(paste0(outdir,tmpsample,"_plot_Final.pdf"),width = 7,height = 7)
 p <- plotSpatialFeature2(tmp_mer, "cluster", colGeometryName = "cellSeg", 
                          shape = 21, alpha = 0.5,
                          rowGeometryName="txSpots",rowGeometryFeatures=tmpmarkers,
                          size=tmpsize,tx_fixed=list(shape=21,size=tmpsize,stroke=0),
                          dark = FALSE,bbox=box_list[[tmpsample]],tx_cols=tmpcols)+
      scale_fill_manual(values=cluster_cols[levels(tmp_mer$cluster)])
 print(p)
dev.off()

I have changed the function plot_var_sf like this

.plot_var_sf2 <- function(df, annot_df, img_df, channel, type, type_annot, feature_aes,
                         feature_fixed, annot_aes, annot_fixed, tx_fixed, divergent,
                         diverge_center,annot_divergent, annot_diverge_center,
                         ncol_sample, scattermore, pointsize,
                         bins, summary_fun, hex, maxcell, show_axes, dark, palette,
                         tx_df,tx_cols) {
    # Add annotGeometry if present
    if (!is.null(annot_df)) {
        annot_fixed <- .get_applicable(type_annot, annot_fixed)
        annot_fixed <- annot_fixed[setdiff(names(annot_fixed), names(annot_aes))]
        if ("color" %in% names(annot_aes) && annot_fixed$linewidth == 0) {
            annot_fixed$linewidth <- 0.5
        }
        geom_annot <- do.call(geom_sf, c(
            list(mapping = aes(!!!syms(annot_aes)), data = annot_df),
            annot_fixed
        ))
        pal_annot <- .get_pal(annot_df, annot_aes, 2, annot_divergent,
                              annot_diverge_center, dark = dark)
    }
    # Deal with gene symbols that are not legal R object names
    df_names_orig <- names(df)
    names(df) <- make.names(names(df))
    feature_aes <- lapply(feature_aes, make.names)
    inds <- !df_names_orig %in% names(df)
    ind_plot <- which(names(df) == unlist(feature_aes)) # should be only one
    if (any(inds) && length(ind_plot))
        name_show <- df_names_orig[ind_plot] else name_show <- waiver()

    p <- ggplot()
    data <- NULL
    if (!is.null(img_df)) {
        # Check if it's RGB
        img <- img_df$data[[1]]
        img_dim <- dim(img)
        p <- p + geom_spi_rgb(data = img_df, aes(spi = data),
                              palette = palette)
    }
    # Filled polygon annotations go beneath feature plot
    is_annot_filled <- !is.null(annot_df) &&
        ("fill" %in% names(c(annot_aes, annot_fixed))) &&
        type_annot %in% c("POLYGON", "MULTIPOLYGON")
    if ("fill" %in% names(annot_fixed)) {
        is_annot_filled <- is_annot_filled && !is.na(annot_fixed[["fill"]])
    }
    if (is_annot_filled) {
        p <- p + geom_annot
        if (!is.null(pal_annot)) p <- p + pal_annot
    }

    feature_fixed <- .get_applicable(type, feature_fixed)
    feature_fixed <- feature_fixed[setdiff(names(feature_fixed),
                                           names(feature_aes))]

    if ("fill" %in% names(feature_aes) && is_annot_filled) {
        p <- p + new_scale_fill()
    }
    if (scattermore) {
        geom_use <- do.call(scattermore::geom_scattermore,
                            c(list(mapping = aes(!!!syms(c(list(x = "X", y = "Y"),
                                                           feature_aes))),
                                   data = df, pointsize = pointsize),
                              feature_fixed))
    } else if (!is.null(bins)) {
        names(feature_aes)[names(feature_aes) == "color"] <- "z"
        name_show <- paste("Aggregated", feature_aes[["z"]], sep = "\n")
        feature_fixed <- feature_fixed["alpha"]
        stat_fun <- if (hex) stat_summary_hex else stat_summary_2d
        geom_use <- do.call(stat_fun,
                            c(list(mapping = aes(!!!syms(c(list(x = "X", y = "Y"),
                                                           feature_aes))),
                                   data = df, bins = bins, fun = summary_fun),
                              feature_fixed))
    } else {
        geom_use <- do.call(geom_sf, c(
            list(mapping = aes(!!!syms(feature_aes)), data = df),
            feature_fixed
        ))
    }
    p <- p + geom_use
    if (scattermore || !is.null(bins)) p <- p + coord_equal()
    # Palette
    pal <- .get_pal(df, feature_aes, 1, divergent, diverge_center,
                    name = name_show, dark = dark)
    if (!is.null(pal)) p <- p + pal

    # Line and point annotations go above feature plot
    if (!is.null(annot_df) && !is_annot_filled) {
        if (!is.null(pal_annot) && "color" %in% names(feature_aes)) {
            p <- p + new_scale_color()
        }
        p <- p + geom_annot
        if (!is.null(pal_annot)) {
            p <- p + pal_annot
        }
    }
    if (!is.null(tx_df)) {
        tx_fixed <- .get_applicable("MULTIPOINT", tx_fixed)
        if (length(unique(tx_df$gene)) > 1L) {
            # tx_fixed <- tx_fixed[names(tx_fixed) != "shape"]
            # geom_tx <- do.call(geom_sf, c(
            #     list(mapping = aes(shape = gene), data = tx_df),
            #     tx_fixed
            # )
             geom_tx <- list(
              do.call(geom_sf, c(
                list(mapping = aes(color = gene,fill=gene), data = tx_df),
                tx_fixed
              )), 
             scale_fill_manual(values = tx_cols)  
              scale_color_manual(values = tx_cols)  
            )
            
        } else geom_tx <- do.call(geom_sf, c(
            list(data = tx_df), tx_fixed
        ))
            p <- p + geom_tx
    }
    if ("sample_id" %in% names(df) && length(unique(df$sample_id)) > 1L) {
        p <- p +
            facet_wrap(~sample_id, ncol = ncol_sample)
    }

    if (dark) {
        p <- p + .dark_theme(show_axes)
    } else {
        if (show_axes) p <- p + theme_bw() else p <- p + theme_void()
    }
    p
}

Do you have any advices? Thanks in advance!