Merge pull request #120 from NICHD-BSPC/sally-edits-may2025

Wk 7 Lesson 04: Update plots, add catch-up script for plotting

Author: esallychang

img/maplot_solution.png

@@ -306,31 +306,32 @@ ggplot(res_OE_df_plotting, aes(x = log2FoldChange, y = -log10(padj))) +
 
 ### Selecting Your Own Gene
 
-What if you want to select your own gene of interest? One way is creating a new data frame where you have added a value in the `genelabels` column just for that one symbol.
+What if you want to select your own gene of interest? One way is creating a new data frame where you have added a value in the `genelabels` column just for that one symbol. In this example, I am using a gene symbol called `HOX3A`.
 
-``` r
-res_tableOE_df[res_tableOE_tb$symbol == "Mygene", "genelabels"]
+Preparing for this will take a few steps in Base R:
+
+-   Add an additional empty column called `my_genelabels` to our current plotting dataframe, to put those gene names we want to use to label the plot
 
-res_tableOE_mygene <- res_tableOE_df %>% 
-  dplyr::mutate(threshold_OE = padj < 0.1 & abs(log2FoldChange) >= 0.58)
+-   Assign value just to the one relevant empty `my_genelabel` slot! You can use a strategy like this in other data frames and with other genes - or whole custom lists of genes.
 
-res_tableOE_mygene <- res_tableOE_mygene %>% dplyr::mutate(genelabels = "")
+``` r
+# Create a new empty column
+res_OE_df_plotting$my_genelabels <- ""
 
-# Assign value just to that one empty genelabel slot! You can use a strategy like this in other data frames and with other genes - or whole custom lists of genes. 
-res_tableOE_mygene[res_tableOE_mygene$symbol == "HOXA3", "genelabels"] = "HOXA3"
+# Assign value just to that one empty genelabel slot where symbol == HOXA3
+res_OE_df_plotting[res_OE_df_plotting$symbol == "HOXA3", "genelabels"] = "HOXA3"
 ```
 
 **Now let's check this out in a volcano plot:**
 
 ``` r
-ggplot(res_tableOE_mygene, aes(x = log2FoldChange, y = -log10(padj))) +
-    geom_point(aes(colour = threshold_OE)) +
-    geom_text_repel(aes(label = genelabels)) +
+ggplot(res_OE_df_plotting, aes(x = log2FoldChange, y = -log10(padj))) +
+    geom_point(aes(colour = significant_OE)) +
+    geom_text_repel(aes(label = my_genelabels)) +
     ggtitle("Mov10 overexpression") +
     xlab("log2 fold change") + 
     ylab("-log10 adjusted p-value") +
-    theme(legend.position = "none",
-          plot.title = element_text(size = rel(1.5), hjust = 0.5),
+    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
           axis.title = element_text(size = rel(1.25))) 
 ```
 
@@ -339,15 +340,14 @@ But where is our gene?? Let's actually increase that overlaps parameter...
 ``` r
 options(ggrepel.max.overlaps = Inf)
 
-ggplot(res_tableOE_mygene, aes(x = log2FoldChange, y = -log10(padj))) +
-    geom_point(aes(colour = threshold_OE)) +
-    geom_text_repel(aes(label = genelabels)) +
+ggplot(res_OE_df_plotting, aes(x = log2FoldChange, y = -log10(padj))) +
+    geom_point(aes(colour = significant_OE)) +
+    geom_text_repel(aes(label = my_genelabels)) +
     ggtitle("Mov10 overexpression") +
     xlab("log2 fold change") + 
     ylab("-log10 adjusted p-value") +
-    theme(legend.position = "none",
-          plot.title = element_text(size = rel(1.5), hjust = 0.5),
-          axis.title = element_text(size = rel(1.25))) 
+    theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
+          axis.title = element_text(size = rel(1.25)))  
 ```
 
 <img src="../img/hox3a_labeled_volcano.png" width="500"/>
@@ -360,19 +360,35 @@ Switching from a volcano plot to an MA plot is pretty easy! Let's fill in this s
 
 ``` r
 ggplot(res_tableOE_plotting, aes(x = , y = )) +
-  geom_point(aes(colour = threshold_OE)) +
+  geom_point(aes(colour = significant_OE)) +
   geom_text_repel(aes(label = genelabels)) +
-  ggtitle("Mov10 overexpression") +
+  ggtitle("MA Plot:Mov10 overexpression") +
   xlab("") + 
   ylab("") +
-  theme(legend.position = "none",
-        plot.title = element_text(size = rel(1.5), hjust = 0.5),
+  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
         axis.title = element_text(size = rel(1.25)))
 ```
 
 If we do this correctly, we should have a something like the following - depending on what you have saved as your `genelabels` variable and how big your plot window is:
 
-<img src="../img/mov10_ggplot_maplot.png" width="450"/>
+<img src="../img/maplot_solution.png" alt="MA plot with top 10 DE genes labeled for overexpression comparison" width="600"/>
+
+<details>
+
+<summary>Code Solution:</summary>
+
+``` r
+ggplot(res_OE_df_plotting, aes(x = log(baseMean) , y = log2FoldChange)) +
+  geom_point(aes(colour = significant_OE)) +
+  geom_text_repel(aes(label = genelabels)) +
+  ggtitle("MA Plot for Mov10 overexpression contrast") +
+  xlab("log of meant counts") + 
+  ylab("log2FoldChange") +
+  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
+        axis.title = element_text(size = rel(1.25)))
+```
+
+</details>
 
 ## **ASSIGNMENT OPTION 2**:
 
@@ -382,12 +398,12 @@ Take the above steps (starting with finding the Ensembl Gene ID) for a gene symb
 
 ------------------------------------------------------------------------
 
-> ### An R package for visualization of DGE results
->
-> The Bioconductor package [`DEGreport`](https://bioconductor.org/packages/release/bioc/html/DEGreport.html) can use the DESeq2 results output to make the top20 genes and the volcano plots generated above by writing much fewer lines of code. The caveat of these functions is you lose the ability to customize plots as we have demonstrated above.
->
-> If you are interested, the example code below shows how you can use DEGreport to create similar plots. **Note that this is example code, do not run.**
->
+### An R package for visualization of DGE results
+
+The Bioconductor package [`DEGreport`](https://bioconductor.org/packages/release/bioc/html/DEGreport.html) can use the DESeq2 results output to make the top20 genes and the volcano plots generated above by writing much fewer lines of code. The caveat of these functions is you lose the ability to customize plots as we have demonstrated above.
+
+If you are interested, the example code below shows how you can use DEGreport to create similar plots. **Note that this is example code, do not run.**
+
 > ``` r
 > ## load degreport
 > library(degreport)
@@ -403,6 +419,35 @@ Take the above steps (starting with finding the Ensembl Gene ID) for a gene symb
 > DEGreport::degPlotWide(dds = dds, genes = row.names(res)[1:5], group = "condition")
 > ```
 
+## Summary script from this lesson
+
+Today, we didn't create any new data objects required for the next lesson. However, if you want to save the commands we used to generate the custom data frames for plotting, you can save relevant lines from below:
+
+``` r
+## Obtain logical vector where TRUE values denote padj values < 0.05 and fold change > 1.5 in either direction (meaning log2FC >= 0.58)
+significant_OE <- res_OE_df$padj < 0.1 & abs(res_OE_df$log2FoldChange) >= 0.58
+
+## Add this vector as a new column to create version of res_OE_df for custom plots
+res_OE_df_plotting <- cbind(res_OE_df, significant_OE)
+
+# Create an empty column
+res_OE_df_plotting$genelabels <- ""
+
+# Sort by padj values
+res_OE_df_plotting <- res_OE_df_plotting[order(res_OE_df_plotting$padj), ]
+
+## Populate the genelabels column with contents of the gene symbols column for the first 10 rows, i.e. the top 10 most significantly expressed genes
+res_OE_df_plotting$genelabels[1:10] <- as.character(res_OE_df_plotting$symbol[1:10])
+
+# Create a new empty column
+res_OE_df_plotting$my_genelabels <- ""
+
+# Assign value just to that one empty genelabel slot where symbol == HOXA3
+res_OE_df_plotting[res_OE_df_plotting$symbol == "HOXA3", "genelabels"] = "HOXA3"
+```
+
+You may also want to save the commands for creating the plots themselves!
+
 ------------------------------------------------------------------------
 
 *This lesson has been developed by members of the teaching team at the [Harvard Chan Bioinformatics Core (HBC)](http://bioinformatics.sph.harvard.edu/). These are open access materials distributed under the terms of the [Creative Commons Attribution license](https://creativecommons.org/licenses/by/4.0/) (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.*