diff --git a/scRNA-seq/05-clustering_markers_scRNA.Rmd b/scRNA-seq/05-clustering_markers_scRNA.Rmd index 3a3ba0241..0739552fd 100644 --- a/scRNA-seq/05-clustering_markers_scRNA.Rmd +++ b/scRNA-seq/05-clustering_markers_scRNA.Rmd @@ -36,7 +36,11 @@ set.seed(12345) data_dir <- file.path("data", "hodgkins") # normalized data file -normalized_rds <- file.path(data_dir, "normalized", "normalized_hodgkins_sce.rds") +normalized_rds <- file.path( + data_dir, + "normalized", + "normalized_hodgkins_sce.rds" +) # Output directory for markers marker_dir <- file.path("analysis", "hodgkins", "markers") @@ -170,10 +174,12 @@ Now we can plot the results of our graph-based clustering. This time we will also use the `text_by` argument to include the cluster ids directly on the plot. ```{r plot_nnclust, live = TRUE} -plotReducedDim(hodgkins_sce, - "UMAP", - color_by = "nncluster", - text_by = "nncluster") +plotReducedDim( + hodgkins_sce, + "UMAP", + color_by = "nncluster", + text_by = "nncluster" +) ``` - How do these results compare to the k-means clustering result? @@ -212,9 +218,11 @@ For now, we will use the `findMarkers()` function to rank the genes in each clus ```{r find_markers, live = TRUE} # use `findMarkers()` to calculate how well each gene # differentiates each cluster from *all* other clusters -markers <- scran::findMarkers(hodgkins_sce, - groups = hodgkins_sce$nncluster, - pval.type = "all") +markers <- scran::findMarkers( + hodgkins_sce, + groups = hodgkins_sce$nncluster, + pval.type = "all" +) ``` Next we can look at one of those tables. @@ -242,14 +250,14 @@ The results are returned as a new list. purrr::map( as.list(markers[1:3]), # select the first 3 clusters and convert to a 'regular' list for purrr head # the function to apply (note no parenthesis) - ) +) ``` This returns a list of data frames, which isn't quite what we want. There is no built-in function that will give us just the first few _row names_, so we will have to define one. -As of version 4.1, R introduced a new approach to defining _anonymous functions_ - that is, functions you can quickly define "on-the-fly" without formally assigning them to a function name. -They are handy when you need to do a very short task that requires a function, but it isn't really a function you need beyond this context. +We'll do this with an "anonymous" function (sometimes called a lambda expression). +Anonymous functions are handy when you need to do a very short task that requires a function, but it isn't really a function you need beyond this context. This new anonymous syntax looks like this: `\(x)...` (or for slightly longer code, use curly braces as in `\(x) {...}`). This defines a function that takes one argument, `x`, with `...` indicating where you would put the expression to calculate. @@ -261,7 +269,7 @@ purrr::map( # convert markers to a 'regular' list for purrr as.list(markers), # our custom function: - \(x) head( rownames(x) ) + \(x) head(rownames(x)) ) ``` @@ -288,12 +296,13 @@ In addition to writing the tables out, we are saving the data frames we created ```{r write_tables} marker_df_list <- purrr::imap( as.list(markers), # convert markers to a 'regular' list for purrr - # purrr function: x is the list element, y is the element name (number here) + # anonymous function: x is the list element, y is the element name (number here) \(table, id) { as.data.frame(table) |> # first convert to a data frame tibble::rownames_to_column("gene") |> # make genes a column dplyr::arrange(FDR) |> # sort to be sure small FDR genes are first - readr::write_tsv( # write each data frame to a file + # write each data frame to a file + readr::write_tsv( file.path( marker_dir, # construct the output path sprintf("cluster%02d_markers.tsv", as.integer(id)) # format cluster numbers in file names with leading zeros @@ -332,8 +341,7 @@ gene_id <- gene_info$ID[rank] symbol <- gene_info$Symbol[rank] # Plot UMAP results colored by expression -plotReducedDim(hodgkins_sce, "UMAP", - color_by = gene_id) + +plotReducedDim(hodgkins_sce, "UMAP", color_by = gene_id) + # label the guide with the gene symbol guides(color = guide_colorbar(title = symbol)) ``` diff --git a/scRNA-seq/06-celltype_annotation.Rmd b/scRNA-seq/06-celltype_annotation.Rmd index 83e6e0fec..76b907c48 100644 --- a/scRNA-seq/06-celltype_annotation.Rmd +++ b/scRNA-seq/06-celltype_annotation.Rmd @@ -152,7 +152,7 @@ While we have the ADT data, it is _not_ being used for this stage of the analysi ```{r cluster cells, live=TRUE} # perform clustering nn_clusters <- scran::clusterCells( - sce, + sce, use.dimred = "PCA", # PCA input BLUSPARAM = bluster::NNGraphParam() # graph clustering & parameters ) @@ -470,9 +470,11 @@ Still annoyingly, however, when we change the palette, the legend title defaults ```{r plot celltype umap palette} scater::plotUMAP(sce, color_by = "celltype_main") + - scale_color_brewer(name = "Cell type", # legend title - palette = "Dark2", # color palette - na.value = "gray80") # use light gray for NA values + scale_color_brewer( + name = "Cell type", # legend title + palette = "Dark2", # color palette + na.value = "gray80" + ) # use light gray for NA values ``` We seem to have a pretty good set of cell type assignments, with most falling into groupings consistent with what we see in the UMAP plot. @@ -519,7 +521,7 @@ singler_finemodel <- SingleR::trainSingleR( monaco_ref, # reference dataset labels = monaco_ref$label.fine, # labels for training dataset # consider only genes in the sce object - test.genes = rownames(sce), + test.genes = rownames(sce), # use DE to select genes (default) genes = "de", # parallel processing @@ -569,28 +571,32 @@ After that, we will use `fct_relevel` to put the remaining factor levels in the collapsed_labels <- singler_result_fine$pruned.labels |> forcats::fct_collapse( "Monocytes" = c( - "Classical monocytes", - "Intermediate monocytes", - "Non classical monocytes"), + "Classical monocytes", + "Intermediate monocytes", + "Non classical monocytes" + ), "Dendritic cells" = c( - "Myeloid dendritic cells", - "Plasmacytoid dendritic cells"), + "Myeloid dendritic cells", + "Plasmacytoid dendritic cells" + ), "T cells" = c( - "MAIT cells", - "Non-Vd2 gd T cells", - "Vd2 gd T cells"), + "MAIT cells", + "Non-Vd2 gd T cells", + "Vd2 gd T cells" + ), "Helper T cells" = c( - "Th1 cells", - "Th1/Th17 cells", - "Th17 cells", - "Th2 cells", - "Follicular helper T cells"), + "Th1 cells", + "Th1/Th17 cells", + "Th17 cells", + "Th2 cells", + "Follicular helper T cells" + ), "B cells" = c( - "Naive B cells", - "Switched memory B cells", - "Non-switched memory B cells", - "Exhausted B cells", - "Plasmablasts" + "Naive B cells", + "Switched memory B cells", + "Non-switched memory B cells", + "Exhausted B cells", + "Plasmablasts" ) ) |> # order for plotting @@ -646,7 +652,7 @@ Since the total number of cells differs across clusters, we'd like to convert th We'll do this by going through the table column by column and dividing each value by the sum for that cluster. This will give us normalized values where the values in each column now sum to 1. To do that, we will use the `apply` function, which allows us to operate on a matrix row by row or column by column, applying a function to each "slice". -Since the function we want to apply is very short, we will use R's new (as of version 4.1) anonymous function shorthand: +Since the function we want to apply is very short, we will use R's anonymous function shorthand: `\(x) ...` can be used to define a function that that takes as input values `x` (where the `...` is where you would put the expression to calculate). Here we will apply the expression `x/sum(x)`, which will divide each element of a vector `x` by the sum of its values. @@ -655,7 +661,7 @@ Here we will apply the expression `x/sum(x)`, which will divide each element of type_cluster_tab <- apply( type_cluster_tab, 2, # apply function to columns - \(x) x/sum(x) # function to apply + \(x) x / sum(x) # function to apply ) # print the normalized values type_cluster_tab[1:5, 1:5] diff --git a/scripts/link-data.sh b/scripts/link-data.sh index 3e7ada669..79f430071 100644 --- a/scripts/link-data.sh +++ b/scripts/link-data.sh @@ -76,7 +76,7 @@ link_locs=( spatial/data/crc-hd/crc_single_cell_normalized_sce.rds spatial/data/crc-hd/single_cell spatial/data/crc-hd/visium_hd - spatial/data/crc-v1/crc_normalized_spe.rds + spatial/data/crc-v1/normalized/crc_normalized_spe.rds spatial/data/osteo/GSM8478586/normalized/osteo_normalized_spe.rds spatial/data/ovarian-carcinoma/spaceranger spatial/data/reference diff --git a/spatial/02-spatial_clustering.Rmd b/spatial/02-spatial_clustering.Rmd index f871d4f19..57bc66b6f 100644 --- a/spatial/02-spatial_clustering.Rmd +++ b/spatial/02-spatial_clustering.Rmd @@ -189,7 +189,7 @@ gene_variance <- scran::modelGeneVar(spe) # calculate the top 2000 HVGs hvg_vector <- scran::getTopHVGs(gene_variance, n = 2000) -# calculate PCE using the top HVGs, saving back to the SPE object +# calculate PCA using the top HVGs, saving back to the SPE object spe <- scater::runPCA(spe, subset_row = hvg_vector) ``` diff --git a/spatial/03-deconvolution.Rmd b/spatial/03-deconvolution.Rmd index 2e0bec15f..7ce299b14 100644 --- a/spatial/03-deconvolution.Rmd +++ b/spatial/03-deconvolution.Rmd @@ -445,8 +445,8 @@ We'll save the normalized values to a new assay called `proportions`. assay(rctd_full, "proportions") <- assay(rctd_full, "weights") |> # convert to regular matrix for normalization as.matrix() |> - # new in R 4.1, a function to normalize! - # use 2 for columns + # normalize columns so they sum to 1 + # use 2 for columns (1 would do this for rows) proportions(margin = 2) |> # convert back to sparse matrix as("CsparseMatrix")