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

Week 6 lesson 01-04 scripts

Author: esallychang

lessons/wk6_lesson02_count_normalization.md

@@ -1,21 +1,18 @@
 ---
 title: "Count normalization with DESeq2"
 author: "Harvard HPC Staff, Adapted by Sally Chang @ NICHD"
-date: "Last Modified April 2025"
+date: "Last Modified May 2025"
 ---
 
 Approximate time: 60 minutes
 
-### NOTE: 
-To make names more generalized for the next course, `/data/Bspc-training/shared/rnaseq_jan2025` is now `/data/Bspc-training/shared/rnaseq_mov10` . Make sure to edit any scripts that refer to the shared data! 
-
 ## Learning Objectives
 
 -   Explore different types of normalization methods
 -   Become familiar with the `DESeqDataSet` object
 -   Understand how to normalize counts using DESeq2
 
-### Opening the project using RStudio (HPC on Demand)
+### Preparing for this lesson: 
 
 Before we get into the details of the analysis, let's start by:
 
@@ -25,6 +22,23 @@ Before we get into the details of the analysis, let's start by:
 
 -   Using the Project menu in the top right corner, or the Files Pane window (clicking rnaseq -\> DEanalysis), to navigate to and open `DEanalysis.Rproj`, which you set up as an Assignment last week.
 
+If you missed the last lesson, or need to make sure you have the right packages and data loaded, please run the following commands:
+
+``` 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)
+```
+
 ## Normalization
 
 The first step in the DE analysis workflow is count normalization, which is necessary to make accurate comparisons of gene expression between samples.

lessons/wk6_lesson03_dge_qc_analysis.md

@@ -157,9 +157,9 @@ Now that we have a good understanding of the QC steps normally employed for RNA-
 
 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`
+-   Opening up RStudio using [HPC on Demand](https://hpcondemand.nih.gov/pun/sys/dashboard/), using default values except for Starting Directory: `/data/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.
+-   To check whether or not you are in the correct working directory, use `getwd()`. Something like `/vf/users/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`
 
@@ -223,7 +223,9 @@ plotPCA(rld, intgroup="sampletype")
 ```
 
 <p align="center">
-<img src="../img/mov10_pca.png" width="600" alt="mov 10 experiement PCA"/>
+
+<img src="../img/mov10_pca.png" alt="mov 10 experiement PCA" width="600"/>
+
 </p>
 
 ------------------------------------------------------------------------
@@ -296,7 +298,9 @@ pheatmap(rld_cor, annotation = meta)
 When you plot using `pheatmap()` the hierarchical clustering information is used to place similar samples together and this information is represented by the tree structure along the axes. The `annotation` argument accepts a dataframe as input, in our case it is the `meta` data frame.
 
 <p align="center">
-<img src="../img/mov10_default_heatmap.png" width="600" alt="mov10 heatmap with default settings"/>
+
+<img src="../img/mov10_default_heatmap.png" alt="mov10 heatmap with default settings" width="600"/>
+
 </p>
 
 Overall, we observe pretty high correlations across the board ( \> 0.999) suggesting no outlying sample(s). Also, similar to the PCA plot you see the samples clustering together by sample group. Together, these plots suggest to us that the data are of good quality and we have the green light to proceed to differential expression analysis.

lessons/wk6_lesson04_design_formulas.md

@@ -11,6 +11,28 @@ Approximate time: 30 minutes
 -   Demonstrate the use of the design formula with simple and complex designs
 -   Construct R code to execute the differential expression analysis workflow with DESeq2
 
+## 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) 
+```
+
 # Differential expression analysis with DESeq2
 
 The final step in the differential expression analysis workflow is **fitting the raw counts to the NB model and performing the statistical test** for differentially expressed genes. In this step we essentially want to determine whether the mean expression levels of different sample groups are significantly different.
@@ -168,3 +190,28 @@ Given the the metadata table you have sent me for your own experiment, do the fo
 
 1.  Write a design formula for your experiment, in the format of `design = ~ sex + age + treatment` . Make sure to include any interaction terms or terms that you want to "regress" out. There are additional recommendations for complex designs in the [DESeq2 vignette](https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions).
 2.  Briefly explain (in 1-2 sentences) the reasoning for this design formula.
+
+## Your DE script
+
+In this lesson, we took the additional step of running the actual DESeq2 analysis. 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)
+```

scripts/de_analysis_wk6_lesson01.R

@@ -11,3 +11,8 @@ data <- read.table("data/mov10_AllSamples_featurecounts.Rmatrix.txt", header=T,
 
 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)