updates to Deseq and A+V

Author: AbbiEdwards

…/10_DE_analysis_with_DESeq2_part1_IM.Rmd …/10_DE_analysis_with_DESeq2_part1_IM.RmdMarkdowns/10_DE_analysis_with_DESeq2_part1_IM.Rmd renamed to InstructorNotes/10_DE_analysis_with_DESeq2_part1_IM.Rmd Markdowns/10_DE_analysis_with_DESeq2_part1_IM.Rmd renamed to InstructorNotes/10_DE_analysis_with_DESeq2_part1_IM.Rmd

@@ -3,8 +3,8 @@ title: "Introduction to Bulk RNAseq data analysis"
 subtitle: Differential Expression of RNA-seq data
 date: '`r format(Sys.time(), "Last modified: %d %b %Y")`'
 output:
-  pdf_document: default
   html_document: default
+  pdf_document: default
 ---
 
 ```{r setup, echo=FALSE, message=FALSE}
@@ -13,9 +13,42 @@ options(tibble.print_max = 4, tibble.print_min = 4, max.print=40,
 library(tidyverse)
 library(DESeq2)
 ```
-
 ### Exercise 1
 
+> Now we have made our results table using our simple model, let have a look at which
+> genes are changing and how many pass our 0.05 threshold. Why might this not be
+> straigtforward?
+>
+>
+
+```{r eval=FALSE}
+sum(results.simple$padj < 0.05)
+```
+
+```{r eval=FALSE}
+sum(is.na(results.simple$padj))
+```
+>
+> a) how many genes are significantly (with an FDR < 0.05) up-regulated?
+>
+>
+>
+
+```{r eval=FALSE}
+sum(results.simple$padj < 0.05 & results.simple$log2FoldChange > 0, na.rm = TRUE)
+```
+>
+> b) how many genes are significantly (with an FDR < 0.05) down-regulated?
+>
+>
+
+```{r eval=FALSE}
+sum(results.simple$padj < 0.05 & results.simple$log2FoldChange < 0, na.rm = TRUE)
+
+```
+
+### Exercise 2
+
 > So far we have fitted a simple model considering just "Status", but in reality 
 > we want to model the effects of both "Status" and "Time Point".
 > 
@@ -84,7 +117,7 @@ sum(results.additive$padj < 0.05, na.rm = TRUE)
 ```
 
 
-### Exercise 2
+### Exercise 3
 
 > If we want a different contrast we can just pass the `results` function the 
 > **name** of the contrast, as given by `resultsNames(ddsObj)`.
@@ -100,7 +133,7 @@ results.d33vd11 <- results(ddsObj, name= "TimePoint_d33_vs_d11", alpha=0.05)
 sum(results.d33vd11$padj < 0.05, na.rm = TRUE)
 ```
 
-### Exercise 3
+### Exercise 4
 >
 > When we looked at the PCA it did seem that an interaction model might be
 > warranted. Let's test that.  
@@ -143,7 +176,7 @@ results.Interaction_v_Additive <- results(ddsObj.LRT)
 table(results.Interaction_v_Additive$padj < 0.05)
 ```
 
-### Exercise 4
+### Exercise 5
 >
 > Let's investigate the uninfected mice
 >  

…10_DE_analysis_with_DESeq2_part1_IM.html …10_DE_analysis_with_DESeq2_part1_IM.htmlMarkdowns/10_DE_analysis_with_DESeq2_part1_IM.html renamed to InstructorNotes/10_DE_analysis_with_DESeq2_part1_IM.html Markdowns/10_DE_analysis_with_DESeq2_part1_IM.html renamed to InstructorNotes/10_DE_analysis_with_DESeq2_part1_IM.html

@@ -165,14 +165,34 @@
 
 <h1 class="title toc-ignore">Introduction to Bulk RNAseq data analysis</h1>
 <h3 class="subtitle">Differential Expression of RNA-seq data</h3>
-<h4 class="date">Last modified: 19 Apr 2021</h4>
+<h4 class="date">Last modified: 28 Nov 2021</h4>
 
 </div>
 
 
 <div id="exercise-1" class="section level3">
 <h3>Exercise 1</h3>
 <blockquote>
+<p>Now we have made our results table using our simple model, let have a look at which genes are changing and how many pass our 0.05 threshold. Why might this not be straigtforward?</p>
+</blockquote>
+<pre class="r"><code>sum(results.simple$padj &lt; 0.05)</code></pre>
+<pre class="r"><code>sum(is.na(results.simple$padj))</code></pre>
+<blockquote>
+<ol style="list-style-type: lower-alpha">
+<li>how many genes are significantly (with an FDR &lt; 0.05) up-regulated?</li>
+</ol>
+</blockquote>
+<pre class="r"><code>sum(results.simple$padj &lt; 0.05 &amp; results.simple$log2FoldChange &gt; 0, na.rm = TRUE)</code></pre>
+<blockquote>
+<ol start="2" style="list-style-type: lower-alpha">
+<li>how many genes are significantly (with an FDR &lt; 0.05) down-regulated?</li>
+</ol>
+</blockquote>
+<pre class="r"><code>sum(results.simple$padj &lt; 0.05 &amp; results.simple$log2FoldChange &lt; 0, na.rm = TRUE)</code></pre>
+</div>
+<div id="exercise-2" class="section level3">
+<h3>Exercise 2</h3>
+<blockquote>
 <p>So far we have fitted a simple model considering just “Status”, but in reality we want to model the effects of both “Status” and “Time Point”.</p>
 <p>Let’s start with the model with only main effects - an additive model with no interaction. The main assumption here is that the effects of Status and the effects of Time Point are indepedent.</p>
 <p>Recapitulate the above steps to generate a new DESeq2 object with the additive model. Then we will extract the results table as above.</p>
@@ -248,8 +268,8 @@ <h4 id="filter-the-data-set">Filter the data set:</h4>
 <pre class="r"><code>sum(results.additive$padj &lt; 0.05, na.rm = TRUE)</code></pre>
 <pre><code>## [1] 2766</code></pre>
 </div>
-<div id="exercise-2" class="section level3">
-<h3>Exercise 2</h3>
+<div id="exercise-3" class="section level3">
+<h3>Exercise 3</h3>
 <blockquote>
 <p>If we want a different contrast we can just pass the <code>results</code> function the <strong>name</strong> of the contrast, as given by <code>resultsNames(ddsObj)</code>. Look at the help page for the <code>results</code> command to see how to do this.</p>
 <ol style="list-style-type: decimal">
@@ -265,8 +285,8 @@ <h3>Exercise 2</h3>
 <pre class="r"><code>sum(results.d33vd11$padj &lt; 0.05, na.rm = TRUE)</code></pre>
 <pre><code>## [1] 109</code></pre>
 </div>
-<div id="exercise-3" class="section level3">
-<h3>Exercise 3</h3>
+<div id="exercise-4" class="section level3">
+<h3>Exercise 4</h3>
 <blockquote>
 <p>When we looked at the PCA it did seem that an interaction model might be warranted. Let’s test that.</p>
 <ol style="list-style-type: decimal">
@@ -323,8 +343,8 @@ <h3>Exercise 3</h3>
 ## FALSE  TRUE 
 ## 16474   455</code></pre>
 </div>
-<div id="exercise-4" class="section level3">
-<h3>Exercise 4</h3>
+<div id="exercise-5" class="section level3">
+<h3>Exercise 5</h3>
 <blockquote>
 <p>Let’s investigate the uninfected mice</p>
 <ol style="list-style-type: decimal">

Markdowns/10_DE_analysis_with_DESeq2.Solutions.Rmd

@@ -3,10 +3,10 @@ title: "Introduction to Bulk RNAseq data analysis"
 subtitle: Differential Expression of RNA-seq data - Part 2
 date: '`r format(Sys.time(), "Last modified: %d %b %Y")`'
 output:
-  html_document:
-    toc: yes
   pdf_document:
     toc: yes
+  html_document:
+    toc: yes
 bibliography: ref.bib
 ---
 
@@ -15,7 +15,33 @@ options(tibble.print_max = 4, tibble.print_min = 4, max.print=40,
         tibble.max_extra_cols=2)
 ```
 
-### Exercise 1
+# Recap from last week
+
+```{r message = FALSE}
+library(DESeq2)
+library(tidyverse)
+```
+
+```{r loadData}
+txi <- readRDS("RObjects/txi.rds")
+sampleinfo <- read_tsv("data/samplesheet_corrected.tsv", col_types="cccc") %>%
+  mutate(Status = fct_relevel(Status, "Uninfected"))
+
+simple.model <- as.formula(~ Status)
+
+ddsObj.raw <- DESeqDataSetFromTximport(txi = txi,
+                                       colData = sampleinfo,
+                                       design = simple.model)
+
+keep <- rowSums(counts(ddsObj.raw)) > 5
+ddsObj.filt <- ddsObj.raw[keep,]
+
+ddsObj <- DESeq(ddsObj.filt)
+
+results.simple <- results(ddsObj, alpha=0.05)
+```
+
+### Exercise 2
 
 > So far we have fitted a simple model considering just "Status", but in reality 
 > we want to model the effects of both "Status" and "Time Point".
@@ -115,7 +141,7 @@ topGenesIvU <- as.data.frame(results.InfectedvUninfected) %>%
 topGenesIvU
 ```
 
-### Exercise 2
+### Exercise 3
 
 > If we want a different contrast we can just pass the `results` function the 
 > **name** of the contrast, as given by `resultsNames(ddsObj)`.
@@ -191,7 +217,7 @@ to all genes. Curiously then, this suggests that overall the simple model
 is more appropriate than the additive model. Let's look into the interaction
 model.
 
-### Exercise 3
+### Exercise 4
 >
 > When we looked at the PCA it did seem that an interaction model might be
 > warranted. Let's test that.  
@@ -283,7 +309,7 @@ sum(results.interaction.33$padj < 0.05, na.rm = TRUE)
 We can see that there is a strong difference in the effects of infection on 
 gene expression between days 11 and 33.
 
-### Exercise 4
+### Exercise 5
 >
 > Let's investigate the uninfected mice
 >