Live and rendered notebooks

intro-to-R-tidyverse/03-intro_to_tidyverse-live.Rmd

@@ -444,7 +444,7 @@ stats_df |>
 #### Save to TSV files
 
 Let's write some of the data frames we created to a file.
-To do this, we can use the `readr` library of `_write()` functions.
+To do this, we can use the `readr` library of `write_()` functions.
 The first argument of `write_tsv()` is the data we want to write, and the second argument is a character string that describes the path to the new file we would like to create.
 Remember that we created a `results` directory to put our output in, but if we want to save our data to a directory other than our working directory, we need to specify this.
 This is what we will use the `file.path()` function for.

scRNA-seq/03-normalizing_scRNA-live.Rmd

@@ -146,6 +146,10 @@ One way to determine whether our normalization yields biologically relevant resu
 Because plotting expression for thousands genes together isn't practical, we will reduce the dimensions of our data using Principal Components Analysis (PCA).
 
 We will also make the same plot with our *unnormalized* data, to visualize the effect of normalization on our sample.
+We'll do this comparison twice:
+
+- Once coloring the points by their total UMI count
+- Once coloring the points based on their cell labels
 
 Before plotting the unnormalized data, we will log transform the raw counts to make their scaling more comparable to the normalized data.
 To do this we will use the `log1p()` function, which is specifically designed for the case where we want to add 1 to all of our values before taking the log, as we do here.
@@ -166,27 +170,29 @@ log_pca <- counts(bladder_sce) |> # get the raw counts
 Note that we are using `scater::calculatePCA()` two different ways here: once on the full `bladder_sce` object, and once on just the `counts` matrix.
 When we use `calculatePCA()` on the object, it automatically uses the log normalized matrix from inside the object.
 
-Next we will arrange the PCA scores for plotting, adding a column with the cell-type data so we can color each point of the plot.
+Next we will arrange the PCA scores for plotting, adding a column for each of the total UMI counts and the cell type labels so we can color each point of the plot.
 
 ```{r pca_df}
 # Set up the PCA scores for plotting
 norm_pca_scores <- data.frame(norm_pca,
                               geo_accession = rownames(norm_pca),
+                              total_umi = bladder_sce$sum,
                               cell_type = bladder_sce$cell_ontology_class)
 log_pca_scores <- data.frame(log_pca,
                              geo_accession = rownames(log_pca),
+                             total_umi = bladder_sce$sum,
                              cell_type = bladder_sce$cell_ontology_class)
 ```
 
-Now we will plot the unnormalized PCA scores with their cell labels:
+First, we will plot the unnormalized PCA scores with their total UMI counts:
 
 ```{r pca_plot}
 # Now plot counts pca
-ggplot(log_pca_scores, aes(x = PC1, y = PC2, color = cell_type)) +
+ggplot(log_pca_scores, aes(x = PC1, y = PC2, color = total_umi)) +
   geom_point() +
   labs(title = "Log counts (unnormalized) PCA scores",
-       color = "Cell Type") +
-  scale_color_brewer(palette = "Dark2", na.value = "grey70") + # add a visually distinct color palette
+       color = "Total UMI count")  +
+  scale_color_viridis_c() +
   theme_bw()
 ```
 
@@ -200,7 +206,31 @@ Let's plot the `norm_pca_scores` data:
 
 ```
 
-Do you see an effect from the normalization in the comparison between these plots?
+Do you see an effect from the normalization when comparing these plots?
+
+
+
+Now, let's plot these two sets of PCA scores again, but colored by cell type.
+Do you see an effect from the normalization when comparing these plots?
+
+```{r celltype_pca_plots}
+# First, plot the normalized pca
+ggplot(norm_pca_scores, aes(x = PC1, y = PC2, color = cell_type)) +
+  geom_point() +
+  labs(title = "Normalized log counts PCA scores",
+       color = "Cell Type") +
+  scale_color_brewer(palette = "Dark2", na.value = "grey70") + # add a visually distinct color palette
+  theme_bw()
+
+# Next, plot log count pca
+ggplot(log_pca_scores, aes(x = PC1, y = PC2, color = cell_type)) +
+  geom_point() +
+  labs(title = "Log counts (unnormalized) PCA scores",
+       color = "Cell Type") +
+  scale_color_brewer(palette = "Dark2", na.value = "grey70") + # add a visually distinct color palette
+  theme_bw()
+```
+
 
 
 ## Save the normalized data to tsv file

scRNA-seq/04-dimension_reduction_scRNA-live.Rmd

@@ -206,13 +206,53 @@ dim(filtered_sce)
 Now we will perform the same normalization steps we did in a previous dataset, using `scran::computeSumFactors()` and `scater::logNormCounts()`.
 You might recall that there is a bit of randomness in some of these calculations, so we should be sure to have used `set.seed()` earlier in the notebook for reproducibility.
 
-```{r normalize}
+```{r sumfactors}
 # Cluster similar cells
 qclust <- scran::quickCluster(filtered_sce)
 
 # Compute sum factors for each cell cluster grouping.
+filtered_sce <- scran::computeSumFactors(filtered_sce, clusters = qclust, positive = FALSE)
+```
+
+It turns out in this case we end up with some negative size factors.
+This is usually an indication that our filtering was not stringent enough, and there remain a number of cells or genes with nearly zero counts.
+This probably happened when we removed the infrequently-expressed genes; cells which had high counts from those particular genes (and few others) could have had their total counts dramatically reduced.
+
+To account for this, we will recalculate the per-cell stats and filter out low counts.
+Unfortunately, to do this, we need to first remove the previously calculated statistics, which we will do by setting them to `NULL`.
+
+```{r reQC}
+# remove previous calculations
+filtered_sce$sum <- NULL
+filtered_sce$detected <- NULL
+filtered_sce$total <- NULL
+filtered_sce$subsets_mito_sum <- NULL
+filtered_sce$subsets_mito_detected <- NULL
+filtered_sce$subsets_mito_sum <- NULL
+
+# recalculate cell stats
+filtered_sce <- scater::addPerCellQC(filtered_sce, subsets = list(mito = mito_genes))
+
+# print the number of cells with fewer than 500 UMIs
+sum(filtered_sce$sum < 500)
+```
+
+Now we can filter again.
+In this case, we will keep cells with at least 500 UMIs after removing the lowly expressed genes.
+Then we will redo the size factor calculation, hopefully with no more warnings.
+
+
+```{r refilter}
+filtered_sce <- filtered_sce[, filtered_sce$sum >= 500]
+
+qclust <- scran::quickCluster(filtered_sce)
+
 filtered_sce <- scran::computeSumFactors(filtered_sce, clusters = qclust)
+```
 
+Looks good! Now we'll do the normalization.
+
+```{r normalize}
 # Normalize and log transform.
 normalized_sce <- scater::logNormCounts(filtered_sce)
 ```