scPredict error

[reesyxan]

When running scPredict from the vignette pipeline, I get an error with scPredict:

query <- scPredict(query, reference)
● Matching reference with new dataset...
─ 2000 features present in reference loadings
─ 2000 features shared between reference and new dataset
─ 100% of features in the reference are present in new dataset
Error in GetAssayData():
! assay must be one of "RNA", not "data".
Backtrace:

1. scPred::scPredict(query, reference)
2. scPred::project_query(...)
3. SeuratObject:::GetAssayData.Seurat(assay = "data")

I have not kept up with Seurat updates, did they change assay names recently that might cause this to fail? I verified that GetAssayData(query) does return assay data and it is of type "RNA". Same for reference.

[bioinformaticsperl]

I have encountered the same problem.

[rustlab1]

same issue

[rikfor]

Same issure here

[yossawat-suw]

I also experience the similar problem.
The guest the problem came from the function being called within scPredict which was project_query which in turn call GetAssayData.

thje orginal line within project_query is

    new_data <- GetAssayData(new,  "data")[shared_features, 
            ]

I suppose it want to call layer 'data' from assay 'RNA'. Unfortunately, the new seurat version function GetAssayData is
GetAssayData(object, assay = NULL, layer = NULL, slot = deprecated(), ...)
making the second argument as assay, not the layer as the project_query want

I did work around by editing the project_query function in that specific line to

        new_data <- GetAssayData(new, layer = "data")[shared_features, 
            ]

and assign new function as project_query_edited

Also, I need to edited the scPredict to call project_query_edited instead of project_query. In addition, i need to manually declare function .make_names which will be somewhere during the process

In summary, I create 2 new edited function, and (dont know why I have to) create 1 more function
.make_names

here is the code i run

.make_names <- function(x){
  x <- gsub("\\+", "_plus", x)
  x <- gsub("\\-", "_minus", x)
  x <- make.names(x)
}

scPredict_edited <- function (new, reference, threshold = 0.55, max.iter.harmony = 20, 
    recompute_alignment = TRUE, seed = 66) 
{
    if (!(is(reference, "Seurat") | is(reference, "scPred"))) 
        stop("'object' must be of class 'scPred' or 'Seurat'")
    if (is(reference, "Seurat")) {
        spmodel <- reference@misc$scPred
    }
    else {
        spmodel <- reference
    }
    if (is.null(spmodel)) 
        stop("No feature space has been determined!")
    if (!length(spmodel@train)) 
        stop("No models have been trained!")
    if (!is(new, "Seurat")) 
        stop("New data must be a Seurat object")
    new <- project_query_edited(new, reference = spmodel, max.iter.harmony = max.iter.harmony, 
        recompute_alignment = recompute_alignment, seed = seed)
    new_embeddings_aligned <- Embeddings(new[["scpred"]])
    colnames(new_embeddings_aligned) <- colnames(spmodel@cell_embeddings)
    cellTypeModelNames <- names(spmodel@features)
    .predictCellClass <- function(cellType, spmodel, testEmbeddings) {
        features <- as.character(spmodel@features[[cellType]]$feature)
        model <- spmodel@train[[cellType]]
        prediction <- predict(model, newdata = scPred:::subsetMatrix(testEmbeddings, 
            features), type = "prob")
        rownames(prediction) <- rownames(testEmbeddings)
        prediction[, 1, drop = FALSE]
    }
    cat(crayon::green(cli::symbol$record, " Classifying cells...\n"))
    res <- sapply(cellTypeModelNames, .predictCellClass, spmodel, 
        new_embeddings_aligned)
    res <- as.data.frame(res)
    colnames(res) <- cellTypeModelNames
    rownames(res) <- colnames(new)
    classes <- cellTypeModelNames
    if (length(cellTypeModelNames) == 1) {
        metadata <- get_metadata(spmodel)
        cellClasses <- levels(metadata$pvar)
        res_comp <- 1 - res[, 1]
        negClass <- cellClasses[cellClasses != names(res)]
        res[[negClass]] <- res_comp
    }
    max_props <- as.data.frame(t(apply(res, 1, function(x) c(index = which.max(x), 
        max = x[which.max(x)]))))
    names(max_props) <- c("index", "max")
    max_props$generic_class <- names(res)[max_props$index]
    res <- cbind(res, max_props)
    pred <- ifelse(res$max > threshold, res$generic_class, "unassigned")
    names(pred) <- colnames(new)
    res$prediction <- pred
    res$index <- NULL
    res$no_rejection <- res$generic_class
    res$generic_class <- NULL
    names(res) <- .make_names(paste0("scpred_", names(res)))
    new <- AddMetaData(new, res)
    cat(crayon::green("DONE!\n"))
    new
}

project_query_edited <- function (new, reference, max.iter.harmony = 20, recompute_alignment = TRUE, 
    seed = 66, ...) 
{
    if (!(is(reference, "Seurat") | is(reference, "scPred"))) 
        stop("'object' must be of class 'scPred' or 'Seurat'")
    if (is(reference, "Seurat")) {
        spmodel <- reference@misc$scPred
    }
    else {
        spmodel <- reference
    }
    if (is.null(spmodel)) 
        stop("No feature space has been determined!")
    if (!is(new, "Seurat")) 
        stop("New data must be a Seurat object")
    if ("scpred" %in% names(new@reductions)) {
        if (recompute_alignment) {
            alignment <- TRUE
            cat(crayon::yellow(cli::symbol$figure_dash, "Data has already been aligned to a reference.\n"), 
                sep = "")
            cat(crayon::yellow(cli::symbol$sup_plus, "Skip data alignment using `recompute.alignment = FALSE`.\n"), 
                sep = "")
        }
        else {
            alignment <- FALSE
        }
    }
    else {
        alignment <- TRUE
    }
    if (alignment) {
        cat(crayon::green(cli::symbol$record, " Matching reference with new dataset...\n"))
        ref_loadings <- spmodel@feature_loadings
        ref_embeddings <- spmodel@cell_embeddings
        new_features <- rownames(new)
        reference_features <- rownames(ref_loadings)
        shared_features <- intersect(reference_features, new_features)
        cat(crayon::cyan("\t", cli::symbol$line, paste(length(reference_features), 
            "features present in reference loadings\n")))
        cat(crayon::cyan("\t", cli::symbol$line, paste(length(shared_features), 
            "features shared between reference and new dataset\n")))
        cat(crayon::cyan("\t", cli::symbol$line, paste0(round(length(shared_features)/length(reference_features) * 
            100, 2), "% of features in the reference are present in new dataset\n")))
        ref_loadings <- ref_loadings[shared_features, ]
        new_data <- GetAssayData(new, layer = "data")[shared_features, 
            ]
        means <- spmodel@scaling$means
        stdevs <- spmodel@scaling$stdevs
        new_data <- Matrix::t(new_data)
        names(means) <- names(stdevs) <- rownames(spmodel@scaling)
        means <- means[shared_features]
        stdevs <- stdevs[shared_features]
        i <- stdevs == 0
        if (any(i)) {
            warning(paste0(sum(i), " features have zero variance but are present in the feature loadings. \nDid you subset or integrated this data before?"))
            cat(crayon::yellow("Removing zero-variance genes from projection\n"))
            new_data <- new_data[, !i]
            ref_loadings <- ref_loadings[!i, ]
            means <- means[!i]
            stdevs <- stdevs[!i]
        }
        scaled_data <- scale(new_data, means, stdevs)
        new_embeddings <- scaled_data %*% ref_loadings
        dataset <- factor(c(rep("reference", nrow(ref_embeddings)), 
            rep("new", nrow(new_embeddings))), levels = c("reference", 
            "new"))
        rownames(ref_embeddings) <- paste0("ref_", rownames(ref_embeddings))
        rownames(new_embeddings) <- paste0("new_", rownames(new_embeddings))
        eigenspace <- as.data.frame(rbind(ref_embeddings, new_embeddings))
        meta_data <- data.frame(rownames(eigenspace), dataset = dataset)
        cat(crayon::green(cli::symbol$record, " Aligning new data to reference...\n"))
        set.seed(seed)
        harmony_embeddings <- harmony::HarmonyMatrix(eigenspace, meta_data, 
            "dataset", do_pca = FALSE, reference_values = "reference", 
            max.iter.harmony = max.iter.harmony, ...)
        new_embeddings_aligned <- harmony_embeddings[dataset == 
            "new", , drop = FALSE]
    }
    else {
        new_embeddings_aligned <- Embeddings(new, reduction = "scpred")
        colnames(new_embeddings_aligned) <- gsub("scpred_", spmodel@reduction_key, 
            colnames(new_embeddings_aligned))
    }
    rownames(new_embeddings_aligned) <- gsub("^new_", "", rownames(new_embeddings_aligned))
    new@reductions[["scpred"]] <- CreateDimReducObject(embeddings = new_embeddings_aligned, 
        key = "scpred_", assay = DefaultAssay(object = new))
    if (recompute_alignment) {
        rownames(new_embeddings) <- gsub("^new_", "", rownames(new_embeddings))
        new@reductions[["scpred_projection"]] <- CreateDimReducObject(embeddings = new_embeddings, 
            key = "Projection_", assay = DefaultAssay(object = new))
    }
    new
}

After this it should be fine. At least in my case.

> sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 23.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0

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       

time zone: Europe/London
tzcode source: system (glibc)

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

other attached packages:
[1] scPred_1.9.2

loaded via a namespace (and not attached):
 [1] dotCall64_1.1-1      beeswarm_0.4.0       gtable_0.3.4         spam_2.10-0         
 [5] xfun_0.41            ggplot2_3.4.4        recipes_1.0.9        lattice_0.22-5      
 [9] vctrs_0.6.5          tools_4.3.2          generics_0.1.3       stats4_4.3.2        
[13] parallel_4.3.2       tibble_3.2.1         fansi_1.0.6          pkgconfig_2.0.3     
[17] ModelMetrics_1.2.2.2 Matrix_1.6-4         data.table_1.14.10   lifecycle_1.0.4     
[21] compiler_4.3.2       stringr_1.5.1        munsell_0.5.0        codetools_0.2-19    
[25] SeuratObject_5.0.1   vipor_0.4.7          class_7.3-22         prodlim_2023.08.28  
[29] crayon_1.5.2         pillar_1.9.0         MASS_7.3-60          gower_1.0.1         
[33] iterators_1.0.14     rpart_4.1.23         foreach_1.5.2        nlme_3.1-164        
[37] parallelly_1.36.0    lava_1.7.3           tidyselect_1.2.0     digest_0.6.34       
[41] stringi_1.8.3        future_1.33.1        dplyr_1.1.4          reshape2_1.4.4      
[45] purrr_1.0.2          listenv_0.9.0        splines_4.3.2        cowplot_1.1.2       
[49] grid_4.3.2           harmony_1.2.0        colorspace_2.1-0     cli_3.6.2           
[53] magrittr_2.0.3       survival_3.5-7       utf8_1.2.4           future.apply_1.11.1 
[57] withr_2.5.2          scales_1.3.0         sp_2.1-2             lubridate_1.9.3     
[61] ggbeeswarm_0.7.2     timechange_0.2.0     globals_0.16.2       nnet_7.3-19         
[65] timeDate_4032.109    progressr_0.14.0     pbapply_1.7-2        knitr_1.45          
[69] hardhat_1.3.0        caret_6.0-94         rlang_1.1.3          Rcpp_1.0.12         
[73] glue_1.7.0           pROC_1.18.5          ipred_0.9-14         rstudioapi_0.15.0   
[77] R6_2.5.1             plyr_1.8.9  

[seekning]

after running your function ,it failed aslo

[seekning]

I also experience the similar problem.我也遇到过类似的问题。 The guest the problem came from the function being called within scPredict which was project_query which in turn call GetAssayData.来宾的问题来自在 scPredict 中调用的函数，该函数project_query而 scPredict 又调用 GetAssayData。

thje orginal line within project_query isproject_query内的原始线是

    new_data <- GetAssayData(new,  "data")[shared_features, 
            ]

I suppose it want to call layer 'data' from assay 'RNA'. Unfortunately, the new seurat version function GetAssayData is我想它想从测定“RNA”中调用层“数据”。不幸的是，新的修拉版本函数 GetAssayData 是 GetAssayData(object, assay = NULL, layer = NULL, slot = deprecated(), ...) making the second argument as assay, not the layer as the project_query want将第二个参数作为 Assay，而不是 project_query 想要的层

I did work around by editing the project_query function in that specific line to我确实通过编辑该特定行中的project_query函数来解决

        new_data <- GetAssayData(new, layer = "data")[shared_features, 
            ]

and assign new function as project_query_edited并将新函数指定为project_query_edited

Also, I need to edited the scPredict to call project_query_edited instead of project_query. In addition, i need to manually declare function .make_names which will be somewhere during the process此外，我需要编辑 scPredict 以调用 project_query_edited 而不是 project_query。此外，我需要手动声明函数.make_names该函数将在此过程中的某个地方

In summary, I create 2 new edited function, and (dont know why I have to) create 1 more function总之，我创建了 2 个新的编辑函数，并且（不知道为什么我必须）再创建 1 个函数 .make_names

here is the code i run这是我运行的代码

.make_names <- function(x){
  x <- gsub("\\+", "_plus", x)
  x <- gsub("\\-", "_minus", x)
  x <- make.names(x)
}

scPredict_edited <- function (new, reference, threshold = 0.55, max.iter.harmony = 20, 
    recompute_alignment = TRUE, seed = 66) 
{
    if (!(is(reference, "Seurat") | is(reference, "scPred"))) 
        stop("'object' must be of class 'scPred' or 'Seurat'")
    if (is(reference, "Seurat")) {
        spmodel <- reference@misc$scPred
    }
    else {
        spmodel <- reference
    }
    if (is.null(spmodel)) 
        stop("No feature space has been determined!")
    if (!length(spmodel@train)) 
        stop("No models have been trained!")
    if (!is(new, "Seurat")) 
        stop("New data must be a Seurat object")
    new <- project_query_edited(new, reference = spmodel, max.iter.harmony = max.iter.harmony, 
        recompute_alignment = recompute_alignment, seed = seed)
    new_embeddings_aligned <- Embeddings(new[["scpred"]])
    colnames(new_embeddings_aligned) <- colnames(spmodel@cell_embeddings)
    cellTypeModelNames <- names(spmodel@features)
    .predictCellClass <- function(cellType, spmodel, testEmbeddings) {
        features <- as.character(spmodel@features[[cellType]]$feature)
        model <- spmodel@train[[cellType]]
        prediction <- predict(model, newdata = scPred:::subsetMatrix(testEmbeddings, 
            features), type = "prob")
        rownames(prediction) <- rownames(testEmbeddings)
        prediction[, 1, drop = FALSE]
    }
    cat(crayon::green(cli::symbol$record, " Classifying cells...\n"))
    res <- sapply(cellTypeModelNames, .predictCellClass, spmodel, 
        new_embeddings_aligned)
    res <- as.data.frame(res)
    colnames(res) <- cellTypeModelNames
    rownames(res) <- colnames(new)
    classes <- cellTypeModelNames
    if (length(cellTypeModelNames) == 1) {
        metadata <- get_metadata(spmodel)
        cellClasses <- levels(metadata$pvar)
        res_comp <- 1 - res[, 1]
        negClass <- cellClasses[cellClasses != names(res)]
        res[[negClass]] <- res_comp
    }
    max_props <- as.data.frame(t(apply(res, 1, function(x) c(index = which.max(x), 
        max = x[which.max(x)]))))
    names(max_props) <- c("index", "max")
    max_props$generic_class <- names(res)[max_props$index]
    res <- cbind(res, max_props)
    pred <- ifelse(res$max > threshold, res$generic_class, "unassigned")
    names(pred) <- colnames(new)
    res$prediction <- pred
    res$index <- NULL
    res$no_rejection <- res$generic_class
    res$generic_class <- NULL
    names(res) <- .make_names(paste0("scpred_", names(res)))
    new <- AddMetaData(new, res)
    cat(crayon::green("DONE!\n"))
    new
}

project_query_edited <- function (new, reference, max.iter.harmony = 20, recompute_alignment = TRUE, 
    seed = 66, ...) 
{
    if (!(is(reference, "Seurat") | is(reference, "scPred"))) 
        stop("'object' must be of class 'scPred' or 'Seurat'")
    if (is(reference, "Seurat")) {
        spmodel <- reference@misc$scPred
    }
    else {
        spmodel <- reference
    }
    if (is.null(spmodel)) 
        stop("No feature space has been determined!")
    if (!is(new, "Seurat")) 
        stop("New data must be a Seurat object")
    if ("scpred" %in% names(new@reductions)) {
        if (recompute_alignment) {
            alignment <- TRUE
            cat(crayon::yellow(cli::symbol$figure_dash, "Data has already been aligned to a reference.\n"), 
                sep = "")
            cat(crayon::yellow(cli::symbol$sup_plus, "Skip data alignment using `recompute.alignment = FALSE`.\n"), 
                sep = "")
        }
        else {
            alignment <- FALSE
        }
    }
    else {
        alignment <- TRUE
    }
    if (alignment) {
        cat(crayon::green(cli::symbol$record, " Matching reference with new dataset...\n"))
        ref_loadings <- spmodel@feature_loadings
        ref_embeddings <- spmodel@cell_embeddings
        new_features <- rownames(new)
        reference_features <- rownames(ref_loadings)
        shared_features <- intersect(reference_features, new_features)
        cat(crayon::cyan("\t", cli::symbol$line, paste(length(reference_features), 
            "features present in reference loadings\n")))
        cat(crayon::cyan("\t", cli::symbol$line, paste(length(shared_features), 
            "features shared between reference and new dataset\n")))
        cat(crayon::cyan("\t", cli::symbol$line, paste0(round(length(shared_features)/length(reference_features) * 
            100, 2), "% of features in the reference are present in new dataset\n")))
        ref_loadings <- ref_loadings[shared_features, ]
        new_data <- GetAssayData(new, layer = "data")[shared_features, 
            ]
        means <- spmodel@scaling$means
        stdevs <- spmodel@scaling$stdevs
        new_data <- Matrix::t(new_data)
        names(means) <- names(stdevs) <- rownames(spmodel@scaling)
        means <- means[shared_features]
        stdevs <- stdevs[shared_features]
        i <- stdevs == 0
        if (any(i)) {
            warning(paste0(sum(i), " features have zero variance but are present in the feature loadings. \nDid you subset or integrated this data before?"))
            cat(crayon::yellow("Removing zero-variance genes from projection\n"))
            new_data <- new_data[, !i]
            ref_loadings <- ref_loadings[!i, ]
            means <- means[!i]
            stdevs <- stdevs[!i]
        }
        scaled_data <- scale(new_data, means, stdevs)
        new_embeddings <- scaled_data %*% ref_loadings
        dataset <- factor(c(rep("reference", nrow(ref_embeddings)), 
            rep("new", nrow(new_embeddings))), levels = c("reference", 
            "new"))
        rownames(ref_embeddings) <- paste0("ref_", rownames(ref_embeddings))
        rownames(new_embeddings) <- paste0("new_", rownames(new_embeddings))
        eigenspace <- as.data.frame(rbind(ref_embeddings, new_embeddings))
        meta_data <- data.frame(rownames(eigenspace), dataset = dataset)
        cat(crayon::green(cli::symbol$record, " Aligning new data to reference...\n"))
        set.seed(seed)
        harmony_embeddings <- harmony::HarmonyMatrix(eigenspace, meta_data, 
            "dataset", do_pca = FALSE, reference_values = "reference", 
            max.iter.harmony = max.iter.harmony, ...)
        new_embeddings_aligned <- harmony_embeddings[dataset == 
            "new", , drop = FALSE]
    }
    else {
        new_embeddings_aligned <- Embeddings(new, reduction = "scpred")
        colnames(new_embeddings_aligned) <- gsub("scpred_", spmodel@reduction_key, 
            colnames(new_embeddings_aligned))
    }
    rownames(new_embeddings_aligned) <- gsub("^new_", "", rownames(new_embeddings_aligned))
    new@reductions[["scpred"]] <- CreateDimReducObject(embeddings = new_embeddings_aligned, 
        key = "scpred_", assay = DefaultAssay(object = new))
    if (recompute_alignment) {
        rownames(new_embeddings) <- gsub("^new_", "", rownames(new_embeddings))
        new@reductions[["scpred_projection"]] <- CreateDimReducObject(embeddings = new_embeddings, 
            key = "Projection_", assay = DefaultAssay(object = new))
    }
    new
}

After this it should be fine. At least in my case.在这之后应该没问题。至少在我的情况下。

> sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 23.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0

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       

time zone: Europe/London
tzcode source: system (glibc)

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

other attached packages:
[1] scPred_1.9.2

loaded via a namespace (and not attached):
 [1] dotCall64_1.1-1      beeswarm_0.4.0       gtable_0.3.4         spam_2.10-0         
 [5] xfun_0.41            ggplot2_3.4.4        recipes_1.0.9        lattice_0.22-5      
 [9] vctrs_0.6.5          tools_4.3.2          generics_0.1.3       stats4_4.3.2        
[13] parallel_4.3.2       tibble_3.2.1         fansi_1.0.6          pkgconfig_2.0.3     
[17] ModelMetrics_1.2.2.2 Matrix_1.6-4         data.table_1.14.10   lifecycle_1.0.4     
[21] compiler_4.3.2       stringr_1.5.1        munsell_0.5.0        codetools_0.2-19    
[25] SeuratObject_5.0.1   vipor_0.4.7          class_7.3-22         prodlim_2023.08.28  
[29] crayon_1.5.2         pillar_1.9.0         MASS_7.3-60          gower_1.0.1         
[33] iterators_1.0.14     rpart_4.1.23         foreach_1.5.2        nlme_3.1-164        
[37] parallelly_1.36.0    lava_1.7.3           tidyselect_1.2.0     digest_0.6.34       
[41] stringi_1.8.3        future_1.33.1        dplyr_1.1.4          reshape2_1.4.4      
[45] purrr_1.0.2          listenv_0.9.0        splines_4.3.2        cowplot_1.1.2       
[49] grid_4.3.2           harmony_1.2.0        colorspace_2.1-0     cli_3.6.2           
[53] magrittr_2.0.3       survival_3.5-7       utf8_1.2.4           future.apply_1.11.1 
[57] withr_2.5.2          scales_1.3.0         sp_2.1-2             lubridate_1.9.3     
[61] ggbeeswarm_0.7.2     timechange_0.2.0     globals_0.16.2       nnet_7.3-19         
[65] timeDate_4032.109    progressr_0.14.0     pbapply_1.7-2        knitr_1.45          
[69] hardhat_1.3.0        caret_6.0-94         rlang_1.1.3          Rcpp_1.0.12         
[73] glue_1.7.0           pROC_1.18.5          ipred_0.9-14         rstudioapi_0.15.0   
[77] R6_2.5.1             plyr_1.8.9  

after running your function ,it failed aslo

[HaReL1]

Easier solution: Install @tfguinan repo with the fixed function:
devtools::install_github("tfguinan/scPred")

[iichelhadi]

I resolved the issue by editing the project_query line 45 from
new_data <- GetAssayData(new, "data")[shared_features, ]
to
new_data <- GetAssayData(new, assay = "RNA", slot = "data")[shared_features, ]

[chisquare09]

I have update with the newest version of scPred but still get the error:
Error in GetAssayData():
! assay must be one of "RNA", not "data".

In my case, I want to create a Seurat object by using data from my .csv file, the Seurat object create is correct and the output is:
`> print(reference)
An object of class Seurat
100 features across 1508 samples within 1 assay
Active assay: RNA (100 features, 100 variable features)
3 layers present: counts, data, scale.data
2 dimensional reductions calculated: pca, umap

print(query)
An object of class Seurat
100 features across 378 samples within 1 assay
Active assay: RNA (100 features, 0 variable features)
2 layers present: counts, data`

Here is my code:
`library(Seurat)
library(scPred)
library(magrittr)
library(readr)
library(Matrix)

Load data

data <- read.csv(file='baron-prepare-log_count_100pca.csv', header=FALSE)
set.seed(123)

data_y <- data[,1]
data_x <- data[,-1]

Train-test split

sample_size <- floor(0.8 * nrow(data_x))
train_indices <- sample(seq_len(nrow(data_x)), size = sample_size)

reference_x <- data_x[train_indices, ]
reference_y <- data_y[train_indices]

query_x <- data_x[-train_indices, ]
query_y <- data_y[-train_indices]

Convert data to sparse matrix format

reference_x <- as(t(reference_x), "dgCMatrix")
query_x <- as(t(query_x), "dgCMatrix")

Create Seurat objects

reference <- CreateSeuratObject(counts = reference_x, assay = "RNA")
query <- CreateSeuratObject(counts = query_x, assay = "RNA")

Assign proper cell names

colnames(reference) <- paste0("RefCell_", seq_len(ncol(reference)))
colnames(query) <- paste0("QueryCell_", seq_len(ncol(query)))

Add cell type metadata

reference$cell_type <- as.factor(reference_y)
query$cell_type <- as.factor(query_y)

Manually copy "counts" to "data" instead of running NormalizeData()

counts_matrix <- GetAssayData(reference,assay = "RNA" ,layer = "counts")
reference <- SetAssayData(reference,assay = "RNA", layer = "data", new.data = counts_matrix)

Feature selection & scaling

VariableFeatures(reference) <- rownames(reference) # Use all genes
reference <- ScaleData(reference, features = rownames(reference))

Run PCA

reference <- RunPCA(reference, features = rownames(reference))

Run UMAP

reference <- RunUMAP(reference, dims = 1:30)

View Seurat object summary

print(reference)

Plot UMAP

DimPlot(reference, group.by = "cell_type", label = TRUE, repel = TRUE)

Train scPred classifier

reference <- getFeatureSpace(reference, "cell_type")
reference <- trainModel(reference)
get_probabilities(reference) %>% head()

get_scpred(reference)

Manually copy "counts" to "data" for query

query_counts_matrix <- GetAssayData(query, assay = "RNA",layer = "counts")
query <- SetAssayData(query, assay = "RNA",layer = "data", new.data = query_counts_matrix)

Run scPred

query <- scPredict(query, reference)

View Seurat object summary

print(query)

Plot scPred results

DimPlot(query, group.by = "scpred_prediction", reduction = "scpred")`

[hk20013106]

Easier solution: Install @tfguinan repo with the fixed function: devtools::install_github("tfguinan/scPred")

thank you so much!