You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: lessons/wk7_lesson02_wald_test.md
+38-24Lines changed: 38 additions & 24 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -1,7 +1,7 @@
1
1
---
2
2
title: "Exploring DESeq2 results: Wald test"
3
3
author: "Harvard HPC Staff, adapted by Sally Chang at NICHD"
4
-
edited: "Last Modified March 2025"
4
+
edited: "Last Modified May 2025"
5
5
---
6
6
7
7
Approximate time: 60 minutes
@@ -12,6 +12,36 @@ Approximate time: 60 minutes
12
12
- Summarize the different levels of gene filtering
13
13
- Explain log fold change shrinkage
14
14
15
+
## Catch-Up Script
16
+
17
+
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.
18
+
19
+
```r
20
+
# Setup
21
+
# Bioconductor and CRAN libraries used - already installed on Biowulf
Remember to get your HPC On Demand session going, if applicable, and open your `DEAnalysis` R project!
44
+
15
45
# Exploring Results (Wald test)
16
46
17
47
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
56
86
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:
57
87
58
88
```r
59
-
## Define contrasts for MOV10 overexpression
89
+
## Define contrasts for MOV10 overexpression - please run this!
> ### Does it matter what I choose to be my base level?
64
94
>
65
95
> 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.
66
96
67
-
## Getting Ready to work in R
68
-
69
-
1.Get your HPC On Demand session going:
70
-
71
-
- 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`
72
-
73
-
- 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.
74
-
75
-
- Using the Project menu in the top right corner, or the Files Pane window (clicking rnaseq -\> DEanalysis), to navigate to and open `DEanalysis.Rproj`
76
-
77
-
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!
78
-
79
-
3. Run the actual DESeq2 analysis if needed `dds <- DESeq(dds)`.
80
-
81
97
## The results table
82
98
83
99
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:
@@ -329,11 +345,11 @@ Now that we have results for the overexpression results, do the same for the **C
329
345
3. Shrink the LFC estimates using `lfcShrink()` and assign it back to `res_tableKD`.
330
346
4. Store the code you used to do the above exercises in an RScript named `control_vs_kd.R` .
331
347
332
-
## R Script updates:
348
+
## R Script updates:
333
349
334
350
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:
335
351
336
-
```r
352
+
```r
337
353
# Gene-level differential expression analysis using DESeq2
*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.*
374
-
375
-
*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.*|
391
+
|*Some materials and hands-on activities were adapted from [RNA-seq workflow](http://www.bioconductor.org/help/workflows/rnaseqGene/#de) on the Bioconductor website*|
0 commit comments