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

esallychang · web-flow · commit 26e0caee1415 · 2025-05-09T15:41:29.000-04:00
Week 7 01-02: Catch-up scripts and remove dplyr pipes
diff --git a/lessons/wk6_lesson05_deseq2_analysis.md b/lessons/wk6_lesson05_deseq2_analysis.md
@@ -1,7 +1,7 @@
 ---
 title: "Gene-level differential expression analysis with DESeq2"
 author: "Harvard HPC Staff, Adapted by Sally Chang at NICHD"
-date: "Last Modified March 2025"
+date: "Last Modified May 2025"
 ---
 
 Approximate time: 60 minutes
@@ -13,6 +13,31 @@ Approximate time: 60 minutes
 -   Inspect gene-level dispersion estimates
 -   Recognize the importance of dispersion during differential expression analysis
 
+## Catch-Up Script
+
+If you need to be completely caught up, you can copy and paste the following into an R Script and run it. If you don't already have the files in your `/data` directory, please see [Wk 5 Lesson 01](../wk5_lesson01_introR_Rstudio.md) for instructions on where to obtain the input files.
+
+``` r
+# Setup
+# Bioconductor and CRAN libraries used - already installed on Biowulf
+library(tidyverse)
+library(RColorBrewer)
+library(DESeq2)
+library(pheatmap)
+library(BiocManager)
+
+# Load in data
+data <- read.table("data/mov10_AllSamples_featurecounts.Rmatrix.txt", header=T, row.names=1)
+
+meta <- read.table("data/mov10_AllSamples_metadata.txt", header=T, row.names=1)
+
+# Create DESeq2Dataset object
+dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = ~ sampletype) 
+
+# Run DESeq2 on DESeq2Dataset object
+dds <- DESeq(dds)
+```
+
 ## DESeq2 differential gene expression analysis workflow
 
 Previously, we created the DESeq2 object using the appropriate design formula and running DESeq2 using the two lines of code:
diff --git a/lessons/wk7_lesson01_hypothesis_testing.md b/lessons/wk7_lesson01_hypothesis_testing.md
@@ -1,7 +1,7 @@
 ---
 title: "Hypothesis testing and multiple testing"
 author: "Harvard HPC Staff, Adapted by Sally Chang @ NICHD"
-date: "Last Modified February 2025"
+date: "Last Modified May 2025"
 ---
 
 Approximate time: 60 minutes
@@ -13,13 +13,40 @@ Approximate time: 60 minutes
 -   Recognize the importance of multiple test correction
 -   Identify different methods for multiple test correction
 
+## Catch-Up Script
+
+If you need to be completely caught up, you can copy and paste the following into an R Script and run it. If you don't already have the files in your `/data` directory, please see [Wk 5 Lesson 01](../wk5_lesson01_introR_Rstudio.md) for instructions on where to obtain the input files.
+
+``` r
+# Setup
+# Bioconductor and CRAN libraries used - already installed on Biowulf
+library(tidyverse)
+library(RColorBrewer)
+library(DESeq2)
+library(pheatmap)
+library(BiocManager)
+
+# Load in data
+data <- read.table("data/mov10_AllSamples_featurecounts.Rmatrix.txt", header=T, row.names=1)
+
+meta <- read.table("data/mov10_AllSamples_metadata.txt", header=T, row.names=1)
+
+# Create DESeq2Dataset object
+dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = ~ sampletype) 
+
+# Run DESeq2 on DESeq2Dataset object
+dds <- DESeq(dds)
+```
+
+Remember to get your HPC On Demand session going, if applicable, and open your `DEAnalysis` R project!
+
 # DESeq2: Model fitting and Hypothesis testing
 
 The final step in the DESeq2 workflow is taking the counts for each gene and fitting it to the model and testing for differential expression.
 
 <p align="center">
 
-<img src="../img/deseq_workflow_full_2018.png" width="385" alt="deseq full workflow"/>
+<img src="../img/deseq_workflow_full_2018.png" alt="deseq full workflow" width="385"/>
 
 </p>
 
@@ -151,6 +178,34 @@ DESeq2 helps reduce the number of genes tested by removing those genes unlikely
 
 By setting the FDR cutoff to \< 0.05, we're saying that the proportion of false positives we expect amongst our differentially expressed genes is 5%. For example, if you call 500 genes as differentially expressed with an FDR cutoff of 0.05, you expect 25 of them to be false positives.
 
+## Your DE script
+
+In this lesson, we took the additional step of running an additional DESeq2 analysis using a Likelihood Ratio Test to create the `dds_lrt` object . Your `de_script.R` should now contain the following commands to re-create necessary data objects (click to show):
+
+``` r
+# Setup
+# Bioconductor and CRAN libraries used - already installed on Biowulf
+library(tidyverse)
+library(RColorBrewer)
+library(DESeq2)
+library(pheatmap)
+library(BiocManager)
+
+# Load in data
+data <- read.table("data/mov10_AllSamples_featurecounts.Rmatrix.txt", header=T, row.names=1)
+
+meta <- read.table("data/mov10_AllSamples_metadata.txt", header=T, row.names=1)
+
+# Create DESeq2Dataset object
+dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = ~ sampletype) 
+
+# Run DESeq2 on DESeq2Dataset object
+dds <- DESeq(dds)
+
+# Likelihood ratio test
+dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)
+```
+
 ------------------------------------------------------------------------
 
 *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.*
diff --git a/lessons/wk7_lesson02_wald_test.md b/lessons/wk7_lesson02_wald_test.md
@@ -1,7 +1,7 @@
 ---
 title: "Exploring DESeq2 results: Wald test"
 author: "Harvard HPC Staff, adapted by Sally Chang at NICHD"
-edited: "Last Modified March 2025"
+edited: "Last Modified May 2025"
 ---
 
 Approximate time: 60 minutes
@@ -12,6 +12,36 @@ Approximate time: 60 minutes
 -   Summarize the different levels of gene filtering
 -   Explain log fold change shrinkage
 
+## Catch-Up Script
+
+If you need to be completely caught up, you can copy and paste the following into an R Script and run it. If you don't already have the files in your `/data` directory, please see [Wk 5 Lesson 01](../wk5_lesson01_introR_Rstudio.md) for instructions on where to obtain the input files.
+
+``` r
+# Setup
+# Bioconductor and CRAN libraries used - already installed on Biowulf
+library(tidyverse)
+library(RColorBrewer)
+library(DESeq2)
+library(pheatmap)
+library(BiocManager)
+
+# Load in data
+data <- read.table("data/mov10_AllSamples_featurecounts.Rmatrix.txt", header=T, row.names=1)
+
+meta <- read.table("data/mov10_AllSamples_metadata.txt", header=T, row.names=1)
+
+# Create DESeq2Dataset object
+dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = ~ sampletype) 
+
+# Run DESeq2 on DESeq2Dataset object
+dds <- DESeq(dds)
+
+# Likelihood ratio test
+dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)
+```
+
+Remember to get your HPC On Demand session going, if applicable, and open your `DEAnalysis` R project!
+
 # Exploring Results (Wald test)
 
 By default DESeq2 uses the Wald test to identify genes that are differentially expressed between two sample classes. Given the factor(s) used in the design formula, and how many factor levels are present, we can extract results for a number of different comparisons. Here, we will walk through how to obtain results from the `dds` object and provide some explanations on how to interpret them.
@@ -56,28 +86,14 @@ Alternatively, if you **only had two factor levels you could do nothing** and no
 To start, we want to evaluate **expression changes between the MOV10 overexpression samples and the control samples**. As such we will use the first method for specifying contrasts and create a character vector:
 
 ``` r
-## Define contrasts for MOV10 overexpression
+## Define contrasts for MOV10 overexpression - please run this!
 contrast_oe <- c("sampletype", "MOV10_overexpression", "control")
 ```
 
 > ### Does it matter what I choose to be my base level?
 >
 > Yes, it does matter. **Deciding what level is the base level will determine how to interpret the fold change that is reported.** So for example, if we observe a log2 fold change of -2 this would mean the gene expression is lower in factor level of interest relative to the base level. Thus, if leaving it up to DESeq2 to decide on the contrasts be sure to check that the alphabetical order coincides with the fold change direction you are anticipating.
 
-## Getting Ready to work in R
-
-1.Get your HPC On Demand session going:
-
--   Opening up RStudio using [HPC on Demand](https://hpcondemand.nih.gov/pun/sys/dashboard/), using default values except for Starting Directory: `/data/Bspc-training/YOUR_USERNAME/rnaseq`
-
--   To check whether or not you are in the correct working directory, use `getwd()`. Something like `/vf/users/Bspc-training/changes/rnaseq` should come up.
-
--   Using the Project menu in the top right corner, or the Files Pane window (clicking rnaseq -\> DEanalysis), to navigate to and open `DEanalysis.Rproj`
-
-2.  We are assuming that you have the `dds` object in your environment and your packages are loaded - run your `de_setup.R` script if needed!
-
-3.  Run the actual DESeq2 analysis if needed `dds <- DESeq(dds)`.
-
 ## The results table
 
 Now that we have our contrast created, we can use it as input to the `results()` function. Let's take a quick look at the help manual for the function:
@@ -102,17 +118,13 @@ The results table that is returned to us is **a `DESeqResults` object**, which i
 class(res_tableOE)
 ```
 
-Now let's take a look at **what information is stored** in the results:
+Now let's take a look at **what information is stored** in the results, using nested functions that convert `res_table0E` into a data frame that we can then View:
 
 ``` r
 # What is stored in results?
-res_tableOE %>% 
-data.frame() %>% 
-View()
+View(data.frame(res_tableOE))
 ```
 
-> **Discussion:** The `%>%` acts as a pipe symbol in R. This functionality comes as part of the [`dplyr`](https://dplyr.tidyverse.org/) package, which was loaded as part of the `tidyverse` that we loaded at the beginning of our lessons. Knowing this, what exactly is the code above doing?
-
 We have six columns of information reported for each gene (row). We can use the `mcols()` function to extract information on what the values stored in each column represent:
 
 ``` r
@@ -151,13 +163,14 @@ The missing values represent genes that have undergone filtering as part of the
 If within a row, all samples have zero counts there is no expression information and therefore these genes are not tested.
 
 ``` r
-# Filter genes by zero expression
-res_tableOE[which(res_tableOE$baseMean == 0),] %>% 
-data.frame() %>% 
-View()
+# Filter genes by zero expression and view using the same type of nested command as above
+View(data.frame(res_tableOE[which(res_tableOE$baseMean == 0),]))
+
+# You could also count the number of rows that are left in the filtered data frame
+nrow(data.frame(res_tableOE[which(res_tableOE$baseMean == 0),]))
 ```
 
-> **The baseMean column for these genes will be zero, and the log2 fold change estimates, p-value and adjusted p-value will all be set to NA. *How would you adjust the command above to count the number of rows matching this condition*?**
+> **The baseMean column for these genes will be zero, and the log2 fold change estimates, p-value and adjusted p-value will all be set to NA.**
 
 **2. Genes with an extreme count outlier**
 
@@ -166,11 +179,9 @@ The `DESeq()` function calculates, for every gene and for every sample, a diagno
 ``` r
 # Filter genes that have an extreme outlier by looking for those rows that have a non-zero base mean but no values for p-value and adjusted p-value. Do we actually have any of these?
 
-res_tableOE[which(is.na(res_tableOE$pvalue) & 
+View(data.frame(res_tableOE[which(is.na(res_tableOE$pvalue) & 
                     is.na(res_tableOE$padj) &
-                    res_tableOE$baseMean > 0),] %>% 
-  data.frame() %>% 
-  View()
+                    res_tableOE$baseMean > 0),]))
 ```
 
 > **If a gene contains a sample with an extreme count outlier then the p-value and adjusted p-value will be set to NA.**
@@ -191,11 +202,9 @@ At a user-specified value (`alpha = 0.1`), DESeq2 evaluates the change in the nu
 
 ``` r
 # Filter genes below the low mean threshold
-res_tableOE[which(!is.na(res_tableOE$pvalue) & 
+View(data.frame(res_tableOE[which(!is.na(res_tableOE$pvalue) & 
                     is.na(res_tableOE$padj) & 
-                    res_tableOE$baseMean > 0),] %>% 
-  data.frame() %>% 
-  View()
+                    res_tableOE$baseMean > 0),]))
 ```
 
 > **If a gene is filtered by independent filtering, then only the adjusted p-value will be set to NA.**
@@ -218,23 +227,9 @@ log2 (normalized_counts_group1 / normalized_counts_group2)
 
 The problem is, these fold change estimates are not entirely accurate as they do not account for the large dispersion we observe with low read counts. To address this, the **log2 fold changes need to be adjusted**.
 
-**This is where we stopped on Tuesday of Week 7!**
-
 ------------------------------------------------------------------------
 
-### More accurate LFC estimates: Picking up again from Tuesday
-
-1.Get your HPC On Demand session going:
-
--   Opening up RStudio using [HPC on Demand](https://hpcondemand.nih.gov/pun/sys/dashboard/), using default values except for Starting Directory and **INCREASE MEMORY TO 8G**: `/data/Bspc-training/YOUR_USERNAME/rnaseq`
-
--   To check whether or not you are in the correct working directory, use `getwd()`. Something like `/vf/users/Bspc-training/changes/rnaseq` should come up.
-
--   Using the Project menu in the top right corner, or the Files Pane window (clicking rnaseq -\> DEanalysis), to navigate to and open `DEanalysis.Rproj`
-
-2.  We are assuming that you have the `dds` object in your environment and your packages are loaded - run your `de_setup.R` script if needed!
-
-3.  Run the actual DESeq2 analysis if needed `dds <- DESeq(dds)`.
+## More accurate LFC estimates
 
 To generate more accurate log2 foldchange (LFC) estimates, DESeq2 allows for the **shrinkage of the LFC estimates toward zero** when the information for a gene is low, which could include:
 
@@ -329,11 +324,11 @@ Now that we have results for the overexpression results, do the same for the **C
 3.  Shrink the LFC estimates using `lfcShrink()` and assign it back to `res_tableKD`.
 4.  Store the code you used to do the above exercises in an RScript named `control_vs_kd.R` .
 
-## R Script updates: 
+## R Script updates:
 
 Here is what your `de_setup.R` script should look like now to regenerate all the R objects we will need to the next lesson:
 
-```r
+``` r
 # Gene-level differential expression analysis using DESeq2
 
 # Setup
@@ -369,9 +364,7 @@ res_tableOE <- lfcShrink(dds, coef="sampletype_MOV10_overexpression_vs_control",
 
 Now, on to [Week 7, lesson 03!](../lessons/wk7_lesson03_summarizing_dge_results.md)
 
-------------------------------------------------------------------------
-*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.*
-
-*Some materials and hands-on activities were adapted from [RNA-seq workflow](http://www.bioconductor.org/help/workflows/rnaseqGene/#de) on the Bioconductor website*
-
-------------------------------------------------------------------------
+|                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
+|------------------------------------------------------------------------|
+| *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.* |
+| *Some materials and hands-on activities were adapted from [RNA-seq workflow](http://www.bioconductor.org/help/workflows/rnaseqGene/#de) on the Bioconductor website*                                                                                                                                                                                                                                                                                               |
