A detailed explanation of scFeatures’ features
18 January 2024
Contents
1 Introduction
scFeatures is a tool for generating multi-view representations of samples in a single-cell dataset. This vignette provides an overview of scFeatures. It uses the main function to generate features and then illustrates case studies of using the generated features for classification, survival analysis and association study.
library(scFeatures)2 Running scFeatures
scFeatures can be run using one line of code scfeatures_result <- scFeatures(data), which generates a list of dataframes containing all feature types in the form of samples x features.
data("example_scrnaseq" , package = "scFeatures")
data <- example_scrnaseq
scfeatures_result <- scFeatures(data = data@assays$RNA@data, sample = data$sample, celltype = data$celltype,
feature_types = "gene_mean_celltype" ,
type = "scrna",
ncores = 1,
species = "Homo sapiens")By default, the above function generates all feature types. To reduce the computational time for the demonstrate, here we generate only the selected feature type “gene mean celltype”. More information on the function customisation can be obtained by typing ?scFeatures()
2.1 Classification of conditions using the generated features
To build disease prediction model from the generated features we utilise ClassifyR.
The output from scFeatures is a matrix of sample x feature, ie, the row corresponds to each sample, the column corresponds to the feature, and can be directly used as the X. The order of the rows is in the order of unique(data$sample).
Here we use the feature type gene mean celltype as an example to build classification model on the disease condition.
feature_gene_mean_celltype <- scfeatures_result$gene_mean_celltype
# inspect the first 5 rows and first 5 columns
feature_gene_mean_celltype[1:5, 1:5]
#> Naive T Cells--SNRPD2 Naive T Cells--NOSIP Naive T Cells--IL2RG
#> Pre_P8 1.6774074 1.856023 1.834704
#> Pre_P6 3.4815573 3.217231 3.583760
#> Pre_P27 0.0000000 1.486346 1.742069
#> Pre_P7 1.1761242 1.282083 2.050024
#> Pre_P20 0.6999054 1.315393 2.781787
#> Naive T Cells--ALDOA Naive T Cells--LUC7L3
#> Pre_P8 1.598921 1.8056758
#> Pre_P6 3.493135 0.0000000
#> Pre_P27 3.518687 1.5080699
#> Pre_P7 2.575957 0.9163035
#> Pre_P20 1.626403 1.2659969
# inspect the dimension of the matrix
dim(feature_gene_mean_celltype)
#> [1] 19 2217We recommend using ClassifyR::crossValidate to do cross-validated classification with the extracted feaures.
library(ClassifyR)
# X is the feature type generated
# y is the condition for classification
X <- feature_gene_mean_celltype
y <- data@meta.data[!duplicated(data$sample), ]
y <- y[match(rownames(X), y$sample), ]$condition
# run the classification model using random forest
result <- ClassifyR::crossValidate(
X, y,
classifier = "randomForest", nCores = 2,
nFolds = 3, nRepeats = 5
)
ClassifyR::performancePlot(results = result)
It is expected that the classification accuracy is low. This is because we are using a small subset of data containing only 3523 genes and 519 cells. The dataset is unlikely to contain enough information to distinguish responders and non-responders.
2.2 Survival analysis using the generated features
Suppose we want to use the features to perform survival analysis. In here, since the patient outcomes are responder and non-responder, and do not contain survival information, we randomly “generate” the survival outcome for the purpose of demonstration.
We use a standard hierarchical clustering to split the patients into 2 groups based on the generated features.
library(survival)
library(survminer)
X <- feature_gene_mean_celltype
X <- t(X)
# run hierarchical clustering
hclust_res <- hclust(
as.dist(1 - cor(X, method = "pearson")),
method = "ward.D2"
)
set.seed(1)
# generate some survival outcome, including the survival days and the censoring outcome
survival_day <- sample(1:100, ncol(X))
censoring <- sample(0:1, ncol(X), replace = TRUE)
cluster_res <- cutree(hclust_res, k = 2)
metadata <- data.frame( cluster = factor(cluster_res),
survival_day = survival_day,
censoring = censoring)
# plot survival curve
fit <- survfit(
Surv(survival_day, censoring) ~ cluster,
data = metadata
)
ggsurv <- ggsurvplot(fit,
conf.int = FALSE, risk.table = TRUE,
risk.table.col = "strata", pval = TRUE
)
ggsurv
The p-value is very high, indicating there is not enough evidence to claim there is a survival difference between the two groups. This is as expected, because we randomly assigned survival status to each of the patient.
3 Association study of the features with the conditions
scFeatures provides a function that automatically run association study of the features with the conditions and produce an HTML file with the visualisation of the features and the association result.
For this, we would first need to generate the features using scFeatures and then store the result in a named list format.
For demonstration purpose, we provide an example of this features list. The code below show the steps of generating the HTML output from the features list.
# here we use the demo data from the package
data("scfeatures_result" , package = "scFeatures")
# here we use the current working directory to save the html output
# modify this to save the html file to other directory
output_folder <- tempdir()
run_association_study_report(scfeatures_result, output_folder )
#> /usr/bin/pandoc +RTS -K512m -RTS output_report.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpQ6QhnW/output_report.html --lua-filter /home/biocbuild/bbs-3.18-bioc/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /home/biocbuild/bbs-3.18-bioc/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --table-of-contents --toc-depth 3 --variable toc_float=1 --variable toc_selectors=h1,h2,h3 --variable toc_collapsed=1 --variable toc_print=1 --template /home/biocbuild/bbs-3.18-bioc/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --number-sections --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpQ6QhnW/rmarkdown-str38ba7c662311ac.html --variable code_folding=hide --variable code_menu=1Inside the directory defined in the output_folder, you will see the html report output with the name output_report.html.
4 sessionInfo()
sessionInfo()
#> R version 4.3.2 Patched (2023-11-13 r85521)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] grid stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] data.table_1.14.10 enrichplot_1.22.0
#> [3] DOSE_3.28.2 clusterProfiler_4.10.0
#> [5] msigdbr_7.5.1 EnsDb.Hsapiens.v79_2.99.0
#> [7] ensembldb_2.26.0 AnnotationFilter_1.26.0
#> [9] GenomicFeatures_1.54.1 org.Hs.eg.db_3.18.0
#> [11] AnnotationDbi_1.64.1 plotly_4.10.4
#> [13] igraph_1.6.0 tidyr_1.3.0
#> [15] DT_0.31 limma_3.58.1
#> [17] pheatmap_1.0.12 dplyr_1.1.4
#> [19] reshape2_1.4.4 survminer_0.4.9
#> [21] ggpubr_0.6.0 ggplot2_3.4.4
#> [23] ClassifyR_3.6.2 survival_3.5-7
#> [25] BiocParallel_1.36.0 MultiAssayExperiment_1.28.0
#> [27] SummarizedExperiment_1.32.0 Biobase_2.62.0
#> [29] GenomicRanges_1.54.1 GenomeInfoDb_1.38.5
#> [31] IRanges_2.36.0 MatrixGenerics_1.14.0
#> [33] matrixStats_1.2.0 generics_0.1.3
#> [35] SeuratObject_5.0.1 sp_2.1-2
#> [37] scFeatures_1.3.2 S4Vectors_0.40.2
#> [39] BiocGenerics_0.48.1 BiocStyle_2.30.0
#>
#> loaded via a namespace (and not attached):
#> [1] R.methodsS3_1.8.2 GSEABase_1.64.0
#> [3] progress_1.2.3 EnsDb.Mmusculus.v79_2.99.0
#> [5] goftest_1.2-3 Biostrings_2.70.1
#> [7] HDF5Array_1.30.0 vctrs_0.6.5
#> [9] spatstat.random_3.2-2 digest_0.6.34
#> [11] png_0.1-8 shape_1.4.6
#> [13] ggrepel_0.9.5 deldir_2.0-2
#> [15] parallelly_1.36.0 magick_2.8.2
#> [17] MASS_7.3-60.0.1 httpuv_1.6.13
#> [19] foreach_1.5.2 qvalue_2.34.0
#> [21] withr_3.0.0 ggfun_0.1.3
#> [23] xfun_0.41 ellipsis_0.3.2
#> [25] memoise_2.0.1 proxyC_0.3.4
#> [27] commonmark_1.9.0 gson_0.1.0
#> [29] tidytree_0.4.6 zoo_1.8-12
#> [31] GlobalOptions_0.1.2 gtools_3.9.5
#> [33] SingleCellSignalR_1.14.0 R.oo_1.25.0
#> [35] prettyunits_1.2.0 KEGGREST_1.42.0
#> [37] promises_1.2.1 httr_1.4.7
#> [39] rstatix_0.7.2 restfulr_0.0.15
#> [41] globals_0.16.2 rhdf5filters_1.14.1
#> [43] rhdf5_2.46.1 babelgene_22.9
#> [45] curl_5.2.0 zlibbioc_1.48.0
#> [47] ScaledMatrix_1.10.0 ggraph_2.1.0
#> [49] polyclip_1.10-6 GenomeInfoDbData_1.2.11
#> [51] SparseArray_1.2.3 xtable_1.8-4
#> [53] stringr_1.5.1 evaluate_0.23
#> [55] S4Arrays_1.2.0 BiocFileCache_2.10.1
#> [57] hms_1.1.3 bookdown_0.37
#> [59] irlba_2.3.5.1 colorspace_2.1-0
#> [61] filelock_1.0.3 spatstat.data_3.0-4
#> [63] magrittr_2.0.3 ggtree_3.10.0
#> [65] viridis_0.6.4 later_1.3.2
#> [67] lattice_0.22-5 spatstat.geom_3.2-7
#> [69] future.apply_1.11.1 genefilter_1.84.0
#> [71] shadowtext_0.1.2 XML_3.99-0.16
#> [73] scuttle_1.12.0 cowplot_1.1.2
#> [75] ggupset_0.3.0 pillar_1.9.0
#> [77] nlme_3.1-164 iterators_1.0.14
#> [79] caTools_1.18.2 compiler_4.3.2
#> [81] beachmat_2.18.0 stringi_1.8.3
#> [83] tensor_1.5 GenomicAlignments_1.38.2
#> [85] plyr_1.8.9 crayon_1.5.2
#> [87] abind_1.4-5 BiocIO_1.12.0
#> [89] gridGraphics_0.5-1 ggtext_0.1.2
#> [91] locfit_1.5-9.8 graphlayouts_1.0.2
#> [93] bit_4.0.5 fastmatch_1.1-4
#> [95] codetools_0.2-19 BiocSingular_1.18.0
#> [97] crosstalk_1.2.1 bslib_0.6.1
#> [99] multtest_2.58.0 mime_0.12
#> [101] splines_4.3.2 markdown_1.12
#> [103] circlize_0.4.15 Rcpp_1.0.12
#> [105] dbplyr_2.4.0 sparseMatrixStats_1.14.0
#> [107] HDO.db_0.99.1 gridtext_0.1.5
#> [109] knitr_1.45 blob_1.2.4
#> [111] utf8_1.2.4 fs_1.6.3
#> [113] listenv_0.9.0 DelayedMatrixStats_1.24.0
#> [115] GSVA_1.50.0 ggplotify_0.1.2
#> [117] ggsignif_0.6.4 tibble_3.2.1
#> [119] Matrix_1.6-5 statmod_1.5.0
#> [121] tweenr_2.0.2 pkgconfig_2.0.3
#> [123] tools_4.3.2 cachem_1.0.8
#> [125] RSQLite_2.3.4 viridisLite_0.4.2
#> [127] DBI_1.2.1 fastmap_1.1.1
#> [129] rmarkdown_2.25 scales_1.3.0
#> [131] Rsamtools_2.18.0 broom_1.0.5
#> [133] sass_0.4.8 patchwork_1.2.0
#> [135] BiocManager_1.30.22 dotCall64_1.1-1
#> [137] graph_1.80.0 carData_3.0-5
#> [139] farver_2.1.1 scatterpie_0.2.1
#> [141] tidygraph_1.3.0 yaml_2.3.8
#> [143] rtracklayer_1.62.0 cli_3.6.2
#> [145] purrr_1.0.2 lifecycle_1.0.4
#> [147] bluster_1.12.0 backports_1.4.1
#> [149] annotate_1.80.0 gtable_0.3.4
#> [151] rjson_0.2.21 progressr_0.14.0
#> [153] parallel_4.3.2 ape_5.7-1
#> [155] jsonlite_1.8.8 edgeR_4.0.10
#> [157] bitops_1.0-7 bit64_4.0.5
#> [159] Rtsne_0.17 yulab.utils_0.1.3
#> [161] spatstat.utils_3.0-4 BiocNeighbors_1.20.2
#> [163] ranger_0.16.0 RcppParallel_5.1.7
#> [165] jquerylib_0.1.4 highr_0.10
#> [167] metapod_1.10.1 GOSemSim_2.28.1
#> [169] dqrng_0.3.2 survMisc_0.5.6
#> [171] R.utils_2.12.3 lazyeval_0.2.2
#> [173] shiny_1.8.0 htmltools_0.5.7
#> [175] KMsurv_0.1-5 GO.db_3.18.0
#> [177] rappdirs_0.3.3 glue_1.7.0
#> [179] spam_2.10-0 XVector_0.42.0
#> [181] RCurl_1.98-1.14 treeio_1.26.0
#> [183] scran_1.30.1 gridExtra_2.3
#> [185] AUCell_1.24.0 R6_2.5.1
#> [187] SingleCellExperiment_1.24.0 gplots_3.1.3
#> [189] km.ci_0.5-6 labeling_0.4.3
#> [191] cluster_2.1.6 Rhdf5lib_1.24.1
#> [193] aplot_0.2.2 DelayedArray_0.28.0
#> [195] tidyselect_1.2.0 ProtGenerics_1.34.0
#> [197] ggforce_0.4.1 xml2_1.3.6
#> [199] car_3.1-2 future_1.33.1
#> [201] rsvd_1.0.5 munsell_0.5.0
#> [203] KernSmooth_2.23-22 htmlwidgets_1.6.4
#> [205] fgsea_1.28.0 RColorBrewer_1.1-3
#> [207] biomaRt_2.58.0 rlang_1.1.3
#> [209] spatstat.sparse_3.0-3 spatstat.explore_3.2-5
#> [211] fansi_1.0.6