Changes between Version 19 and Version 20 of SOPs/rna-seq-diff-expressions

Timestamp:

09/11/14 15:53:13 (10 years ago)

Author:

gbell

Comment:

--

SOPs/rna-seq-diff-expressions

@@ -112 +112 @@
     * Fold change thresholds may be used in addition, if desired, to select genes that are changing an amount that is biologically meaningful.
 
- * **Remove batch effects**
-    * With [[http://www.bioconductor.org/packages/release/bioc/html/edgeR.html|edgeR]], the treatments can be adjusted for differences between the batch with addictive model formula by design = model.matrix( ~Batch+Treatment). For example, with target.txt like this:
+ * **Accounting for a batch effect in a differential expression model**
+    * With [[http://www.bioconductor.org/packages/release/bioc/html/edgeR.html|edgeR]], the comparison can be adjusted for differences between batches (such as date of sample preparation) with an addictive model formula using the equation 
+
+  '''design = model.matrix(~Batch+Treatment)'''. 
+
+  For example, with a targets file (target.txt) like this:
 
 
    ||= Sample =||= genotype =||= date =|| 
-   ||  KO.1  || KO || old ||
-   ||  KO.2  || KO      || old ||
-   ||  WT.1  || WT      || old ||
-   ||  KO.3  || KO      || new ||
-   ||  KO.4  || KO      || new ||
-   ||  WT.2  || WT      || new ||
+   ||  KO.1  ||  KO  ||  old  ||
+   ||  WT.1  ||  WT  ||  old  ||
+   ||  KO.2  ||  KO  ||  new  ||
+   ||  WT.2  ||  WT  ||  new  ||
 
 
-  Sample R codes for reducing the date effect (old/new) are
+  Sample R code for reducing the date effect (old vs. new) is
 
  {{{
-   target   <- read.delim("target.txt", header=T)
-   genotype <- factor(target$genotype, levels=c("WT", "KO"))
-   mydate   <- factor(target$date, levels=c("old", "new"))
-   Design   <-model.matrix(~mydate+genotype)
-   colnames(Design)
-# [1] "(Intercept)" "mydatenew"   "genotypeKO" 
+        target   <- read.delim("target.txt", header=T)
+        genotype <- factor(target$genotype, levels=c("WT", "KO"))
+        mydate   <- factor(target$date, levels=c("old", "new"))
+        Design   <- model.matrix(~mydate+genotype)
+        colnames(Design)
+        # [1] "(Intercept)" "mydatenew"   "genotypeKO" 
 
-   Y2  <- calcNormFactors(Y) # Y is the DGEList
-   Y2  <- estimateGLMCommonDisp(Y2, Design, verbose=TRUE)
-   Y2  <- estimateGLMTrendedDisp(Y2, Design)
-   Y2  <- estimateGLMTagwiseDisp(Y2, Design)
-   Fit <- glmFit(Y2, Design)
-   Lrt <- glmLRT(Fit, coef="genotypeKO") # coef is pointed to the column name in Design
+        Y2  <- calcNormFactors(Y) # Y is the DGEList
+        Y2  <- estimateGLMCommonDisp(Y2, Design, verbose=TRUE)
+        Y2  <- estimateGLMTrendedDisp(Y2, Design)
+        Y2  <- estimateGLMTagwiseDisp(Y2, Design)
+        Fit <- glmFit(Y2, Design)
+        Lrt <- glmLRT(Fit, coef="genotypeKO") # where 'coef' points to the column name in 'Design'
  }}}
 
-  Lrt$table has the log2FC(Fold change) and P-value.
+  The data structure Lrt$table contains the values for log2FC (log 2 fold change) and p-value.
    
-  Detail information can be found in [[http://www.bioconductor.org/packages/release/bioc/html/edgeR.html|edgeR]]
+  Detailed information can be found in the [[http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf|edgeR User Guide]] (after a search for "batch effects").