Compare revisions

Maxime BORRY · Maxime BORRY · Maxime BORRY · Maxime BORRY · Maxime BORRY · Maxime BORRY
--- a/R/comparResumPlot.r
+++ b/R/comparResumPlot.r
@@ -16,5 +16,5 @@
 comparResumPlot = function(preparedData, savePlot = TRUE, legend = TRUE, directoryPath = directoryPath){
    print("Plotting Summary comparison Plot for all conditions")
    generateSummaryPlots(comparisonSummaryData(preparedData), savePlot = savePlot, legend = legend, file.name = paste(directoryPath, "_comparison_sumplot", sep = ""), format = "pdf")
 }
--- a/R/generateSummaryPlots.r
+++ b/R/generateSummaryPlots.r
@@ -18,10 +18,6 @@
 generateSummaryPlots <- function(plot.data, file.name = "resumPlot", Xlab="-log10(p-value)",
        Ylab="Average Absolute logFC", format = NULL, firstN = 10,  savePlot = TRUE, legend = TRUE){
-        if(!require(ggplot2)){
-            install.packages("ggplot2")
-            require(ggplot2)
-        }
    tryCatch({
        plot.data.sig = plot.data[plot.data[, "rank"] <= firstN, ]
@@ -44,23 +40,23 @@ generateSummaryPlots <- function(plot.data, file.name = "resumPlot", Xlab="-log1
    #       print(dim(plot.data))
        if (savePlot == TRUE){
            # plot rank-based coloured bubbles
-            p = qplot(x.data, y.data, data=plot.data, size=gsSize,asp=1,
+            p = ggplot2::qplot(x.data, y.data, data=plot.data, size=gsSize,asp=1,
                    colour=rank,
                    xlab = Xlab, ylab = Ylab,
                    xlim=xlimits,
                    ylim=ylimits)
            # customize bubbles colour
-            p = p + scale_colour_gradient(guide="colourbar", low="#56B1F7",
+            p = p + ggplot2::scale_colour_gradient(guide="colourbar", low="#56B1F7",
    high="#000000",
                    limits=c(1,100), na.value="black", name="Rank")
            # customize bubble size
-            p = p + scale_size("Cardinality", range=c(2,20))
+            p = p + ggplot2::scale_size("Gene set size", range=c(2,20))
            if (is.null(format) || tolower(format) == "pdf"){
                pdf(paste0(file.name, ".rank.pdf"), width = 10, height = 7,
                        useDingbats = FALSE)
                # label the bubbles of the top 10 gene sets
-                print(p + geom_text(size=5, mapping=aes(x=x.data, y=y.data,
+                print(p + ggplot2::geom_text(size=5, mapping=ggplot2::aes(x=x.data, y=y.data,
                                label=id),
                                data=plot.data.sig,
                                colour=sig.cols, vjust=-1, hjust=1) )
@@ -68,7 +64,7 @@ generateSummaryPlots <- function(plot.data, file.name = "resumPlot", Xlab="-log1
            }
            if (is.null(format) || tolower(format) == "png"){
                png(paste0(file.name, ".rank.png"), width = 800, height = 700)
-                print(p + geom_text(size=5, mapping=aes(x=x.data, y=y.data,
+                print(p + ggplot2::geom_text(size=5, mapping=ggplot2::aes(x=x.data, y=y.data,
        label=id),
                                data=plot.data.sig,
        colour=sig.cols, vjust=-1, hjust=1) )
@@ -81,7 +77,7 @@ generateSummaryPlots <- function(plot.data, file.name = "resumPlot", Xlab="-log1
        # plot direction-based coloured bubbles
        top.10.ids = as.character(plot.data[plot.data[, "rank"] <= firstN,
 "id"])
-        sig.ids = setdiff(plot.data[rank(-plot.data[,"sig"], na.last =
+        sig.ids = ggplot2::setdiff(plot.data[rank(-plot.data[,"sig"], na.last =
 TRUE) <= 5, "id"], top.10.ids)
        sig.cols = c(rep("black", length(top.10.ids)), rep("blue",
 length(sig.ids)))
@@ -92,18 +88,18 @@ plot.data[, "id"]), ]
                xlab = Xlab, ylab = Ylab,
                xlim=xlimits,
                ylim=ylimits)
-        p = p + scale_colour_gradient(guide="colourbar", low="#56B1F7",
+        p = p + ggplot2::scale_colour_gradient(guide="colourbar", low="#56B1F7",
 high="#E35F5F",
                limits=c(-1,1), na.value="black",
 name="Regulation Direction") # low="#5FE377"
-        p = p + scale_size("significance", range=c(2,20))
+        p = p + ggplot2::scale_size("significance", range=c(2,20))
        if (savePlot == TRUE){
            if (is.null(format) || tolower(format) == "pdf"){
                pdf(paste0(file.name, "direction_.pdf"), width = 10, height = 7,
                        useDingbats = FALSE)
-                print(p + geom_text(size=5, mapping=aes(x=x.data, y=y.data,
+                print(p + ggplot2::geom_text(size=5, mapping=ggplot2::aes(x=x.data, y=y.data,
                                label=id),
                                data=plot.data.sig,
                                colour=sig.cols, vjust=-1, hjust=1) )
@@ -111,14 +107,14 @@ name="Regulation Direction") # low="#5FE377"
            }
            if (is.null(format) || tolower(format) == "png"){
                png(paste0(file.name, "_direction.png"), width = 800, height = 700)
-                print(p + geom_text(size=5, mapping=aes(x=x.data, y=y.data,
+                print(p + ggplot2::geom_text(size=5, mapping=aes(x=x.data, y=y.data,
                                label=id),
                                data=plot.data.sig,
                                colour=sig.cols, vjust=-1, hjust=1) )
                dev.off()
            }
        } else if (savePlot == FALSE && legend == TRUE){
-            print(p + geom_text(size=5, mapping=aes(x=x.data, y=y.data,
+            print(p + ggplot2::geom_text(size=5, mapping=aes(x=x.data, y=y.data,
                            label=id),
                            data=plot.data.sig,
                            colour=sig.cols, vjust=-1, hjust=1))

--- a/R/main.r
+++ b/R/main.r
-if(!require(parallel)){
-  install.packages("parallel")
-  require(parallel)
-}
 cGSEA = function(ElistObject,
    contrastMatrix,
    geneSetCollection = "H",
@@ -42,10 +37,10 @@ cGSEAcore = function(ElistObject, contrast, geneset = "H", specie = "Mus musculu
  if(camera == TRUE){
    print("RUNNING CAMERA")
    try({
-        camerastart = Sys.time()
+        camerastart = base::Sys.time()
        res$camera = runcamera(voom.results = ElistObject, contrast = contrast, genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        camerastop = Sys.time()
+        camerastop = base::Sys.time()
-        time$camera = difftime(camerastop, camerastart, units = c("secs"))
+        time$camera = base::difftime(camerastop, camerastart, units = c("secs"))
        })
@@ -54,100 +49,100 @@ cGSEAcore = function(ElistObject, contrast, geneset = "H", specie = "Mus musculu
  if (gage == TRUE){
    print("RUNNING GAGE")
    try({
-        gagestart = Sys.time()
+        gagestart = base::Sys.time()
        res$gage = rungage(voom.results = ElistObject, contrast = contrast, genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        gagestop = Sys.time()
+        gagestop = base::Sys.time()
-        time$gage = difftime(gagestop,gagestart, units = c("secs"))
+        time$gage = base::difftime(gagestop,gagestart, units = c("secs"))
        })
  }
  if (globaltest == TRUE){
    print("RUNNING GLOBALTEST")
    try({
-        globalteststart = Sys.time()
+        globalteststart = base::Sys.time()
        res$globaltest = runglobaltest(voom.results = ElistObject, contrast = contrast, genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        globalteststop = Sys.time()
+        globalteststop = base::Sys.time()
-        time$globaltest = difftime(globalteststop,globalteststart, units = c("secs"))
+        time$globaltest = base::difftime(globalteststop,globalteststart, units = c("secs"))
        })
  }
  if (gsva == TRUE){
    print("RUNNING GSVA")
    try({
-        gsvastart = Sys.time()
+        gsvastart = base::Sys.time()
        res$gsva = rungsva(method = "gsva", voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        gsvastop = Sys.time()
+        gsvastop = base::Sys.time()
-        time$gsva = difftime(gsvastop,gsvastart, units = c("secs"))
+        time$gsva = base::difftime(gsvastop,gsvastart, units = c("secs"))
        })
  }
  if (ssgsea == TRUE){
    print("RUNNING SSGSEA")
    try({
-        ssgseastart = Sys.time()
+        ssgseastart = base::Sys.time()
        res$ssgsea = rungsva(method = "ssgsea", voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        ssgseastop = Sys.time()
+        ssgseastop = base::Sys.time()
-        time$ssgsea = difftime(ssgseastop,ssgseastart, units = c("secs"))
+        time$ssgsea = base::difftime(ssgseastop,ssgseastart, units = c("secs"))
        })
  }
  if (zscore == TRUE){
    print("RUNNING ZSCORE")
    try({
-        zscorestart = Sys.time()
+        zscorestart = base::Sys.time()
        res$zscore = rungsva(method = "zscore", voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        zscorestop = Sys.time()
+        zscorestop = base::Sys.time()
-        time$zscore = difftime(zscorestop,zscorestart, units = c("secs"))
+        time$zscore = base::difftime(zscorestop,zscorestart, units = c("secs"))
        })
  }
  if (plage == TRUE){
    print("RUNNING PLAGE ")
    try({
-        plagestart = Sys.time()
+        plagestart = base::Sys.time()
        res$plage = rungsva(method = "plage", voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        plagestop = Sys.time()
+        plagestop = base::Sys.time()
-        time$plage = difftime(plagestop,plagestart, units = c("secs"))
+        time$plage = base::difftime(plagestop,plagestart, units = c("secs"))
        })
  }
  if (ora == TRUE){
    print("RUNNING ORA")
    try({
-        orastart = Sys.time()
+        orastart = base::Sys.time()
        res$ora = runora(voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        orastop = Sys.time()
+        orastop = base::Sys.time()
-        time$ora = difftime(orastop, orastart, units = c("secs"))
+        time$ora = base::difftime(orastop, orastart, units = c("secs"))
        })
  }
  if (padog == TRUE){
    print("RUNNING PADOG")
    try({
-        padogstart = Sys.time()
+        padogstart = base::Sys.time()
        res$padog = runpadog(voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        padogstop = Sys.time()
+        padogstop = base::Sys.time()
-        time$padog = difftime(padogstop, padogstart, units = c("secs"))
+        time$padog = base::difftime(padogstop, padogstart, units = c("secs"))
        })
  }
  if (roast == TRUE){
    print("RUNNING ROAST")
    try({
-        roaststart = Sys.time()
+        roaststart = base::Sys.time()
        res$roast = runroast(voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        roaststop = Sys.time()
+        roaststop = base::Sys.time()
-        time$roast = difftime(roaststop,roaststart, units = c("secs"))
+        time$roast = base::difftime(roaststop,roaststart, units = c("secs"))
        })
  }
  if (safe == TRUE){
    print("RUNNING SAFE")
    try({
-        safestart = Sys.time()
+        safestart = base::Sys.time()
        res$safe = runsafe(voom.results = ElistObject, contrast = contrast,  genesetIdx = genesetIdx, num.workers = num.workers, verbose = verbose)
-        safestop = Sys.time()
+        safestop = base::Sys.time()
-        time$safe = difftime(safestop, safestart, units = c("secs"))
+        time$safe = base::difftime(safestop, safestart, units = c("secs"))
        })
  }

--- a/R/plots.r
+++ b/R/plots.r
 runTimePlot = function(preparedData, savePlot = TRUE, directoryPath = directoryPath){
  print("Plotting runtime")
  pdf(paste(directoryPath,"cGSEA_methods_runtime.pdf",sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  barplot(unlist(preparedData$time), las = 2, ylab = "second", main = "RunTime")
+  graphics::barplot(unlist(preparedData$time), las = 2, ylab = "second", main = "RunTime")
  dev.off()
 }
@@ -9,85 +9,66 @@ runTimePlot = function(preparedData, savePlot = TRUE, directoryPath = directoryP
 clusteringPlot = function(preparedData, contrCondi, savePlot = TRUE, directoryPath = directoryPath){
  print(paste("Plotting cluster Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_clustering.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  plot(preparedData$clustering[[contrCondi]], hang = -0.1)
+  graphics::plot(preparedData$clustering[[contrCondi]], hang = -0.1)
  dev.off()
 }
 heatmapPlot = function(preparedData, contrCondi, savePlot = TRUE, directoryPath = directoryPath){
-    if(!require(pheatmap)){
-        install.packages("pheatmap")
-        require(pheatmap)
-    }
    if(savePlot == TRUE){
        print(paste("Plotting heatmap for condition", contrCondi))
        pdf(paste(directoryPath, contrCondi, "_pvalue_heatmap.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-        pheatmap(preparedData$heatmap[[contrCondi]], fontsize_row = 5)
+        pheatmap::pheatmap(preparedData$heatmap[[contrCondi]], fontsize_row = 5)
        dev.off()
    } else {
-        pheatmap(preparedData$heatmap[[contrCondi]])
+        pheatmap::pheatmap(preparedData$heatmap[[contrCondi]])
    }
 }
 pcaPlot = function(preparedData, contrCondi, savePlot = TRUE, directoryPath = directoryPath){
-  if(!require(FactoMineR)){
-    install.packages("FactoMineR")
-    require(FactoMineR)
-  }
-  res.pca = PCA(preparedData$PCA[[contrCondi]], scale.unit = T, axes = c(1,2), graph = FALSE)
+  res.pca = FactoMineR::PCA(preparedData$PCA[[contrCondi]], scale.unit = T, axes = c(1,2), graph = FALSE)
  eigen = res.pca$eig$`percentage of variance`
-  names(eigen) = rownames(res.pca$eig)
+  base::names(eigen) = base::rownames(res.pca$eig)
  print(paste("Plotting Eigen values Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_eigen_fall.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  barplot(eigen, las = 2, ylab = "%")
+  graphics::barplot(eigen, las = 2, ylab = "%")
  dev.off()
  print(paste("Plotting PCA Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_pca.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  plot(res.pca, choix = 'ind')
+  graphics::plot(res.pca, choix = 'ind')
  dev.off()
 }
 corPlot = function(preparedData, contrCondi, savePlot = TRUE, directoryPath = directoryPath){
-  if(!require(corrplot)){
-  install.packages("corrplot")
-  require(corrplot)
-  }
  print(paste("Plotting correlation Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_correlation.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  corrplot(preparedData$correlation[[contrCondi]], method = "circle", type = "full", order = "hclust",mar=c(10, 4, 4, 2) + 0.1)
+  corrplot::corrplot(preparedData$correlation[[contrCondi]], method = "circle", type = "full", order = "hclust",mar=c(10, 4, 4, 2) + 0.1)
  dev.off()
 }
 upsetrPlot = function(preparedData, contrCondi, savePlot = TRUE, directoryPath = directoryPath){
-  if(!require(UpSetR)){
-    install.packages("UpSetR")
-    require(UpSetR)
-  }
  print(paste("Plotting UpsetR Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_upsetr.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  upset(fromList(preparedData$snailPlot[[contrCondi]]),order.by = "freq")
+  UpSetR::upset(fromList(preparedData$snailPlot[[contrCondi]]),order.by = "freq")
  dev.off()
 }
 snailPlot = function(preparedData, contrCondi, savePlot = TRUE, min.intersection.size = 1, directoryPath = directoryPath){
-  if(!require(SuperExactTest)){
-    install.packages("SuperExactTest")
-    require(SuperExactTest)
-  }
-  snail = supertest(preparedData$snailPlot[[contrCondi]], n = length(levels(preparedData$snailPlot[[contrCondi]][[1]])))
+  snail = SuperExactTest::supertest(preparedData$snailPlot[[contrCondi]], n = length(levels(preparedData$snailPlot[[contrCondi]][[1]])))
  print(paste("Plotting Snail Plot for condition", contrCondi))
  if (savePlot == TRUE){
      pdf(paste(directoryPath, contrCondi, "_snailplot.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-      plot.msets(snail, sort.by="size", keep.empty.intersections=FALSE, min.intersection.size = min.intersection.size)
+      SuperExactTest::plot.msets(snail, sort.by="size", keep.empty.intersections=FALSE, min.intersection.size = min.intersection.size)
      dev.off()
  } else {
-      plot.msets(snail, sort.by="size", keep.empty.intersections=FALSE, min.intersection.size = min.intersection.size)
+      SuperExactTest::plot.msets(snail, sort.by="size", keep.empty.intersections=FALSE, min.intersection.size = min.intersection.size)
  }
 }
@@ -96,11 +77,11 @@ resumPlot1 = function(preparedData, contrCondi, savePlot = TRUE, directoryPath =
  print(paste("Plotting simple Summary Plot for condition", contrCondi))
  pdf(paste(directoryPath, contrCondi, "_simple_sumplot.pdf", sep = ""), width = 10, height = 7, useDingbats = FALSE)
-  plot(preparedData$resumPlot1[[contrCondi]][["x"]], preparedData$resumPlot1[[contrCondi]][["y"]], xlab = '-log10(pVal)', ylab = "avg logFC")
+  graphics::plot(preparedData$resumPlot1[[contrCondi]][["x"]], preparedData$resumPlot1[[contrCondi]][["y"]], xlab = '-log10(pVal)', ylab = "avg logFC")
-  abline(h = 0)
+  graphics::abline(h = 0)
-  abline(v = 0)
+  graphics::abline(v = 0)
-  text(preparedData$resumPlot1[[contrCondi]][['x']], preparedData$resumPlot1[[contrCondi]][['y']], labels=preparedData$resumPlot1[[contrCondi]][['labels']], cex= 0.7)
+  graphics::text(preparedData$resumPlot1[[contrCondi]][['x']], preparedData$resumPlot1[[contrCondi]][['y']], labels=preparedData$resumPlot1[[contrCondi]][['labels']], cex= 0.7)
-  abline(v = -log10(0.05), col = "red")
+  graphics::abline(v = -log10(0.05), col = "red")
  dev.off()
 }
@@ -111,9 +92,9 @@ resumPlot2 = function(preparedData, contrCondi, savePlot = TRUE, directoryPath=
 }
 cGSEAMakePlots = function(preparedData, directoryPath){
-  dir.create(file.path(directoryPath,"/plots/"), showWarnings = FALSE)
+  base::dir.create(base::file.path(directoryPath,"/plots/"), showWarnings = FALSE)
  for (condi in names(preparedData$result)){
-    dir.create(file.path(directoryPath,"/plots/", condi), showWarnings = FALSE)
+    base::dir.create(file.path(directoryPath,"/plots/", condi), showWarnings = FALSE)
    clusteringPlot(preparedData = preparedData, contrCondi = condi, directoryPath = paste0(directoryPath,"/plots/",condi,"/"))
    pcaPlot(preparedData = preparedData, contrCondi = condi, directoryPath = paste0(directoryPath,"/plots/",condi,"/"))
    corPlot(preparedData = preparedData, contrCondi = condi, directoryPath = paste0(directoryPath,"/plots/",condi,"/"))

--- a/R/runners.r
+++ b/R/runners.r
--- a/R/utils.r
+++ b/R/utils.r
--- a/README.md
+++ b/README.md
+# --- WARNING -----
+#### Permanantly moved to [GitHub](https://github.com/maxibor/coGSEA)
 # coGSEA
 **co**mparative **G**ene **S**et **E**nrichment **Analysis**
@@ -8,8 +13,7 @@ This version doesn't include the *SetRank* method for reason of package installa
 ## How to install
 ```
-require(devtools)
+devtools::install_github("maxibor/coGSEA")
-devtools::install_git('https://mborry@gitlab.pasteur.fr/mborry/coGSEA.git')
 ```
 ## Input files

--- a/preparing_microarray.Rmd
+++ b/preparing_microarray.Rmd
+---
+title: "Preparing micro-array data for coGSEA analysis"
+output: html_notebook
+---
+## Reading and preparing data for [coGSEA](link) analysis.
+In this document, we will show you how to prepare micro-array data to perform a coGSEA analysis.  
+This example dataset is from a [article](https://www.ncbi.nlm.nih.gov/pubmed/21836821) about Astrocymotas tumors published in 2011 by Liu et al.  
+You can download it on the Gene Expression Omnibus database with the accession number [GSE19728](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19728)  
+#
+#### Loading the necessary packages for this analysis
+```{r, eval}
+library(affy)
+library(hgu133plus2.db)
+library(edgeR)
+library(gage)
+library(coGSEA)
+```
+#### Loading the CEL FILES and Normalizing them usin `rma`
+```{r}
+setwd("~/GitLab/GSE19728/data/")
+celfiles = ReadAffy()
+celfiles = rma(celfiles)
+```
+#### Getting the expression data, and the Probe Names
+```{r}
+intensity = exprs(celfiles)
+intensity = cbind(rownames(intensity), intensity)
+colnames(intensity)[1] = "PROBEID"
+intensity = as.data.frame(intensity)
+intensity$PROBEID = as.character(intensity$PROBEID)
+```
+#### Getting the annotations of Probes IDS to ENTREZ accession numbers
+```{r}
+annots = select(hgu133plus2.db, intensity$PROBEID, "ENTREZID", "PROBEID")
+```
+#### Merge expression matrix and annotation
+```{r}
+res = merge(intensity, annots, by= "PROBEID")
+```
+#### Getting rid of PROBE ids and type casting 
+```{r}
+resmin = res[,2:ncol(res)]
+cname = colnames(resmin)
+resmin = apply(resmin, 2, as.numeric)
+colnames(resmin)= cname
+resmin = as.data.frame(resmin)
+```
+#### Aggregating the PROBES matching the same ENTREZ accession number by averaging them and applying a log transformation
+```{r}
+result = aggregate(. ~ ENTREZID, resmin, mean)
+result$ENTREZID = levels(as.factor(as.character(resmin$ENTREZID)))
+rownames(result) = result$ENTREZID
+result = result[,-1]
+result = log(result)
+```
+#### Visualizing in boxplot to check normalization
+```{r}
+boxplot(result, las = 2)
+```
+#### Changing column names to remove the `.CEL` filename extension
+```{r}
+colnames(result) = gsub(".CEL","", colnames(result))
+print(colnames(result))
+```
+#### Selecting only the samples we are interested in
+```{r}
+result2 = cbind(result$GSM492649_Astrocytomas_N, result$GSM525014, result$GSM525015, result$GSM525016, result$`GSM492662_Astrocytomas_T4-1`, result$`GSM492663_Astrocytomas_T4-2` , result$`GSM492664_Astrocytomas_T4-3`, result$`GSM492665_Astrocytomas_T4-4`, result$`GSM492666_Astrocytomas_T4-5`)
+colnames(result2) = c("GSM492649_Astrocytomas_N", "GSM525014", "GSM525015", "GSM525016","GSM492662_Astrocytomas_T4-1", "GSM492663_Astrocytomas_T4-2", "GSM492664_Astrocytomas_T4-3", "GSM492665_Astrocytomas_T4-4", "GSM492666_Astrocytomas_T4-5" )
+rownames(result2) = rownames(result)
+```
+#### Preparing the design matrix
+```{r}
+Normal = c(rep(1,4),rep(0,5))
+Tumor = c(rep(0,4),rep(1,5))
+design = cbind(Normal, Tumor)
+rownames(design) = colnames(result2)
+```
+#### Preparing the contrast matrix
+```{r}
+contr.matrix = makeContrasts(NormalVSTumor = Normal - Tumor, levels = design)
+```
+#### Preparing expression list object
+```{r}
+temp = new("EList")
+temp$design = design
+temp$E = as.matrix(result2)
+rownames(temp$E) = as.numeric(rownames(temp$E))
+temp$genes$ENTREZ = rownames(result2)
+temp$common.dispersion = estimateDisp(temp$E, design = temp$design)$common.dispersion
+temp$samples = colnames(result2)
+```
+#### Preparing gene set collection
+```{r}
+gs = gage::kegg.gsets(species = "hsa", id.type = "entrez")
+geneset = gs$kg.sets
+```
+#### Function for simplifying gene sets names
+```{r}
+nameshorter = function(names){
+  namemod = c()
+  for (i in seq(1,length(names))){
+    namemod[i] = paste(strsplit(names[i], " ")[[1]][-1], sep = "", collapse = " ")
+    namemod[i] = gsub("/","", names[i])
+    namemod[i] = gsub(" ","_", names[i])
+  }
+  return(namemod)
+}
+```
+#### Simplifying gene sets names
+```{r}
+names(geneset) = nameshorter(names(geneset))
+names(geneset) = gsub("/","_",names(geneset))
+```
+#### Saving necessary objects to RDS files
+```{r}
+saveRDS(contr.matrix, "~/GitLab/GSE19728/contrast.rds")
+saveRDS(temp, "~/GitLab/GSE19728/elist.rds")
+saveRDS(geneset, "~/GitLab/GSE19728/geneset.rds")
+```
+#### Reading necessary objects (generated above) from RDS files
+```{r}
+elist = readRDS("~/GitLab/GSE19728/elist.rds")
+contrast = readRDS("~/GitLab/GSE19728/contrast.rds")
+geneset = readRDS("~/GitLab/GSE19728/geneset.rds")
+```
+#### Running coGSEA analysis 
+```{r}
+coGSEA(ElistObject = elist, contrastMatrix = contrast, ENTREZGenesIds = elist$genes$ENTREZ, geneSetCollection = geneset,specie = "Homo sapiens", directoryPath = "~/GitLab/GSE19728/results", alpha = 0.05, pvalAdjMethod = "BH", pvalCombMethod = "sumlog",min.intersection.size = 1, GSEA.Methods = c("camera", "gage","globaltest", "gsva", "ssgsea", "zscore", "ora", "padog", "roast","safe"), num.workers = 4, shinyMode = FALSE)
+```
\ No newline at end of file
--- a/preparing_microarray.nb.html
+++ b/preparing_microarray.nb.html
--- a/preparing_rnaseq.Rmd
+++ b/preparing_rnaseq.Rmd
+---
+title: "Preparing micro-array data for coGSEA analysis"
+output: html_notebook
+---
+$$\log{FC} = \frac{log{E_{S4}}}{log{E_{normal}}}$$
+## Reading and preparing data for [coGSEA](link) analysis.
+In this document, we will show you how to prepare RNA-seq data to perform a coGSEA analysis.  
+You can download this example dataset on the Gene Expression Omnibus database with the accession number [GSE63310](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63310)    
+This dataset was analyzed in a very detailed [article](https://f1000research.com/articles/5-1408/v2) on how to do differential expression analysis that we strongly advise you to read.  
+The file you're looking for is : `GSE63310_RAW.tar`  
+#### Loading necessary packages
+```{r}
+library(edgeR)
+library(limma)
+library(Mus.musculus)
+library(coGSEA)
+```
+#### Reading files
+```{r}
+setwd("~/GitLab/GSE63310/data/")
+files <- c(
+"GSM1545535_10_6_5_11.txt",
+"GSM1545536_9_6_5_11.txt",   
+"GSM1545538_purep53.txt",
+"GSM1545539_JMS8-2.txt",
+"GSM1545540_JMS8-3.txt",
+"GSM1545541_JMS8-4.txt",
+"GSM1545542_JMS8-5.txt",
+"GSM1545544_JMS9-P7c.txt",
+"GSM1545545_JMS9-P8c.txt")
+x <- readDGE(files, columns=c(1,3))
+```
+#### Simplifying file names
+```{r}
+colnames(x) = substring(colnames(x), 12, nchar(colnames(x)))
+```
+#### Grouping by sample condition
+```{r}
+group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", "Basal", "ML", "LP"))
+x$samples$group <- group
+```
+#### Grouping by lane
+```{r}
+lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
+x$samples$lane <- lane
+x$samples
+```
+#### Annotation
+```{r}
+geneid <- rownames(x)
+genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),keytype="ENTREZID")
+dim(genes)
+head(genes)
+```
+#### Getting rid of duplicated annoations by keeping only the first one
+```{r}
+genes <- genes[!duplicated(genes$ENTREZID),]
+```
+#### Count per million of reads
+```{r}
+cpm <- cpm(x)
+```
+#### Removing genes lowly expressed  (genes with 0 expression across all samples)
+```{r}
+#Removing genes lowly expressed (genes with 0 expression across all samples)
+table(rowSums(x$counts==0)==9)
+keep.exprs <- rowSums(cpm>1)>=3
+x <- x[keep.exprs,, keep.lib.sizes=FALSE]
+dim(x)
+```
+#### Normlization with edgeR
+```{r}
+x <- calcNormFactors(x, method = "TMM")
+x$samples$norm.factors
+```
+#### Making the design matrix
+```{r}
+design <- model.matrix( ~ 0 + group + lane)
+colnames(design) <- gsub("group", "", colnames(design))
+rownames(design) = colnames(x)
+```
+#### Making the contrast matrix
+```{r}
+contr.matrix <- makeContrasts(
+  BasalvsLP = Basal-LP,
+  BasalvsML = Basal - ML,
+  LPvsML = LP - ML,
+  levels = colnames(design))
+contr.matrix
+```
+#### Applying voom transformation
+```{r}
+v <- voom(x, design, plot=F)
+v$genes$ENTREZ = rownames(v$E)
+```
+#### Saving objects
+```{r}
+saveRDS(v, "~/GitLab/GSE63310/elist.rds")
+saveRDS(contr.matrix, "~/GitLab/GSE63310/contrast.rds")
+```
+#### Reading RDS objects (previously genereated above)
+```{r}
+elist = readRDS("~/GitLab/GSE63310/elist.rds")
+contrast = readRDS("~/GitLab/GSE63310/contrast.rds")
+```
+#### Running coGSEA
+```{r}
+coGSEA(ElistObject = elist, contrastMatrix = contrast, ENTREZGenesIds = elist$genes$ENTREZ, geneSetCollection = "C2_KEGG",specie = "Mus musculus", directoryPath = "~/GitLab/GSE63310/results", alpha = 0.05, pvalAdjMethod = "BH", pvalCombMethod = "sumlog",min.intersection.size = 1, GSEA.Methods = c("camera", "gage","globaltest", "gsva", "ssgsea", "zscore", "ora", "padog", "roast","safe"), num.workers = 4, shinyMode = FALSE)
+```
--- a/preparing_rnaseq.nb.html
+++ b/preparing_rnaseq.nb.html
No results found