--- title: "PBMC 3k analysis" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE ) ``` Let's see how Rmd files are handled on Lamin's end by performing an example Seurat analysis. First we connect to the instance we want to use. This must be done ***BEFORE*** importing `lamindb`. If you don't explicitly connect, the current default instance will be used. ```{r connect} lc <- laminr::import_module("lamin_cli", silent = TRUE) lc$connect("laminlabs/training") ``` Now we can import the `lamindb` module to run commands. ```{r import} library(laminr) ln <- import_module("lamindb", silent = TRUE) ``` From here, we're going to track every move of this session. ```{r track} # Set up a project proj <- ln$Project(name = "Basic Seurat analysis")$save() # Start a tracked run ln$track(project = proj) ``` Great. Now let's get the PBMC 3k dataset in here. ```{r libraries} library(tidyverse) library(Seurat) ``` Let's get this going. ```{r load} counts <- read.table( file = system.file('extdata', 'pbmc_raw.txt', package = 'Seurat'), as.is = TRUE ) cells <- CreateSeuratObject(counts) ``` And now we're going to run it through the usual set of things, using Seurat's best practices. ```{r normalise} cells <- NormalizeData( cells, normalization.method = "LogNormalize", scale.factor = 10000 ) ``` ```{r variable-features} cells <- FindVariableFeatures(cells, selection.method = "vst", nfeatures = 2000) ``` Note that this step makes the Seurat object much more bulky. ```{r scale} all.genes <- rownames(cells) cells <- ScaleData(cells, features = all.genes) ``` From here, we do our PCA. ```{r pca} cells <- RunPCA(cells, features = VariableFeatures(object = cells)) ``` Which sets the stage for the clustering and non-linear dimensionality reduction steps below. We can get a feel for the intrinsic dimensionality of the data below. ```{r elbow} ElbowPlot(cells) ``` We then do clustering, which is Louvain clustering of a SNN graph. ```{r cluster} cells <- FindNeighbors(cells, dims = 1:10) cells <- FindClusters(cells, resolution = 0.5) ``` And from here, we run a UMAP for visualization. ```{r umap} cells <- RunUMAP(cells, dims = 1:10) ``` Now let's look at it. ```{r plot-umap} # note that you can set `label = TRUE` or use the LabelClusters function to # help label individual clusters DimPlot(cells, reduction = "umap") ``` More or less what was expected. And let's color by some features. ```{r plot-features} FeaturePlot(cells, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A")) ``` Cool. It's really ugly when you visualize it in the Rmd interface, but when you knit this thing, it's ok. Ok, and from here, we're going to set this as a Lamin artifact and save it out accordingly. We tried to save it as a h5ad object, but that didn't work. So we're just going to go ahead and output as a rds, which is standard when you're handing Seurat. There are ways to convert it to h5ad that are beyond the scope of this test here. And now we save. ```{r save-rds} saveRDS(cells, "pbmc3k_processed.rds") ``` And now we turn the saved file into a Lamin artifact. ```{r save-artifact} # Use this rather than lamin_save ln$Artifact("pbmc3k_processed.rds", key = "pbmc3k/pbmc3k_processed.rds")$save() ``` And now we end the session. ```{r finish} # optional: close out the tracked run so it shows up neatly in the hub ln$finish() ```