Merge pull request #31 from PF2-pasteur-fr/development

Development
PF2-pasteur-fr · May 2, 2017 · ffd6d18 · ffd6d18
2 parents 651598e + 9890059
commit ffd6d18
Show file tree

Hide file tree

Showing 47 changed files with 95 additions and 90 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: SARTools
 Type: Package
 Title: Statistical Analysis of RNA-Seq Tools
-Version: 1.4.0
-Date: 2016-11-29
+Version: 1.4.1
+Date: 2017-05-02
 Author: Marie-Agnes Dillies and Hugo Varet
 Maintainer: Hugo Varet <[email protected]>
 Depends: R (>= 3.3.0), DESeq2 (>= 1.12.0), edgeR (>= 3.12.0), xtable

diff --git a/NEWS b/NEWS
@@ -1,3 +1,8 @@
+CHANGES IN VERSION 1.4.1
+------------------------
+	o SARTools now accepts count files generated by featureCounts (still one count file per sample)
+	o added sample names and colors to the diagnostic of the size factors with DESeq2
+
 CHANGES IN VERSION 1.4.0
 ------------------------
 	o added a parameter in exportResults.DESeq2() and exportResults.edgeR() to avoid exporting the results

diff --git a/R/diagSizeFactorsPlots.R b/R/diagSizeFactorsPlots.R
@@ -3,12 +3,15 @@
 #' Plots to assess the estimations of the size factors
 #'
 #' @param dds a \code{DESeqDataSet} object
+#' @param group factor vector of the condition from which each sample belongs
+#' @param col colors for the plots
 #' @param outfile TRUE to export the figure in a png file
 #' @param plots vector of plots to generate
 #' @return Two files in the figures directory: diagSizeFactorsHist.png containing one histogram per sample and diagSizeFactorsTC.png for a plot of the size factors vs the total number of reads
 #' @author Marie-Agnes Dillies and Hugo Varet
 
-diagSizeFactorsPlots <- function(dds, outfile=TRUE, plots=c("diag","sf_libsize")){
+diagSizeFactorsPlots <- function(dds, group, col=c("lightblue","orange","MediumVioletRed","SpringGreen"), 
+                                 outfile=TRUE, plots=c("diag","sf_libsize")){
   # histograms
   if ("diag" %in% plots){
     ncol <- ifelse(ncol(counts(dds))<=4, ceiling(sqrt(ncol(counts(dds)))), 3)
@@ -27,13 +30,22 @@ diagSizeFactorsPlots <- function(dds, outfile=TRUE, plots=c("diag","sf_libsize")
     }
     if (outfile) dev.off()
   }
-
+  
   # total read counts vs size factors
   if ("sf_libsize" %in% plots){
-    if (outfile) png(filename="figures/diagSizeFactorsTC.png",width=1800,height=1800,res=300)  
-    plot(sizeFactors(dds), colSums(counts(dds)), pch=19, las=1, xlab="Size factors",
-	 ylab="Total number of reads",main="Diagnostic: size factors vs total number of reads")
-    abline(lm(colSums(counts(dds)) ~ sizeFactors(dds) + 0), lty=2, col="grey")
+    if (outfile) png(filename="figures/diagSizeFactorsTC.png", width=1800, height=1800, res=300)  
+    sf <- sizeFactors(dds)
+    libsize <- colSums(counts(dds))/1e6
+    plot(sf, libsize, pch=16, las=1,
+         col = col[as.integer(group)],
+         xlab="Size factors", ylab="Total number of reads (millions)",
+         main="Diagnostic: size factors vs total number of reads")
+    abs <- range(sf); meanAbs <- mean(abs); abs <- abs(abs[2]-abs[1])/25;
+    ord <- range(libsize); meanOrd <- mean(ord); ord <- abs(ord[2]-ord[1])/25;
+    text(sf - ifelse(sf > meanAbs, abs, -abs), 
+         libsize - ifelse(libsize > meanOrd, ord, -ord),
+         colnames(dds), col=col[as.integer(group)])
+    abline(lm(libsize ~ sf + 0), lty=2, col="grey")
     if (outfile) dev.off()
   }
 }

diff --git a/R/loadCountData.R b/R/loadCountData.R
@@ -4,32 +4,48 @@
 #'
 #' @param target target \code{data.frame} of the project returned by \code{loadTargetFile()}
 #' @param rawDir path to the directory containing the count files
-#' @param header a logical value indicating whether the file contains the names of the variables as its first line
 #' @param skip number of lines of the data file to skip before beginning to read data
 #' @param featuresToRemove vector of feature Ids (or character string common to feature Ids) to remove from the counts
 #' @return The \code{matrix} of raw counts with row names corresponding to the feature Ids and column names to the sample names as provided in the first column of the target.
 #' @details If \code{featuresToRemove} is equal to \code{"rRNA"}, all the features containing the character string "rRNA" will be removed from the counts.
 #' @author Marie-Agnes Dillies and Hugo Varet
 
-loadCountData <- function(target, rawDir="raw", header=FALSE, skip=0,
-                          featuresToRemove=c("alignment_not_unique", "ambiguous", "no_feature", "not_aligned", "too_low_aQual")){
-
+loadCountData <- function(target, rawDir="raw", skip=0, featuresToRemove=c("alignment_not_unique", "ambiguous", "no_feature", "not_aligned", "too_low_aQual")){
+
   labels <- as.character(target[,1])
   files <- as.character(target[,2])
-
-  rawCounts <- read.table(paste(rawDir,files[1],sep="/"), sep="\t", quote="\"", header=header, skip=skip)
-  rawCounts <- rawCounts[,1:2]
+
+  # detect if input count files are from featureCounts or HTSeq-count
+  f1 <- read.table(paste(rawDir,files[1],sep="/"), sep="\t", quote="\"", header=FALSE, nrows=1, stringsAsFactors=FALSE)
+  if (ncol(f1) >= 6 && all(apply(f1[1,1:6], 2, is.character))){
+    # counter featurecounts
+    idCol <- 1
+    countsCol <- 7
+    header <- TRUE
+  } else{
+    if (ncol(f1) >= 2 && is.character(f1[1,1]) & is.numeric(f1[1,2])){
+      # counter htseq-count
+      idCol <- 1
+      countsCol <- 2
+      header <- FALSE
+    } else{
+      stop("Can't determine if count files come from HTSeq-count or featureCounts")
+    }
+  }
+
+  rawCounts <- read.table(paste(rawDir,files[1],sep="/"), sep="\t", quote="\"", header=header, skip=skip, stringsAsFactors=FALSE)
+  rawCounts <- rawCounts[,c(idCol, countsCol)]
   colnames(rawCounts) <- c("Id", labels[1])
   if (any(duplicated(rawCounts$Id))) stop("Duplicated feature names in ", files[1])
   cat("Loading files:\n")
   cat(files[1],": ",length(rawCounts[,labels[1]])," rows and ",sum(rawCounts[,labels[1]]==0)," null count(s)\n",sep="")
-
+  
   for (i in 2:length(files)){
-  	tmp <- read.table(paste(rawDir,files[i],sep="/"), sep="\t", header=header, skip=skip)
-    tmp <- tmp[,1:2]
-  	colnames(tmp) <- c("Id", labels[i])
-  	if (any(duplicated(tmp$Id))) stop("Duplicated feature names in ", files[i])
-	rawCounts <- merge(rawCounts, tmp, by="Id", all=TRUE)
+    tmp <- read.table(paste(rawDir,files[i],sep="/"), sep="\t", header=header, skip=skip, stringsAsFactors=FALSE)
+    tmp <- tmp[,c(idCol, countsCol)]
+    colnames(tmp) <- c("Id", labels[i])
+    if (any(duplicated(tmp$Id))) stop("Duplicated feature names in ", files[i])
+    rawCounts <- merge(rawCounts, tmp, by="Id", all=TRUE)
     cat(files[i],": ",length(tmp[,labels[i]])," rows and ",sum(tmp[,labels[i]]==0)," null count(s)\n",sep="")
   }
 
@@ -45,11 +61,11 @@ loadCountData <- function(target, rawDir="raw", header=FALSE, skip=0,
   for (f in setdiff(featuresToRemove,"")){
     match <- grep(f, rownames(counts))
     if (length(match)>0){
-	  cat(rownames(counts)[match],sep="\n")
-	  counts <- counts[-match,]
-	}
+      cat(rownames(counts)[match],sep="\n")
+      counts <- counts[-match,]
+    }
   }
-
+  
   cat("\nTop of the counts matrix:\n")
   print(head(counts))
   cat("\nBottom of the counts matrix:\n")

diff --git a/R/loadTargetFile.R b/R/loadTargetFile.R
@@ -25,6 +25,7 @@ loadTargetFile <- function(targetFile, varInt, condRef, batch){
   if (any(is.na(cbind(target[,c(varInt, batch)], target[,1:2])))) stop("NA are present in the target file")
   # warning message if batch is numeric
   if (!is.null(batch) && is.numeric(target[,batch])) warning(paste("The", batch, "variable is numeric. Use factor() or rename the levels with letters to convert it into a factor"))
+  if (any(grepl("[[:punct:]]", as.character(target[,varInt])))) stop(paste("The", varInt, "variable contains punctuation characters, please remove them"))
   cat("Target file:\n")
   print(target)
   return(target)

diff --git a/R/summarizeResults.DESeq2.r b/R/summarizeResults.DESeq2.r
@@ -21,7 +21,7 @@ summarizeResults.DESeq2 <- function(out.DESeq2, group, independentFiltering=TRUE
   results <- out.DESeq2$results
 
   # diagnostic of the size factors
-  diagSizeFactorsPlots(dds=dds)
+  diagSizeFactorsPlots(dds=dds, group=group, col=col)
 
   # boxplots before and after normalisation
   countsBoxplots(dds, group=group, col=col)

diff --git a/README.md b/README.md
@@ -16,7 +16,7 @@ To install the SARTools package from GitHub, open a R session and:
 - install DESeq2, edgeR and genefilter with `source("http://bioconductor.org/biocLite.R")` and `biocLite(c("DESeq2", "edgeR", "genefilter"))` (if not installed yet)
 - install devtools with `install.packages("devtools")` (if not installed yet)
 - Note: Ubuntu users may have to install some libraries (libxml2-dev, libcurl4-openssl-dev and libssl-dev) to be able to install DESeq2 and devtools
-- for Windows users only, install [Rtools](http://cran.r-project.org/bin/windows/Rtools/) or check that it is already installed (needed to build the package)
+- for Windows users only, install [Rtools](https://cran.r-project.org/bin/windows/Rtools/) or check that it is already installed (needed to build the package)
 - load the devtools R package with `library(devtools)`
 - run `install_github("PF2-pasteur-fr/SARTools", build_vignettes=TRUE)`
 

diff --git a/man/BCVPlot.Rd b/man/BCVPlot.Rd
diff --git a/man/MAPlot.Rd b/man/MAPlot.Rd
diff --git a/man/MDSPlot.Rd b/man/MDSPlot.Rd
diff --git a/man/PCAPlot.Rd b/man/PCAPlot.Rd
diff --git a/man/SARTools-package.Rd b/man/SARTools-package.Rd
diff --git a/man/SERE.Rd b/man/SERE.Rd
diff --git a/man/barplotNull.Rd b/man/barplotNull.Rd
diff --git a/man/barplotTotal.Rd b/man/barplotTotal.Rd
diff --git a/man/checkParameters.DESeq2.Rd b/man/checkParameters.DESeq2.Rd
diff --git a/man/checkParameters.edgeR.Rd b/man/checkParameters.edgeR.Rd
diff --git a/man/clusterPlot.Rd b/man/clusterPlot.Rd
diff --git a/man/countsBoxplots.Rd b/man/countsBoxplots.Rd
diff --git a/man/densityPlot.Rd b/man/densityPlot.Rd
diff --git a/man/descriptionPlots.Rd b/man/descriptionPlots.Rd
diff --git a/man/diagSizeFactorsPlots.Rd b/man/diagSizeFactorsPlots.Rd
diff --git a/man/dispersionsPlot.Rd b/man/dispersionsPlot.Rd
diff --git a/man/exploreCounts.Rd b/man/exploreCounts.Rd
diff --git a/man/exportResults.DESeq2.Rd b/man/exportResults.DESeq2.Rd
diff --git a/man/exportResults.edgeR.Rd b/man/exportResults.edgeR.Rd
diff --git a/man/loadCountData.Rd b/man/loadCountData.Rd
diff --git a/man/loadTargetFile.Rd b/man/loadTargetFile.Rd
diff --git a/man/majSequences.Rd b/man/majSequences.Rd
diff --git a/man/nDiffTotal.Rd b/man/nDiffTotal.Rd
diff --git a/man/pairwiseScatterPlots.Rd b/man/pairwiseScatterPlots.Rd
diff --git a/man/rawpHist.Rd b/man/rawpHist.Rd
diff --git a/man/removeNull.Rd b/man/removeNull.Rd
diff --git a/man/run.DESeq2.Rd b/man/run.DESeq2.Rd
diff --git a/man/run.edgeR.Rd b/man/run.edgeR.Rd
diff --git a/man/summarizeResults.DESeq2.Rd b/man/summarizeResults.DESeq2.Rd
diff --git a/man/summarizeResults.edgeR.Rd b/man/summarizeResults.edgeR.Rd
diff --git a/man/tabIndepFiltering.Rd b/man/tabIndepFiltering.Rd
diff --git a/man/tabSERE.Rd b/man/tabSERE.Rd
diff --git a/man/volcanoPlot.Rd b/man/volcanoPlot.Rd
diff --git a/man/writeReport.DESeq2.Rd b/man/writeReport.DESeq2.Rd
diff --git a/man/writeReport.edgeR.Rd b/man/writeReport.edgeR.Rd