Fix comparison against CMap perturbations using GSEA

nuno-agostinho · Nov 11, 2018 · fc66327 · fc66327
1 parent ecb24eb
commit fc66327
Show file tree

Hide file tree

Showing 5 changed files with 17 additions and 5 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: cTRAP
 Title: Identification of candidate causal perturbations from differential gene 
     expression data
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(
         person(c("Bernardo", "P."), "de Almeida", role="aut"),
         person("Nuno", "Saraiva-Agostinho",
@@ -25,7 +25,7 @@ Suggests: testthat,
     knitr,
     covr,
     biomaRt
-RoxygenNote: 6.1.0
+RoxygenNote: 6.1.1
 Imports: data.table,
     limma,
     stats,

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,8 @@
+# 1.0.2 (11 November, 2018)
+
+* Fix comparison against CMap perturbations using gene set enrichment analysis 
+(the resulting score was the additive inverse of the real scores)
+
 # 1.0.1 (2 November, 2018)
 
 * Update title, author names, version and README

diff --git a/R/L1000.R b/R/L1000.R
@@ -233,7 +233,7 @@ compareAgainstL1000 <- function(diffExprGenes, perturbations, cellLine,
                               perturbations, method, pAdjustMethod)
         colnameSuffix <- sprintf("_%s_coef", method)
     } else if (method == "gsea") {
-        ordered     <- order(diffExprGenes)
+        ordered     <- order(diffExprGenes, decreasing=TRUE)
         topGenes    <- names(diffExprGenes)[head(ordered, geneSize)]
         bottomGenes <- names(diffExprGenes)[tail(ordered, geneSize)]
         gsc <- loadGSC(matrix(c(

diff --git a/tests/testthat/test_L1000.R b/tests/testthat/test_L1000.R
@@ -50,4 +50,11 @@ test_that("Compare using GSEA", {
                                 method="gsea")
     expect_is(data, "l1000comparison")
     expect_identical(colnames(data), c("genes", "HepG2_WTCS", "Average_WTCS"))
+    expect_identical(head(data$genes),
+                     c("BRD-A14014306", "BRD-A65142661", "BRD-K31030218",
+                       "BRD-K41172353", "BRD-K77508012", "BRD-K84389640"))
+    expect_identical(head(data[order(data$HepG2_WTCS), ]$genes),
+                     c("BRD-A65142661", "caffeic-acid-phenethyl-ester",
+                       "BRD-K94818765", "BRD-K77508012", "BRD-K31030218",
+                       "BRD-K41172353"))
 })
diff --git a/vignettes/comparing_DGE_with_perturbations.Rmd b/vignettes/comparing_DGE_with_perturbations.Rmd
@@ -219,7 +219,7 @@ compareKnockdown$spearman_ordered <- compareKnockdown$spearman[
 compareKnockdown$pearson_ordered <- compareKnockdown$pearson[
     order(compareKnockdown$pearson$HepG2_pearson_coef, decreasing=TRUE)]
 compareKnockdown$gsea_ordered <- compareKnockdown$gsea[
-    order(compareKnockdown$gsea$HepG2_WTCS, decreasing=FALSE)]
+    order(compareKnockdown$gsea$HepG2_WTCS, decreasing=TRUE)]
 
 # Most positively associated perturbations (note that EIF4G1 knockdown is the
 # 6th, 1st and 2nd most positively associated perturbation based on Spearman
@@ -243,7 +243,7 @@ compareSmallMolecule$spearman_ordered <- compareSmallMolecule$spearman[
 compareSmallMolecule$pearson_ordered <- compareSmallMolecule$pearson[
     order(compareSmallMolecule$pearson$HepG2_pearson_coef, decreasing=TRUE)]
 compareSmallMolecule$gsea_ordered <- compareSmallMolecule$gsea[
-    order(compareSmallMolecule$gsea$HepG2_WTCS, decreasing=FALSE)]
+    order(compareSmallMolecule$gsea$HepG2_WTCS, decreasing=TRUE)]
 
 # Most positively associated perturbations
 head(compareSmallMolecule$spearman_ordered)