## ----include=FALSE, cache=FALSE----------------------------------------------- library(CoGAPS) library(BiocParallel) ## ----------------------------------------------------------------------------- packageVersion("CoGAPS") ## ----eval=FALSE--------------------------------------------------------------- # devtools::install_github("FertigLab/CoGAPS") # # # To install via BioConductor: # install.packages("BiocManager") # BiocManager::install("FertigLab/CoGAPS") ## ----------------------------------------------------------------------------- library(CoGAPS) ## ----modsim------------------------------------------------------------------- load('../data/modsimdata.rda') # input to CoGAPS modsimdata ## ----modsim params------------------------------------------------------------ # create new parameters object params <- new("CogapsParams") # view all parameters params # get the value for a specific parameter getParam(params, "nPatterns") # set the value for a specific parameter params <- setParam(params, "nPatterns", 3) getParam(params, "nPatterns") ## ----------------------------------------------------------------------------- # run CoGAPS with specified parameters cogapsresult <- CoGAPS(modsimdata, params, outputFrequency = 10000) cogapsresult ## ----modsim result------------------------------------------------------------ cogapsresult cogapsresult@sampleFactors cogapsresult@featureLoadings # check reference result: load('../data/modsimresult.rda') ## ----load single cell data from zenodo---------------------------------------- # OPTION: download data object from Zenodo options(timeout=50000) # adjust this if you're getting timeout downloading the file url = "https://zenodo.org/record/7709664/files/inputdata.Rds" download.file(url,"inputdata_download.Rds", mode="wb") pdac_data <- readRDS("inputdata_download.Rds") pdac_data ## ----extract counts matrix, eval=FALSE---------------------------------------- # pdac_epi_counts <- as.matrix(pdac_data@assays$originalexp@counts) # norm_pdac_epi_counts <- log1p(pdac_epi_counts) # # head(pdac_epi_counts, n=c(5L, 2L)) # head(norm_pdac_epi_counts, n=c(5L, 2L)) ## ----set params--------------------------------------------------------------- library(CoGAPS) pdac_params <- CogapsParams(nIterations=50000, # 50000 iterations seed=42, # for consistency across stochastic runs nPatterns=8, # each thread will learn 8 patterns sparseOptimization=TRUE, # optimize for sparse data distributed="genome-wide") # parallelize across sets pdac_params ## ----distributed params------------------------------------------------------- pdac_params <- setDistributedParams(pdac_params, nSets=7) pdac_params ## ----run CoGAPS on single-cell data, eval=FALSE------------------------------- # startTime <- Sys.time() # # pdac_epi_result <- CoGAPS(pdac_epi_counts, pdac_params) # endTime <- Sys.time() # # saveRDS(pdac_epi_result, "../data/pdac_epi_cogaps_result.Rds") # # # To save as a .csv file, use the following line: # toCSV(pdac_epi_result, "../data") ## ----load precomputed from zenodo--------------------------------------------- # OPTION: download precomputed CoGAPS result object from Zenodo # Bioc vignette uses this options(timeout=50000) # adjust this if you're getting timeout downloading the file url = "https://zenodo.org/record/7709664/files/cogapsresult.Rds?download=1" download.file(url,"cogapsresult_download.Rds", mode="wb") cogapsresult <- readRDS("cogapsresult_download.Rds") ## ----load saved, eval=FALSE--------------------------------------------------- # library(CoGAPS) # cogapsresult <- readRDS("../data/pdac_epi_cogaps_result.Rds") ## ----pattern assay, eval=FALSE------------------------------------------------ # library(Seurat) # # make sure pattern matrix is in same order as the input data # patterns_in_order <-t(cogapsresult@sampleFactors[colnames(pdac_data),]) # # # add CoGAPS patterns as an assay # pdac_data[["CoGAPS"]] <- CreateAssayObject(counts = patterns_in_order) ## ----pattern UMAP, eval=FALSE------------------------------------------------- # DefaultAssay(pdac_data) <- "CoGAPS" # pattern_names = rownames(pdac_data@assays$CoGAPS) # # library(viridis) # color_palette <- viridis(n=10) # # FeaturePlot(pdac_data, pattern_names, cols=color_palette, reduction = "umap") & NoLegend() ## ----pattern markers---------------------------------------------------------- pm <- patternMarkers(cogapsresult, threshold="cut") ## ----------------------------------------------------------------------------- hallmarks <- getPatternHallmarks(cogapsresult) ## ----------------------------------------------------------------------------- pl_pattern7 <- plotPatternHallmarks(cogapsresult, hallmarks, whichpattern = 7) pl_pattern7 ## ----MANOVA variables--------------------------------------------------------- # create dataframe of interested variables interestedVariables <- cbind(pdac_data@meta.data[["nCount_RNA"]], pdac_data@meta.data[["nFeature_RNA"]]) # call cogaps manova wrapper manova_results <- MANOVA(interestedVariables, cogapsresult)