adding Fig1 R code

sekalylab · Sep 3, 2018 · 92e6209 · 92e6209
1 parent 7f36d13
commit 92e6209
Show file tree

Hide file tree

Showing 14 changed files with 222 additions and 8 deletions.
diff --git a/README.md b/README.md
@@ -5,19 +5,20 @@
 [Fig. 2](#fig-2), [Fig. 3](#fig-3), [Fig. 4](#fig-4), [Fig. 5](#fig-5), [Fig. 6](#fig-6)
 
 ### Fig. 2
-![Fig. 2](figure/20150201_RV144pilot.Fig1.png)
+![Fig. 2](figure/20150201_RV144pilot.Fig2.png)
+Fig. 1: [R code [MD]](code/20160510_RV144pilot.Fig2.code.md), [Input file [RData]](output/rv144pilot.gsSetVehSubstracted.RData) 
 
 ### Fig. 3
-![Fig. 3](figure/20150201_RV144pilot.Fig2.png)
+![Fig. 3](figure/20150201_RV144pilot.Fig3.png)
 
 ### Fig. 4
-![Fig. 4](figure/20150201_RV144pilot.Fig3.png)
+![Fig. 4](figure/20150201_RV144pilot.Fig4.png)
 
 ### Fig. 5
-![Fig. 5](figure/20150201_RV144pilot.Fig4.png)
+![Fig. 5](figure/20150201_RV144pilot.Fig5.png)
 
 ### Fig. 6
-![Fig. 6](figure/20170127_RV144.Fig5.png)
+![Fig. 6](figure/20170127_RV144.Fig6.png)
 
 ## Supplemental material:
 

diff --git a/code/20160509_RV144pilot.preprocessing.code.md b/code/20160509_RV144pilot.preprocessing.code.md
@@ -1,7 +1,7 @@
 ---
 title: RV144 pilot gene-expression analysis
 author: Slim Fourati
-date: May 9 2016
+date: May, 9 2016
 output: github_documents
 ---
 

diff --git a/code/20160510_RV144pilot.Fig2.code.md b/code/20160510_RV144pilot.Fig2.code.md
@@ -0,0 +1,214 @@
+---
+title: R code to generate Figure 2
+author: Slim Fourati
+date: May 10, 2016
+output: github_documents
+---
+
+loading require packages
+
+```r
+suppressPackageStartupMessages(library(package = "knitr"))
+suppressPackageStartupMessages(library(package = "RCurl"))
+suppressPackageStartupMessages(library(package = "Biobase"))
+suppressPackageStartupMessages(library(package = "ggplot2"))
+suppressPackageStartupMessages(library(package = "cluster"))
+suppressPackageStartupMessages(library(package = "pheatmap"))
+suppressPackageStartupMessages(library(package = "readr"))
+suppressPackageStartupMessages(library(package = "dplyr"))
+suppressPackageStartupMessages(library(package = "tidyr"))
+suppressPackageStartupMessages(library(package = "tibble"))
+```
+
+set default options/variables
+
+```r
+workDir <- dirname(getwd())
+opts_chunk$set(tidy = FALSE, fig.path = "../figure/")
+options(stringsAsFactors = FALSE, width = 80)
+```
+
+read SLEA ExpressionSet
+
+```r
+gsSetRawFile <- file.path(workDir,
+                          "output/rv144pilot.gsSetVehSubstracted.RData")
+load(file = gsSetRawFile)
+```
+
+plot barplot showing normalized enrichment score of enriched genesets
+
+```r
+plotDF <- fData(gsSetRaw) %>%
+  filter(`FDR q-val` <= 0.05) %>%
+  arrange(NES, `FDR q-val`) %>%
+  mutate(NAME = factor(NAME, levels = NAME),
+         `FDR q-val` = pmax(`FDR q-val`, 1e-04),
+         `FDR q-val` = -log10(`FDR q-val`) * sign(NES)) %>%
+  select(NAME, NES, `FDR q-val`)
+plotBar <- ggplot(data    = plotDF,
+                  mapping = aes(x = NAME, y = NES, fill = `FDR q-val`)) +
+  geom_bar(stat = "identity") +
+  scale_y_continuous(expand = c(0.01, 0)) +
+  scale_fill_gradient2(name = "-log10(FDR q-val)",
+                       low = "blue", mid = "white", high = "red") +
+  labs(x = NULL, y = "GSEA: NES") +
+  geom_hline(yintercept = 0) +
+  coord_flip() +
+  theme_bw() +
+  theme(axis.line.x  = element_line(color = "black"),
+        axis.line.y  = element_blank(),
+        axis.text.y  = element_text(hjust = 1, vjust = 0.5),
+        axis.ticks.y = element_blank(),
+        panel.border = element_blank(),
+        panel.grid   = element_blank())
+print(plotBar)
+```
+
+![plot of chunk fig2a](../figure/fig2a-1.png)
+
+determine the number of set of correlated genesets
+
+```r
+flag <- pData(gsSetRaw) %>%
+  filter(`stimulation time` %in% 15) %>%
+  select(donor, visit, `Sample name`) %>%
+  spread(visit, `Sample name`) %>%
+  filter(complete.cases(.)) %>%
+  gather(visit, `Sample name`, -donor)
+mat <- exprs(gsSetRaw)[fData(gsSetRaw)$"FDR q-val" <= 0.05,
+                       flag$"Sample name"]
+hclustFun <- function(X, k) {
+  distMat <- as.dist((1-cor(t(X)))/2)
+  hc <- hclust(distMat)
+  group <- cutree(hc, k)
+  return(value = list(cluster = group))
+}
+set.seed(seed = 111)
+fit <- clusGap(mat, hclustFun, K.max = 10)
+plot(fit)
+```
+
+![plot of chunk figs1](../figure/figs1-1.png)
+
+```r
+print(fit, method = "globalSEmax")
+```
+
+```
+## Clustering Gap statistic ["clusGap"] from call:
+## clusGap(x = mat, FUNcluster = hclustFun, K.max = 10)
+## B=100 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
+##  --> Number of clusters (method 'globalSEmax', SE.factor=1): 4
+##           logW   E.logW         gap     SE.sim
+##  [1,] 4.332752 4.286456 -0.04629572 0.09489605
+##  [2,] 4.113765 3.979655 -0.13411004 0.12651150
+##  [3,] 3.903299 3.727155 -0.17614382 0.11567892
+##  [4,] 3.208124 3.504154  0.29603045 0.11518557
+##  [5,] 3.065703 3.267838  0.20213479 0.11898030
+##  [6,] 2.843997 3.035234  0.19123737 0.14227075
+##  [7,] 2.606959 2.743604  0.13664472 0.16986437
+##  [8,] 2.306303 2.389217  0.08291446 0.20388936
+##  [9,] 2.034552 1.882376 -0.15217670 0.24243110
+## [10,] 1.294605 1.106186 -0.18841879 0.33161430
+```
+
+```r
+Kmax <- maxSE(fit$Tab[, "gap"], fit$Tab[, "SE.sim"], method = "globalSEmax")
+```
+
+correlation heatmap
+
+```r
+corMat <- cor(t(mat))
+distMat <- as.dist(1 - corMat)
+pheatmap(mat                      = corMat,
+         color                    = colorRampPalette(c("white",
+                                                       "yellow",
+                                                       "orange",
+                                                       "red"))(100),
+         breaks                   = c(min(corMat),
+                                      seq(from = 0,
+                                          to = 1,
+                                          length.out = 100)),
+         cellwidth                = 8,
+         cellheight               = 8,
+         clustering_distance_rows = distMat,
+         clustering_distance_cols = distMat,
+         treeheight_row           = 0)
+```
+
+![plot of chunk fig2b](../figure/fig2b-1.png)
+
+heatmap geneset by pathway
+
+```r
+colOrder <- apply(mat, 1, rank) %>%
+  rowMeans() %>%
+  data.frame(meanRank = .) %>%
+  rownames_to_column() %>%
+  merge(pData(gsSetRaw), by.y = "Sample name", by.x = "rowname") %>%
+  arrange(treatment, desc(visit), meanRank)
+pheatmap(mat                      = mat[, colOrder$rowname],
+         color                    = colorRampPalette(c("blue",
+                                                       "white",
+                                                       "red"))(100),
+         breaks                   = c(-max(abs(mat)),
+                                      seq(from       = -1 * min(abs(range(mat))),
+                                          to         = min(abs(range(mat))),
+                                          length.out = 99),
+                                      max(abs(mat))),
+         cellwidth                = 3,
+         cellheight               = 8,
+         fontsize                 = 8,
+         cluster_cols             = FALSE,
+         clustering_distance_rows = "correlation",
+         annotation_col           = pData(gsSetRaw)[, c("visit", "treatment")],
+         annotation_color         = list(treatment = c(VACCINE = "green",
+                                                       PLACEBO = "blue"),
+                                         visit     = c(pre  = "purple",
+                                                       post = "orange")),
+         show_colnames            = FALSE)
+```
+
+![plot of chunk fig2c](../figure/fig2c-1.png)
+
+print session info
+
+```r
+sessionInfo()
+```
+
+```
+## R version 3.5.1 (2018-07-02)
+## Platform: x86_64-apple-darwin17.6.0 (64-bit)
+## Running under: macOS High Sierra 10.13.6
+## 
+## Matrix products: default
+## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
+## LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib
+## 
+## locale:
+## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
+## 
+## attached base packages:
+## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
+## [8] base     
+## 
+## other attached packages:
+##  [1] bindrcpp_0.2.2      tibble_1.4.2        tidyr_0.8.1        
+##  [4] dplyr_0.7.6         readr_1.1.1         pheatmap_1.0.10    
+##  [7] cluster_2.0.7-1     ggplot2_3.0.0       Biobase_2.40.0     
+## [10] BiocGenerics_0.26.0 RCurl_1.95-4.11     bitops_1.0-6       
+## [13] knitr_1.20         
+## 
+## loaded via a namespace (and not attached):
+##  [1] Rcpp_0.12.18       highr_0.7          pillar_1.3.0       compiler_3.5.1    
+##  [5] RColorBrewer_1.1-2 plyr_1.8.4         bindr_0.1.1        tools_3.5.1       
+##  [9] digest_0.6.15      evaluate_0.11      gtable_0.2.0       pkgconfig_2.0.1   
+## [13] rlang_0.2.1        withr_2.1.2        stringr_1.3.1      hms_0.4.2         
+## [17] grid_3.5.1         tidyselect_0.2.4   glue_1.3.0         R6_2.2.2          
+## [21] purrr_0.2.5        magrittr_1.5       scales_1.0.0       assertthat_0.2.0  
+## [25] colorspace_1.3-2   labeling_0.3       stringi_1.2.4      lazyeval_0.2.1    
+## [29] munsell_0.5.0      crayon_1.3.4
+```
diff --git a/code/20160510_RV144pilot.geneset_analysis.code.md b/code/20160510_RV144pilot.geneset_analysis.code.md
@@ -1,7 +1,7 @@
 ---
 title: RV144 pilot gene-expression geneset analysis
 author: Slim Fourati
-date: May 10 2016
+date: May 10, 2016
 output: github_documents
 ---
 
@@ -417,4 +417,3 @@ sessionInfo()
 ## [29] colorspace_1.3-2   utf8_1.1.4         stringi_1.2.4      munsell_0.5.0     
 ## [33] crayon_1.3.4
 ```
-\end{document}
diff --git a/figure/20150201_RV144pilot.Fig1.png b/figure/20150201_RV144pilot.Fig1.png
diff --git a/figure/20150201_RV144pilot.Fig2.png b/figure/20150201_RV144pilot.Fig2.png
diff --git a/figure/20150201_RV144pilot.Fig3.png b/figure/20150201_RV144pilot.Fig3.png
diff --git a/figure/20150201_RV144pilot.Fig4.png b/figure/20150201_RV144pilot.Fig4.png
diff --git a/figure/20150201_RV144pilot.Fig5.png b/figure/20150201_RV144pilot.Fig5.png
diff --git a/figure/20170127_RV144.Fig5.png → figure/20170127_RV144.Fig6.png b/figure/20170127_RV144.Fig5.png → figure/20170127_RV144.Fig6.png
diff --git a/figure/fig2a-1.png b/figure/fig2a-1.png
diff --git a/figure/fig2b-1.png b/figure/fig2b-1.png
diff --git a/figure/fig2c-1.png b/figure/fig2c-1.png
diff --git a/figure/figs1-1.png b/figure/figs1-1.png