Merge pull request #825 from AlexsLemonade/auto_render_live

GHA: Automated live (rendered) versions of the notebooks

RNA-seq/02-gastric_cancer_tximeta.nb.html

@@ -3264,8 +3264,8 @@ <h2>Summarize to gene</h2>
 <pre class="r"><code># Summarize to the gene level
 gene_summarized &lt;- summarizeToGene(txi_data)</code></pre>
 <!-- rnb-source-end -->
-<!-- rnb-message-begin eyJkYXRhIjoibG9hZGluZyBleGlzdGluZyBFbnNEYiBjcmVhdGVkOiAyMDI0LTA4LTA4IDE2OjEwOjI3XG4ifQ== -->
-<pre><code>loading existing EnsDb created: 2024-08-08 16:10:27</code></pre>
+<!-- rnb-message-begin eyJkYXRhIjoibG9hZGluZyBleGlzdGluZyBFbnNEYiBjcmVhdGVkOiAyMDI0LTEyLTA0IDIxOjU5OjQ0XG4ifQ== -->
+<pre><code>loading existing EnsDb created: 2024-12-04 21:59:44</code></pre>
 <!-- rnb-message-end -->
 <!-- rnb-message-begin eyJkYXRhIjoib2J0YWluaW5nIHRyYW5zY3JpcHQtdG8tZ2VuZSBtYXBwaW5nIGZyb20gZGF0YWJhc2VcbiJ9 -->
 <pre><code>obtaining transcript-to-gene mapping from database</code></pre>

scRNA-seq-advanced/01-read_filter_normalize_scRNA-live.Rmd

@@ -90,9 +90,7 @@ Finally, we will set up the our output directory, creating it if it does not yet
 normalized_dir <- file.path(data_dir, "normalized")
 
 # create the directory if it does not exist
-if (!dir.exists(normalized_dir)) {
-  dir.create(normalized_dir, recursive = TRUE)
-}
+fs::dir_create(normalized_dir)
 
 # output RDS file for normalized data
 output_sce_file <- file.path(normalized_dir,

scRNA-seq-advanced/02-dataset_integration-live.Rmd

@@ -525,13 +525,12 @@ Now, let's see how this new `merged_UMAP` looks compared to the `UMAP` calculate
 # UMAPs scaled together when calculated from the merged SCE
 scater::plotReducedDim(merged_sce,
                        dimred = "merged_UMAP",
-                       colour_by = "sample",
+                       color_by = "sample",
                        # Some styling to help us see the points:
                        point_size = 0.5,
                        point_alpha = 0.2) +
-  # Modify the legend key so its points are larger and easier to see
-  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
-  # Add a plot title
+  scale_color_brewer(palette = "Dark2", name = "sample") +
+  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
   ggtitle("UMAP calculated on merged_sce")
 ```
 
@@ -596,13 +595,12 @@ scater::plotReducedDim(merged_sce,
                        # plot the fastMNN coordinates
                        dimred = "fastmnn_UMAP",
                        # color by sample
-                       colour_by = "sample",
+                       color_by = "sample",
                        # Some styling to help us see the points:
                        point_size = 0.5,
                        point_alpha = 0.2) +
-  # Modify legend so they key is larger and easier to see
-  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
-  # add plot title
+  scale_color_brewer(palette = "Dark2", name = "sample") +
+  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
   ggtitle("UMAP after integration with fastMNN")
 ```
 
@@ -642,10 +640,12 @@ Let's re-plot this UMAP to highlight cell types:
 scater::plotReducedDim(merged_sce,
                        dimred = "fastmnn_UMAP",
                        # color by broad celltypes
-                       colour_by = "celltype_broad",
+                       color_by = "celltype_broad",
                        point_size = 0.5,
                        point_alpha = 0.2) +
-  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
+  # include argument to specify color of NA values
+  scale_color_brewer(palette = "Dark2", name = "Broad celltype", na.value = "grey80") +
+  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
   ggtitle("UMAP after integration with fastMNN")
 ```
 
@@ -658,12 +658,13 @@ One way we can see all the points a bit better is to facet the plot by sample, u
 ```{r plot fastmnn umap celltypes faceted}
 scater::plotReducedDim(merged_sce,
                        dimred = "fastmnn_UMAP",
-                       colour_by = "celltype_broad",
+                       color_by = "celltype_broad",
                        point_size = 0.5,
                        point_alpha = 0.2,
                        # Allow for faceting by a variable using `other_fields`:
                        other_fields = "sample") +
-  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
+  scale_color_brewer(palette = "Dark2", name = "Broad celltype", na.value = "grey80") +
+  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
   ggtitle("UMAP after integration with fastMNN") +
   # Facet by sample
   facet_wrap(vars(sample)) +
@@ -743,11 +744,12 @@ Let's see how the `harmony` UMAP, colored by sample, looks compared to the `fast
 ```{r plot harmony umap batches}
 scater::plotReducedDim(merged_sce,
                        dimred = "harmony_UMAP",
-                       colour_by = "sample",
+                       color_by = "sample",
                        point_size = 0.5,
                        point_alpha = 0.2) +
-  ggtitle("UMAP after integration with harmony") +
-  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1)))
+  scale_color_brewer(palette = "Dark2", name = "sample") +
+  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
+  ggtitle("UMAP after integration with harmony")
 ```
 
 How do you think this `harmony` UMAP compares to that from `fastMNN` integration?
@@ -757,19 +759,19 @@ Let's see how this UMAP looks colored by cell type, and faceted for visibility:
 ```{r plot harmony umap celltypes}
 scater::plotReducedDim(merged_sce,
                        dimred = "harmony_UMAP",
-                       colour_by = "celltype_broad",
+                       color_by = "celltype_broad",
                        point_size = 0.5,
                        point_alpha = 0.2,
                        # Specify variable for faceting
                        other_fields = "sample") +
+  scale_color_brewer(palette = "Dark2", name = "Broad celltype", na.value = "grey80") +
+  guides(color = guide_legend(override.aes = list(size = 3))) +
   ggtitle("UMAP after integration with harmony") +
-  guides(colour = guide_legend(override.aes = list(size = 3))) +
   facet_wrap(vars(sample))
 ```
 
-With the faceted view, we can now see that the small amounts of `Tumor_Myocyte` from other samples are "pulled" towards those cells in `SCPCL000481`.
-
-What other patterns do you see that are similar or different from the `fastMNN` UMAP?
+What do you now notice in this faceted view that wasn't clear previously?
+Are there other patterns you see that are similar or different from the `fastMNN` UMAP?
 How do you think `fastMNN` vs. `harmony` performed in integrating these samples?
 
 ### Export

scRNA-seq-advanced/03-differential_expression-live.Rmd

@@ -26,7 +26,7 @@ In this notebook, we will work with multiple samples to identify differentially
 ![Single-cell roadmap: Differential expression](diagrams/roadmap_differential_expression.png)
 
 We will continue working with samples from the [`SCPCP000005` project](https://scpca.alexslemonade.org/projects/SCPCP000005), an investigation of pediatric solid tumors led by the Dyer and Chen labs at St. Jude Children's Research Hospital.
-This particular dataset contains 10 different samples that have been integrated using `fastMNN`, following the same procedure we outlined in `03-dataset_integration.Rmd`.
+This particular dataset contains 10 different samples that have been integrated using `fastMNN`, following the same procedure we outlined in `02-dataset_integration.Rmd`.
 These 10 samples represent two different types of rhabdomyosarcoma (RMS): embryonal rhabdomyosarcoma (ERMS) and alveolar rhabdomyosarcoma (ARMS).
 These two subtypes are distinguished by the presence of the `PAX3/PAX7-FOXO1` fusion gene, which is present only in ARMS patients.
 Additionally, cells found in ARMS tumors tend to have an increased mutational burden with cells in a more differentiated state compared to ERMS tumor cells ([Shern _et al._ 2014](https://doi.org/10.1158/2159-8290.CD-13-0639); [Stewart _et al._ 2018](https://doi.org/10.1016/j.ccell.2018.07.012)).
@@ -79,9 +79,7 @@ sample_metadata_file <- file.path(data_dir,
 
 # directory to store output
 deseq_dir <- file.path("analysis", "rms", "deseq")
-if(!dir.exists(deseq_dir)){
-  dir.create(deseq_dir, recursive = TRUE)
-}
+fs::dir_create(deseq_dir)
 
 # results file to output from DE analysis
 deseq_output_file <- file.path(deseq_dir, 
@@ -223,7 +221,7 @@ The samples which could be further classified have a mix of `Tumor_Mesoderm`, `T
 scater::plotReducedDim(integrated_sce,
                        dimred = "fastmnn_UMAP",
                        # color each point by cell type
-                       colour_by = "celltype_broad",
+                       color_by = "celltype_broad",
                        point_size= 0.5, 
                        point_alpha = 0.4)
 ```
@@ -267,6 +265,7 @@ ggplot(tumor_cells_df, aes(x = sample, fill = celltype_broad)) +
       y = "Proportion of cells", 
       fill = "Cell type"
     ) +
+  scale_fill_brewer(palette = "Dark2") +
   theme_bw() +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
   # facet by diagnosis group 
@@ -694,7 +693,7 @@ scater::plotExpression(tumor_sce,
                        # a vector of genes to plot
                        features = genes_to_plot, 
                        x = "diagnosis_group", 
-                       colour_by = "diagnosis_group",
+                       color_by = "diagnosis_group",
                        other_fields = "celltype_broad",
                        point_size = 0.1) +
   # each celltype is its own column
@@ -703,7 +702,7 @@ scater::plotExpression(tumor_sce,
              rows = vars(Feature)) + 
   # change the font size of the facet labels
   theme(strip.text = element_text(size = 7)) + 
-  guides(colour = guide_legend(
+  guides(color = guide_legend(
     title = "Subtype", # update the legend title
     # change the size of the legend colors
     override.aes = list(size = 3, alpha = 1))