scpModel workflow/readSCP

[edemmott]

So, if I use the example data (Leduc_minimal), the modelling workflow works as described in the vignette. However when I try it with my own data it fails.

One difference is when I use readSCP on my own data, it seems to be generating a QFeatures object, whereas the Leduc_minimal dataset is recorded as Large SingleCellExperiment.

Running the standard SCoPE2 workflow (as described in that vignette) works on my data and doesnt mind it being described as a QFeatures rather than SCE object, it crashes at the point of using 'scpModelWorkflow', giving the error:

Error in scpModel(*tmp*, name) :
inherits(object, "SummarizedExperiment") is not TRUE

Environment:
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.5

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0

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

time zone: Europe/London
tzcode source: internal

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

other attached packages:
[1] scpdata_1.13.0 ExperimentHub_2.13.1 AnnotationHub_3.13.3 BiocFileCache_2.13.0
[5] dbplyr_2.5.0 scran_1.33.2 scater_1.33.4 scuttle_1.15.4
[9] SingleCellExperiment_1.27.2 sva_3.53.0 BiocParallel_1.39.0 genefilter_1.87.0
[13] mgcv_1.9-1 nlme_3.1-166 impute_1.79.0 patchwork_1.2.0
[17] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[21] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[25] ggplot2_3.5.1 tidyverse_2.0.0 scp_1.15.1 BiocManager_1.30.25
[29] QFeatures_1.15.2 MultiAssayExperiment_1.31.5 SummarizedExperiment_1.35.1 Biobase_2.65.1
[33] GenomicRanges_1.57.1 GenomeInfoDb_1.41.1 IRanges_2.39.2 S4Vectors_0.43.2
[37] BiocGenerics_0.51.1 MatrixGenerics_1.17.0 matrixStats_1.4.1

loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 rstudioapi_0.16.0 jsonlite_1.8.8 magrittr_2.0.3
[5] ggbeeswarm_0.7.2 farver_2.1.2 zlibbioc_1.51.1 vctrs_0.6.5
[9] memoise_2.0.1 BiocBaseUtils_1.7.3 S4Arrays_1.5.7 curl_5.2.2
[13] BiocNeighbors_1.99.0 SparseArray_1.5.34 plyr_1.8.9 cachem_1.1.0
[17] igraph_2.0.3 lifecycle_1.0.4 pkgconfig_2.0.3 rsvd_1.0.5
[21] Matrix_1.7-0 R6_2.5.1 fastmap_1.2.0 GenomeInfoDbData_1.2.12
[25] clue_0.3-65 fdrtool_1.2.18 colorspace_2.1-1 AnnotationDbi_1.67.0
[29] dqrng_0.4.1 irlba_2.3.5.1 lpsymphony_1.33.1 RSQLite_2.3.7
[33] beachmat_2.21.6 filelock_1.0.3 labeling_0.4.3 fansi_1.0.6
[37] timechange_0.3.0 httr_1.4.7 abind_1.4-8 compiler_4.4.1
[41] bit64_4.0.5 withr_3.0.1 viridis_0.6.5 DBI_1.2.3
[45] MASS_7.3-61 rappdirs_0.3.3 DelayedArray_0.31.11 bluster_1.15.1
[49] tools_4.4.1 vipor_0.4.7 beeswarm_0.4.0 glue_1.7.0
[53] grid_4.4.1 Rtsne_0.17 cluster_2.1.6 reshape2_1.4.4
[57] generics_0.1.3 gtable_0.3.5 tzdb_0.4.0 hms_1.1.3
[61] BiocSingular_1.21.3 ScaledMatrix_1.13.0 metapod_1.13.0 utf8_1.2.4
[65] XVector_0.45.0 BiocVersion_3.20.0 ggrepel_0.9.6 pillar_1.9.0
[69] limma_3.61.9 splines_4.4.1 lattice_0.22-6 survival_3.7-0
[73] bit_4.0.5 annotate_1.83.0 tidyselect_1.2.1 locfit_1.5-9.10
[77] Biostrings_2.73.1 gridExtra_2.3 ProtGenerics_1.37.1 edgeR_4.3.14
[81] IHW_1.33.0 statmod_1.5.0 stringi_1.8.4 UCSC.utils_1.1.0
[85] yaml_2.3.10 lazyeval_0.2.2 codetools_0.2-20 nipals_0.8
[89] MsCoreUtils_1.17.1 cli_3.6.3 xtable_1.8-4 munsell_0.5.1
[93] Rcpp_1.0.13 png_0.1-8 XML_3.99-0.17 parallel_4.4.1
[97] blob_1.2.4 AnnotationFilter_1.29.0 viridisLite_0.4.2 slam_0.1-53
[101] scales_1.3.0 crayon_1.5.3 rlang_1.1.4 cowplot_1.1.3
[105] KEGGREST_1.45.1

[edemmott]

Ok. So looking more at it, the QFeatures objects usually created contain the SCE as individual assays. Probably useful to change the example for the vignette over to the QFeatures object?

[lgatto]

One difference is when I use readSCP on my own data, it seems to be generating a QFeatures object, whereas the Leduc_minimal dataset is recorded as Large SingleCellExperiment.

All the data management from loading to QC, joining, aggregation, ... is indeed handled by QFeatures objects to keep track of all the steps. The modelling only requires the last set of the QFeatures object, that one can extract with getWithColData().

Probably useful to change the example for the vignette over to the QFeatures object?

Not sure I understand what you mean? readSCP() is a version that created a QFeatures object containing SingleCellExperiment sets. The readQFeatures() version documented in the QFeatures package returns a QFeatures object containing SummarizedExperiment sets.

[edemmott]

Sorry, I phrased badly. I think this is because the Leduc example dataset is the extracted single-cell experiment. So the vignette for the scpModel workflow starts from the singlecellexperiment leduc part. Wheres the more complete scp reproduction example using this dataset makes it clear you need to use getwithcoldata to extract the singlecellexperiment from the QFeatures object. The current vignete is just a little incomplete. Perhaps transfer over more of the scp repdocution example?

[lgatto]

Ok, thanks. That indeed makes sense, as some users might come straight to that vignette.

Instead of transferring more code (which leads to more maintenance and a longer, less straight-to-the-point vignette), I will add a sentence or to clarifying that we start immediately from the SingleCellExperiment, and point to the other vignette for more context.

[michabirklbauer]

Sorry to revive an old discussion, but maybe to add to that: Is there any documentation that compares QFeatures, SummarizedExperiment and SingleCellExperiment? They pretty much seem like (almost) the same data structure to me, but they can't simply be type cast between each other?

[leopoldguyot]

Sorry to revive an old discussion, but maybe to add to that: Is there any documentation that compares QFeatures, SummarizedExperiment and SingleCellExperiment? They pretty much seem like (almost) the same data structure to me, but they can't simply be type cast between each other?

Hello, thank you for your interest in scp/QFeatures. Here is a short explanation of the relationship between the 3 classes.

SingleCellExperiment (SCE) extends SummarizedExperiment (SE), so an SCE is essentially a specialized SE with extra single-cell specific slots, such as reduced dimensions, alternative experiments, ...

QFeatures is different: it extends MultiAssayExperiment, not SummarizedExperiment. A QFeatures object contains multiple SummarizedExperiment or SingleCellExperiment objects, each one called a set. Its main benefit is that it can keep track of the links between features across these different sets, for example PSMs to peptides to proteins.

You could technically create a QFeatures object containing only one SE, but that wouldn't be really useful because the main added value is the relationship between multiple sets. If you want to work with one specific set inside a QFeatures object, you can extract it as an SE or SCE with getWithColData().

The vignettes that could probably help you to understand the QFeatures class are scp: QFeatures in a nutshell, scp: Load data using readSCP, and QFeatures: Load data using readQFeatures().

[michabirklbauer]

Thanks a lot for the explanation!

I think what I am having a hard time with wrapping my head around is, that QFeatures has assays which are SummarizedExperiments or SingleCellExperiments - but both of these classes also have assays. At least in theory they support more than one assay per SE/SCE, but I guess when working with QFeatures each QFeatures assay only contains one SE/SCE with exactly one SE/SCE assay? Otherwise the whole data structure essentially becomes 4-dimensional.

Which boils down to QFeatures is the same data structure/different class as SE but with the added benefit of being able to link assays to keep track of different transformation levels - if I understood that correctly?

[leopoldguyot]

I think what I am having a hard time with wrapping my head around is, that QFeatures has assays which are SummarizedExperiments or SingleCellExperiments - but both of these classes also have assays. At least in theory they support more than one assay per SE/SCE, but I guess when working with QFeatures each QFeatures assay only contains one SE/SCE with exactly one SE/SCE assay? Otherwise the whole data structure essentially becomes 4-dimensional.

There is indeed one terminology issue that makes this quite confusing.

In SummarizedExperiment, an assay is a matrix stored inside the SE. In QFeatures, the word “assay” was historically used for
what is now better described as a set: one SummarizedExperiment or SingleCellExperiment inside the QFeatures object. That overlap is confusing, which is why we generally try to use “set” for the SE/SCE objects contained in a QFeatures object. Some older terminology still remains in function names such as joinAssays(), replaceAssay(), etc.

When an SE/SCE is used as a set inside a QFeatures object, only the first assay of that SE/SCE is considered.

Which boils down to QFeatures is the same data structure/different class as SE but with the added benefit of being able to link assays to keep track of different transformation levels - if I understood that correctly?

The reason we don't represent the data as a single SE with multiple assays is that the dimensions often change between processing steps. For example, filtering can remove rows, and aggregation changes the feature level entirely. Since assays inside one SE are expected to have the same dimensions, those different processing levels cannot naturally live as multiple assays in a single SE. The rowData would also be different after an aggregation step.

To summarize it, a QFeatures is a container for multiple SE/SCE sets (historically named assays), where each set represents a different
processing step (QFeatures operations only take into account the first assay of each set), and QFeatures adds the infrastructure to track the relationships between features across those sets.

[michabirklbauer]

Ahhh yes, I failed to consider that dimensions/row data could change in QFeatures - that makes a lot of sense then!

When an SE/SCE is used as a set inside a QFeatures object, only the first assay of that SE/SCE is considered.

Good to know!

I think I have it figured out now, thanks a lot for the thorough explanation! ☺️