diff --git a/DESCRIPTION b/DESCRIPTION
index cccded119582109549b196e352156cc4d0987709..c1791343536172517a4ff7aae82268912b5cb252 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,11 +2,11 @@ Package: ChIPuanaR
Type: Package
Title: Differential analysis of ChIP-Seq data
Version: 0.99.0
-Date: 2019-06-24
+Date: 2020-01-13
Author: Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
Maintainer: Maëlle Daunesse <maelle.daunesse@pasteur.fr>
-Depends: R (>= 3.4.0), DESeq2, limma, kableExtra, knitr
-Imports: stats, utils, graphics, grDevices, rmarkdown (>= 1.4), SummarizedExperiment, S4Vectors
+Depends: R (>= 3.4.0), DESeq2, limma, kableExtra, knitr, ggplot2
+Imports: stats, utils, graphics, ggrepel, grDevices, grid, ggdendro, gridExtra, rmarkdown (>= 1.4), SummarizedExperiment, S4Vectors, scales
biocViews: Software
VignetteBuilder: knitr
Encoding: latin1
@@ -15,4 +15,4 @@ Description: Provide R tools and an environment for the differential
perform statistical analysis/testing with DESeq2 or limma, export
results and create final report.
License: GPL-2
-RoxygenNote: 6.1.1
+RoxygenNote: 7.0.0
diff --git a/NAMESPACE b/NAMESPACE
index 4282e1341dd0ae7a9baf18ece1d13789268e29a1..e95a88d68f988e694e3cde883ce5edf8a31a16c9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,7 @@
exportPattern("^[a-zA-Z]")
import(DESeq2)
+import(ggplot2)
import(kableExtra)
import(knitr)
import(limma)
@@ -9,8 +10,11 @@ importFrom(S4Vectors,mcols)
importFrom(S4Vectors,metadata)
importFrom(SummarizedExperiment,assay)
importFrom(SummarizedExperiment,colData)
+importFrom(ggdendro,ggdendrogram)
+importFrom(ggrepel,geom_text_repel)
importFrom(grDevices,col2rgb)
importFrom(grDevices,dev.off)
+importFrom(grDevices,palette)
importFrom(grDevices,png)
importFrom(graphics,abline)
importFrom(graphics,barplot)
@@ -25,7 +29,16 @@ importFrom(graphics,plot)
importFrom(graphics,points)
importFrom(graphics,text)
importFrom(graphics,title)
+importFrom(grid,gpar)
+importFrom(grid,textGrob)
+importFrom(gridExtra,grid.arrange)
importFrom(rmarkdown,render)
+importFrom(scales,log10_trans)
+importFrom(scales,log_breaks)
+importFrom(scales,math_format)
+importFrom(scales,trans_breaks)
+importFrom(scales,trans_format)
+importFrom(scales,trans_new)
importFrom(stats,density)
importFrom(stats,dist)
importFrom(stats,dnorm)
@@ -45,5 +58,6 @@ importFrom(utils,head)
importFrom(utils,installed.packages)
importFrom(utils,packageVersion)
importFrom(utils,read.table)
+importFrom(utils,stack)
importFrom(utils,tail)
importFrom(utils,write.table)
diff --git a/R/MAPlot.R b/R/MAPlot.R
index 2f141535655173f9d2dbc25ebc684f44c876284d..ecef32f1fb0f1f37809d95f2564bd1f0e8b8a7a5 100755
--- a/R/MAPlot.R
+++ b/R/MAPlot.R
@@ -6,29 +6,42 @@
#' @param results A \code{list} of \code{data.frame} containing features results
#' @param alpha cut-off to apply on each adjusted p-value
#' @param outfile TRUE to export the figure in a png file
+#' @param log2FClim numeric vector containing both upper and lower y-axis limits for all the MA-plots produced (NULL by default to set them automatically)
#' @return One MA-plot per comparison
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
-MAPlot <- function(results, alpha=0.05, outfile=FALSE){
+MAPlot <- function(results, alpha=0.05, outfile=FALSE, log2FClim=NULL){
ncol <- ifelse(length(results)<=4, ceiling(sqrt(length(results))), 3)
nrow <- ceiling(length(results)/ncol)
if (outfile) png(filename="MAPlot.png", width=cairoSizeWrapper(1800*ncol), height=cairoSizeWrapper(1800*nrow), res=300)
- par(mfrow=c(nrow,ncol))
- for (name in names(results)){
- results.name <- results[[name]]
- results.name <- results.name[results.name$baseMean>0,]
- results.name$padj <- ifelse(is.na(results.name$padj),1,results.name$padj)
- results.name$log2FC <- results.name$log2FoldChange
- ylim <- 1.2 * c(-1,2) * quantile(abs(results.name$log2FC[is.finite(results.name$log2FC)]), probs=0.99)
- plot(results.name$baseMean, pmax(ylim[1], pmin(ylim[2], results.name$log2FC)),
- log = "x", cex=0.45, las = 1, ylim = ylim,
- col = ifelse(results.name[,"padj"] < alpha, ifelse(results.name[,"log2FC"] > 0,"red3","royalblue3"), "black"),
- pch = ifelse(results.name$log2FC<ylim[1], 6, ifelse(results.name$log2FC>ylim[2], 2, 20)),
- xlab = "Mean of normalized counts",
- ylab = expression(log[2]~fold~change),
- main = paste0("MA-plot - ",gsub("_"," ",name)))
- abline(h=0, lwd=1, col="lightgray")
- }
+ p <- list()
+ for (name in names(results)){
+ results.name <- results[[name]]
+ results.name <- results.name[which(results.name$baseMean>0),]
+ results.name$padj <- ifelse(is.na(results.name$padj), 1, results.name$padj)
+ results.name$DE <- factor(ifelse(results.name$padj <= alpha, "yes", "no"), levels=c("no", "yes"))
+ py <- results.name$log2FoldChange
+ if (is.null(log2FClim)) ymax <- quantile(abs(py[is.finite(py)]), probs=0.99) else ymax <- log2FClim
+ results.name$log2FoldChange[which(py > ymax)] <- ymax
+ results.name$log2FoldChange[which(py < -ymax)] <- -ymax
+ results.name$outfield <- factor(ifelse(py > ymax, "top", ifelse(py < -ymax, "bottom", "in")),
+ levels=c("bottom", "in", "top"))
+ p[[name]] <- ggplot(data=results.name,
+ aes(x=.data$baseMean, y=.data$log2FoldChange, color=.data$DE, fill=.data$DE, shape=.data$outfield)) +
+ scale_x_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ geom_point(show.legend=FALSE, alpha=0.5, size=0.8) +
+ scale_colour_manual(values=c("no"="black", "yes"="red"), drop=FALSE) +
+ scale_shape_manual(values=c("bottom"=25, "in"=21, "top"=24), drop=FALSE) +
+ scale_fill_manual(values=c("no"="black", "yes"="red"), drop=FALSE) +
+ scale_y_continuous(expand=expand_scale(mult=c(0.03, 0.03))) +
+ xlab("Mean of normalized counts") +
+ ylab(expression(log[2]~fold~change)) +
+ ggtitle(paste0("MA-plot - ", gsub("_"," ",name)))
+ }
+ tmpfun <- function(...) grid.arrange(..., nrow=nrow, ncol=ncol)
+ do.call(tmpfun, p)
if (outfile) dev.off()
}
diff --git a/R/NAMESPACE.R b/R/NAMESPACE.R
index 76ac1422cce184d733326df8682a536d6941bf92..400358fa1750308b265081ff61d71561581fa623 100755
--- a/R/NAMESPACE.R
+++ b/R/NAMESPACE.R
@@ -3,11 +3,17 @@
#' @import limma
#' @import kableExtra
#' @import knitr
-#' @importFrom grDevices col2rgb dev.off png
+#' @import ggplot2
+#' @importFrom grDevices col2rgb dev.off palette png
+#' @importFrom gridExtra grid.arrange
+#' @importFrom ggdendro ggdendrogram
+#' @importFrom ggrepel geom_text_repel
+#' @importFrom grid textGrob gpar
#' @importFrom graphics abline barplot boxplot curve hist legend lines pairs par plot points text title
#' @importFrom stats density dist dnorm formula hclust lm model.matrix p.adjust.methods prcomp quantile relevel sd var p.adjust
-#' @importFrom utils combn head installed.packages packageVersion read.table tail write.table
+#' @importFrom utils combn head installed.packages packageVersion read.table tail write.table stack
#' @importFrom rmarkdown render
+#' @importFrom scales trans_breaks log10_trans trans_new trans_format log_breaks math_format
#' @importFrom SummarizedExperiment assay colData
#' @importFrom S4Vectors mcols metadata
NULL
diff --git a/R/PCAPlot.R b/R/PCAPlot.R
index 8886a56d55ef0f69eda57bfcc11d182c3b477723..7cc9ea0c89ac326be40c14bb2e3e1a5ef7f8652d 100755
--- a/R/PCAPlot.R
+++ b/R/PCAPlot.R
@@ -12,40 +12,35 @@
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
PCAPlot <- function(counts.trans, conditions, n=min(500, nrow(counts.trans)),
- col=c("#8d8888","#7e9fab","#b59696","#9280a6","#2e2d2d","#5bbec9","#d48acc",palette()), outfile=FALSE){
+ col=c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482", "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"), outfile=FALSE){
# PCA on the 500 most variables features
rv = apply(counts.trans, 1, var, na.rm=TRUE)
pca = prcomp(t(counts.trans[order(rv, decreasing = TRUE), ][1:n,]))
prp <- pca$sdev^2 * 100 / sum(pca$sdev^2)
prp <- round(prp[1:3],2)
-
+
# create figure
- if (outfile) png(filename="PCA.png",width=cairoSizeWrapper(1800*2),height=cairoSizeWrapper(1800),res=300)
- par(mfrow=c(1,2))
- # axes 1 et 2
- abs=range(pca$x[,1]); abs=abs(abs[2]-abs[1])/25;
- ord=range(pca$x[,2]); ord=abs(ord[2]-ord[1])/25;
- plot(pca$x[,1], pca$x[,2],
- las = 1, cex = 2, pch = 16, col = col[as.integer(factor(conditions))],
- xlab = paste0("PC1 (",prp[1],"%)"),
- ylab = paste0("PC2 (",prp[2],"%)"),
- main = "PCA - Axes 1 and 2")
- abline(h=0,v=0,lty=2,col="lightgray")
- text(pca$x[,1] - ifelse(pca$x[,1]>0,abs,-abs), pca$x[,2] - ifelse(pca$x[,2]>0,ord,-ord),
- colnames(counts.trans), col=col[as.integer(factor(conditions))])
-
- # axes 1 et 3
- abs=range(pca$x[,1]); abs=abs(abs[2]-abs[1])/25;
- ord=range(pca$x[,3]); ord=abs(ord[2]-ord[1])/25;
- plot(pca$x[,1], pca$x[,3],
- las = 1, cex = 2, pch = 16, col = col[as.integer(factor(conditions))],
- xlab = paste0("PC1 (",prp[1],"%)"),
- ylab = paste0("PC3 (",prp[3],"%)"),
- main = "PCA - Axes 1 and 3")
- abline(h=0,v=0,lty=2,col="lightgray")
- text(pca$x[,1] - ifelse(pca$x[,1]>0,abs,-abs), pca$x[,3] - ifelse(pca$x[,3]>0,ord,-ord),
- colnames(counts.trans), col=col[as.integer(factor(conditions))])
+ if (outfile) png(filename="PCA.png", width=cairoSizeWrapper(1900*2), height=cairoSizeWrapper(1600), res=200)
+
+ tmpFunction <- function(axes=c(1, 2)){
+ index1 <- axes[1]
+ index2 <- axes[2]
+ d <- data.frame(x=pca$x[,index1], y=pca$x[,index2],
+ conditions=conditions, sample=factor(rownames(pca$x), levels=rownames(pca$x)))
+ ggplot(data=d, aes(x=.data$x, y=.data$y, color=conditions, label=sample)) +
+ geom_point(show.legend=TRUE, size=3) +
+ labs(color="") +
+ scale_colour_manual(values=col) +
+ geom_text_repel(show.legend=FALSE, size=5, point.padding=0.2) +
+ xlab(paste0("PC", index1, " (",prp[index1],"%)")) +
+ ylab(paste0("PC", index2, " (",prp[index2],"%)"))
+ }
+ p1 <- tmpFunction(c(1, 2))
+ p2 <- tmpFunction(c(1, 3))
+ grid.arrange(p1, p2, nrow=1, ncol=2,
+ top=textGrob("Principal Component Analysis", x=0.01, hjust=0, gp=gpar(fontsize=15)))
+
if (outfile) dev.off()
-
+
return(invisible(pca$x))
}
diff --git a/R/barplotTotal.R b/R/barplotTotal.R
index 95a2420a7995404ddac57fd65583da2b52d627ff..cd1cffb5cd9d059e0783195665177b634b0f7d2c 100755
--- a/R/barplotTotal.R
+++ b/R/barplotTotal.R
@@ -10,18 +10,19 @@
#' @return A bar plot of the total number of reads per sample
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
-barplotTotal <- function(counts, conditions, col=c("#8d8888","#7e9fab","#b59696","#9280a6","#2e2d2d","#5bbec9","#d48acc",palette()), outfile=FALSE){
+barplotTotal <- function(counts, conditions,
+ col=c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482", "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile=FALSE){
if (outfile) png(filename="BarplotTotal.png",width=min(3600,800+800*ncol(counts)),height=3500,res=300)
- libsize <- colSums(counts)/1e6
- par(mar=c(8,5,4,2))
- barplot(libsize,
- main = "Total read count per sample (million)",
- ylab = "Total read count (million)",
- ylim = c(0, max(libsize)*1.2),
- col = col[as.integer(as.factor(conditions))],
- las = 2,
- cex.names = 0.7,
- cex.main=2)
- legend("topright", levels(as.factor(conditions)), fill=col[1:nlevels(as.factor(conditions))], bty="n")
+ d <- data.frame(tc=colSums(counts)/1e6, sample=factor(colnames(counts), colnames(counts)), conditions)
+ print(ggplot(d, aes(x=.data$sample, y=.data$tc, fill=.data$conditions)) +
+ geom_bar(stat="identity", show.legend=TRUE) +
+ labs(fill="") +
+ scale_fill_manual(values=col) +
+ xlab("Samples") +
+ ylab("Total read count (million)") +
+ scale_y_continuous(expand=expand_scale(mult=c(0.01, 0.05))) +
+ ggtitle("Total read count per sample (million)") +
+ theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5)))
if (outfile) dev.off()
}
diff --git a/R/clusterPlot.R b/R/clusterPlot.R
index dc88e384359028c5a2c9480f37fe9a275665b8d7..5aadf6e031ccbb46e9f0fc88c95129de827f8fe5 100755
--- a/R/clusterPlot.R
+++ b/R/clusterPlot.R
@@ -12,9 +12,12 @@
clusterPlot <- function(counts.trans, conditions, outfile=FALSE){
hc <- hclust(dist(t(counts.trans)), method="ward.D")
if (outfile) png(filename="Cluster.png",width=1800,height=1800,res=300)
- plot(hc, hang=-1, ylab="Height",
- las=2,
- xlab="Method: Euclidean distance - Ward criterion",
- main="Cluster dendrogram")
+ print(ggdendrogram(hc, theme_dendro=FALSE) +
+ xlab("Samples") +
+ ylab("Height") +
+ ggtitle("Cluster dendrogram\nEuclidean distance, Ward criterion") +
+ theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
+ axis.text.y=element_text(angle=0)) +
+ scale_y_continuous(expand=expand_scale(mult=c(0.01, 0.05))))
if (outfile) dev.off()
}
diff --git a/R/countsBoxplots.R b/R/countsBoxplots.R
index 66161f8f25bcd7738a188aa1e24ff6ce1cdb109e..ab72af406dce157daa409021d4da99f75a81961f 100755
--- a/R/countsBoxplots.R
+++ b/R/countsBoxplots.R
@@ -11,7 +11,9 @@
#' @return A boxplots of the raw and normalized counts
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
-countsBoxplots <- function(results, conditions,method,col = c("#8d8888","#97adb5","#c9b6b6"), outfile=FALSE){
+countsBoxplots <- function(results, conditions, method,
+ col = c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482", "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile=FALSE){
# extract counts and norm counts
idx <- which(names(results[[1]]) == "baseMean")
tmp <- results[[1]][,1:(idx-1)]
@@ -19,19 +21,42 @@ countsBoxplots <- function(results, conditions,method,col = c("#8d8888","#97adb5
counts <- removeNull(tmp[,1:nbSamples])
norm.counts <- tmp[,(nbSamples+1):ncol(tmp)]
- if (outfile) png(filename="figures/countsBoxplots.png",width=2*min(2200,1800+800*ncol(norm.counts)/10),height=1800,res=300)
- par(mfrow=c(1,2))
+ if (outfile) png(filename="countsBoxplots.png",width=2*min(2200,1800+800*ncol(norm.counts)/10),height=1800,res=300)
# raw counts
- boxplot(log2(counts+1), col = col[as.integer(factor(conditions))], las = 2,
- main = "Raw counts distribution", ylab = expression(log[2] ~ (raw ~ count + 1)))
- legend("topright", levels(as.factor(conditions)), fill=col[1:nlevels(factor(conditions))], bty="n")
+ d <- stack(as.data.frame(counts))
+ d$conditions <- rep(conditions, each=nrow(counts))
+ p1 <- ggplot(d) +
+ geom_boxplot(aes(x=.data$ind, y=.data$values+1, fill=.data$conditions), show.legend=TRUE) +
+ labs(fill="") +
+ scale_y_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ scale_fill_manual(values=col) +
+ xlab("Samples") +
+ ylab("Raw counts") +
+ ggtitle("Raw counts distribution") +
+ theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
+
# norm counts
if(method=="DEseq2") {
norm.count <- log2(norm.counts+1)
}
+ norm.count <- removeNull(norm.counts)
+ d <- stack(as.data.frame(norm.counts))
+ d$conditions <- rep(conditions, each=nrow(norm.counts))
+ p2 <- ggplot(d) +
+ geom_boxplot(aes(x=.data$ind, y=.data$values+1, fill=.data$conditions), show.legend=TRUE) +
+ labs(fill="") +
+ scale_y_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ scale_fill_manual(values=col) +
+ xlab("Samples") +
+ ylab("Normalized counts") +
+ ggtitle("Normalized counts distribution") +
+ theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
+
+ grid.arrange(p1, p2, nrow=1, ncol=2)
- boxplot(norm.counts, col = col[as.integer(factor(conditions))], las = 2,
- main = "Normalized counts distribution", ylab = expression(log[2] ~ (norm ~ count + 1)))
- legend("topright", levels(as.factor(conditions)), fill=col[1:nlevels(factor(conditions))], bty="n")
if (outfile) dev.off()
}
\ No newline at end of file
diff --git a/R/densityPlot.R b/R/densityPlot.R
index e92e76c054c45354f672b8c9a9056b9bf3a684a8..e4ad71013fc0c5f8ddbd0786ca3769d6f4f5d28e 100755
--- a/R/densityPlot.R
+++ b/R/densityPlot.R
@@ -10,17 +10,22 @@
#' @return A density plot for each sample
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
-densityPlot <- function(counts, conditions, col=c("#8d8888","#7e9fab","#b59696","#9280a6","#2e2d2d","#5bbec9","#d48acc",palette()), outfile=TRUE){
+densityPlot <- function(counts, conditions,
+ col=c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482", "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile=TRUE){
if (outfile) png(filename="Densityplot_DESeq2.png",width=1800,height=1800,res=300)
counts <- removeNull(counts)
- plot(density(log2(counts[,1]+1)), las = 1, lwd = 2,
- main = "Density of counts distribution",
- xlab = expression(log[2] ~ (raw ~ count + 1)),
- ylim = c(0,max(apply(counts,2,function(x){max(density(log2(x+1))$y)}))*1.05),
- col = col[as.integer(conditions)[1]])
- for (i in 2:ncol(counts)){
- lines(density(log2(counts[,i]+1)),col=col[as.integer(conditions)[i]],lwd=2)
- }
- legend("topright", levels(conditions), lty=1, col=col[1:nlevels(conditions)], lwd=2, bty="n")
+ d <- stack(data.frame(counts))
+ d$conditions <- rep(conditions, each=nrow(counts))
+ print(ggplot(d, aes(x=.data$values+1)) +
+ stat_density(aes(group=.data$ind, color=.data$conditions), position="identity", geom="line", show.legend=TRUE) +
+ scale_x_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ labs(color="") +
+ scale_colour_manual(values=col) +
+ xlab("Raw counts") +
+ ylab("Density") +
+ ggtitle("Density of counts distribution"))
if (outfile) dev.off()
}
diff --git a/R/dispersionsPlot.R b/R/dispersionsPlot.R
index fb4217c56e45c25689fa487482c5dffc96371724..e1b37a0974dfdf41bc859431aab0233c4357bb0f 100755
--- a/R/dispersionsPlot.R
+++ b/R/dispersionsPlot.R
@@ -9,14 +9,30 @@
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
dispersionsPlot <- function(dds, outfile=FALSE){
- disp <- mcols(dds)$dispGeneEst
- disp <- disp[!is.na(disp)]
- disp <- disp[disp>1e-8]
- disp <- log(disp)
- mean.disp <- mean(disp,na.rm=TRUE)
- sd.disp <- sd(disp,na.rm=TRUE)
- if (outfile) png(filename="DispersionsPlot.png",width=1800,height=1800,res=300)
- # dispersions plot
- plotDispEsts(dds, main="Dispersions", las=1, xlab="Mean of normalized counts", ylab="Dispersion")
+ d <- as.data.frame(mcols(dds)[,c("baseMean", "dispGeneEst", "dispFit", "dispersion")])
+ d <- d[which(d$baseMean > 0),]
+ d <- data.frame(baseMean=rep(d$baseMean, 3),
+ value=c(d$dispGeneEst, d$dispersion, d$dispFit),
+ variable=factor(rep(c("dispGeneEst", "dispersion", "dispFit"), each=nrow(d)),
+ levels=c("dispGeneEst", "dispersion", "dispFit")))
+ if (outfile) png(filename="DispersionsPlot.png", width=1800, height=1800, res=300)
+ p <- ggplot(d, aes(x=.data$baseMean, y=.data$value, colour=.data$variable)) +
+ geom_point(size=0.1) +
+ scale_x_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format())) +
+ scale_y_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format())) +
+ ylab("Dispersion") +
+ xlab("Mean of normalized counts") +
+ scale_colour_manual(
+ values=c("Black", "#377eb8", "#e41a1c"),
+ breaks=c("dispGeneEst", "dispersion", "dispFit"),
+ labels=c("Estimate", "Final", "Fit"),
+ name="") +
+ guides(colour = guide_legend(override.aes = list(size=2))) +
+ ggtitle("Dispersions")
+ print(p)
if (outfile) dev.off()
}
diff --git a/R/meanVariancePlot.R b/R/meanVariancePlot.R
index b5f348f1f7f6d7f82dedfbefd3eb9cb44121f62f..2105c6794caa1b73ead96597a095ab1bc7cd3c33 100644
--- a/R/meanVariancePlot.R
+++ b/R/meanVariancePlot.R
@@ -16,5 +16,3 @@ meanVariancePlot <- function(voom, outfile=FALSE){
lines(voom$voom.line, col = "red")
if (outfile) dev.off()
}
-
-
diff --git a/R/pairwiseScatterPlots.R b/R/pairwiseScatterPlots.R
index 68bb698e949ac901b4af7329fe07e49bcb615842..f5cff420d36cdb1c1f035ed2a160a23620cf1761 100755
--- a/R/pairwiseScatterPlots.R
+++ b/R/pairwiseScatterPlots.R
@@ -4,33 +4,55 @@
#'
#' @name pairwiseScatterPlots
#' @param counts \code{matrix} of raw counts
-#' @param conditions factor vector of the condition from which each sample belongs
#' @param outfile TRUE to export the figure in a png file
#' @return A pairwise scatter plot with the SERE statistics in the lower panel
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
-pairwiseScatterPlots <- function(counts, conditions, outfile=FALSE){
+pairwiseScatterPlots <- function(counts, outfile=FALSE){
ncol <- ncol(counts)
- if (ncol <= 30){
- if (outfile) png(filename="PairwiseScatter.png", width=cairoSizeWrapper(700*ncol), height=cairoSizeWrapper(700*ncol), res=300)
- # defining panel and lower.panel functions
- panel <- function(x,y,...){points(x, y, pch=".");abline(a=0,b=1,lty=2);}
- lower.panel <- function(x,y,...){
- horizontal <- (par("usr")[1] + par("usr")[2]) / 2;
- vertical <- (par("usr")[3] + par("usr")[4]) / 2;
- text(horizontal, vertical, round(SERE(2^cbind(x,y) - 1), digits=2), cex=sqrt(ncol))
+ if (ncol <= 12){
+ if (outfile) png(filename="PairwiseScatter.png", width=cairoSizeWrapper(850*ncol), height=cairoSizeWrapper(700*ncol), res=200)
+ d <- data.frame(counts+1)
+ p <- list()
+ layout_matrix = matrix(NA, ncol=ncol, nrow=ncol)
+ k <- 1
+ for (i in 1:ncol){
+ for (j in 1:ncol){
+ if (i==j) next
+ if (i > j){
+ p[[k]] <- ggplot(data=cbind(d, z=1), aes_string(x=names(d)[i], y=names(d)[j], z="z")) +
+ stat_summary_2d(fun=function(z) log(sum(z)), bins=60, show.legend=FALSE) +
+ scale_x_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ scale_y_continuous(trans = log10_trans(),
+ breaks = trans_breaks("log10", function(x) 10^x),
+ labels = trans_format("log10", math_format(~10^.x))) +
+ xlab(names(d)[i]) + ylab(names(d)[j]) +
+ geom_abline(slope=1, intercept=0, linetype="dashed", col="lightgrey")
+ } else{
+ text = as.character(round(SERE(counts[,c(i,j)]), digits=2))
+ p[[k]] <- ggplot() + annotate("text", x=4, y=25, size=6, label=text) +
+ theme_bw() + xlab(names(d)[i]) + ylab(names(d)[j]) +
+ theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
+ axis.text.x=element_blank(), axis.ticks.x=element_blank(),
+ axis.text.y=element_blank(), axis.ticks.y=element_blank())
+ }
+ layout_matrix[j,i] <- k
+ k <- k+1
+ }
}
- # use of the paris function
- pairs(log2(counts+1), panel=panel, lower.panel=lower.panel,
- las=1, labels=paste(colnames(counts),conditions,sep="\n"),
- main="Pairwise scatter plot",cex.labels=sqrt(ncol),cex.main=ncol/4)
+ tmpfun <- function(...) grid.arrange(..., ncol=ncol, nrow=ncol,
+ layout_matrix=layout_matrix,
+ top=textGrob("Pairwise scatter plot", x=0.01, just="left", gp=gpar(fontsize=20)))
+ do.call(tmpfun, p)
if (outfile) dev.off()
} else{
- warning("No pairwise scatter-plot produced because of a too high number of samples (>30).")
+ warning("No pairwise scatter-plot produced because of a too high number of samples (>12).")
if (outfile) png(filename="figures/pairwiseScatter.png", width=1900, height=1900, res=300)
par(mar=c(1.5, 1.5, 1.5, 1.5))
plot(0, 0, bty="o", pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="")
- text(0, 0.3, "No pairwise scatter-plot produced because of\na too high number of samples (>30).", pos=1, cex=1.5)
+ text(0, 0.3, "No pairwise scatter-plot produced because of\na too high number of samples (>12).", pos=1, cex=1.5)
if (outfile) dev.off()
}
}
diff --git a/R/rawpHist.R b/R/rawpHist.R
index e252879aa5f07cd63573274d3a814578e84847c3..3433c4cead13a7aac6feddf1abf356877b7c9fb6 100755
--- a/R/rawpHist.R
+++ b/R/rawpHist.R
@@ -9,13 +9,20 @@
#' @author Marie-Agnes Dillies, Hugo Varet and Maëlle Daunesse
rawpHist <- function(result, outfile=FALSE){
- ncol <- ifelse(length(result)<=4, ceiling(sqrt(length(result))), 3)
+ ncol <- min(2, length(result))
nrow <- ceiling(length(result)/ncol)
if (outfile) png(filename="HistRawPvalue.png", width=cairoSizeWrapper(1800*ncol), height=cairoSizeWrapper(1800*nrow), res=300)
- par(mfrow=c(nrow,ncol))
- for (name in names(result)){
- hist(result[[name]][,"pvalue"], nclass=50, xlab="Raw p-value",
- col="skyblue", las=1, main=paste0("Distribution of raw p-values - ",gsub("_"," ",name)))
- }
+ p <- list()
+ for (name in names(result)){
+ result.name <- result[[name]]
+ p[[name]] <- ggplot(data=result.name, aes(x=.data$pvalue)) +
+ geom_histogram(binwidth=0.02) +
+ scale_y_continuous(expand=expand_scale(mult=c(0.01, 0.05))) +
+ xlab("Raw p-value") +
+ ylab("Frequency") +
+ ggtitle(paste0("Distribution of raw p-values - ", gsub("_"," ",name)))
+ }
+ tmpfun <- function(...) grid.arrange(..., nrow=nrow, ncol=ncol)
+ do.call(tmpfun, p)
if (outfile) dev.off()
}
diff --git a/inst/Report_ChIPuanaR.Rmd b/inst/Report_ChIPuanaR.Rmd
index e40b0f750d9314f1d0d7c1045b18d23a275e00b8..94b8bd2b9e14cd71c5f50a551527dabf357c4612 100644
--- a/inst/Report_ChIPuanaR.Rmd
+++ b/inst/Report_ChIPuanaR.Rmd
@@ -122,9 +122,9 @@ We may wish to assess the similarity between samples across conditions. A pairwi
- 1 for technical replicates (technical variability follows a Poisson distribution)
- greater than 1 for biological replicates and samples from different biological conditions (biological variability is higher than technical one, data are over-dispersed with respect to Poisson). The higher the SERE value, the lower the similarity. It is expected to be lower between biological replicates than between samples of different biological conditions. Hence, the SERE statistic can be used to detect inversions between samples.
-```{r pairewiseScatter, echo=FALSE,fig.align="center",fig.cap="Figure 2: Pairwise comparison of samples (not produced when more than 30 samples).", out.width="1200px"}
+```{r pairewiseScatter, echo=FALSE,fig.align="center",fig.cap="Figure 2: Pairwise comparison of samples (not produced when more than 12 samples).", out.width="1200px"}
#PairwiseScatter
-pairwiseScatterPlots(counts = counts, conditions = Conditions, outfile = FALSE)
+pairwiseScatterPlots(counts = counts, outfile = FALSE)
```
# 3. Variability within the experiment: data exploration
@@ -133,9 +133,9 @@ The main variability within the experiment is expected to come from biological d
Figure 3 sample clustering based on normalized data. An euclidean distance is computed between samples, and the dendrogram is built upon the Ward criterion. We expect this dendrogram to group replicates and separate biological conditions.
-```{r clusterplot, echo=FALSE,fig.align="center",fig.cap="Figure 3: Sample clustering based on normalized data.", out.width="600px"}
+```{r clusterplot, echo=FALSE,fig.align="center",fig.cap="Figure 3: Sample clustering based on normalized data.", out.width="600px", warning=FALSE, message=FALSE}
# Cluster plot
-clusterPlot(counts.trans = resAnDif$voom$E,conditions = Conditions,outfile = FALSE)
+clusterPlot(counts.trans = resAnDif$voom$E,conditions = Conditions, outfile = FALSE)
```
Another way of visualizing the experiment variability is to look at the first principal components of the PCA, as shown on the figure 4. On this figure, the first principal component (PC1) is expected to separate samples from the different biological conditions, meaning that the biological variability is the main source of variance in the data.
@@ -201,7 +201,7 @@ cat("The figure 6 shows the result of the dispersion estimation step. The x- and
if (method=="Limma") {
cat("A paragraph about Limma's method.\nLa variabilité systématique des données est modélisée avec une approche linéaire pour la différencier de la variabilité aléatoire. Cette modélisation linéaire est très similaire à l’ANOVA classique ou à la régression multiple, sauf qu’un modèle est adapté à chaque pic.\n")
-meanVariancePlot(voom = resAnDif$voom,outfile = FALSE)
+meanVariancePlot(voom = resAnDif$voom, outfile = FALSE)
cat("Explain figure 6")
}
diff --git a/inst/config.R b/inst/config.R
index a9f05a237398c73fb370dc70f76769db3f5be8b1..fb2b23c263d803db7e2a09a25c607f89953d86ff 100644
--- a/inst/config.R
+++ b/inst/config.R
@@ -1,6 +1,6 @@
-file <- "~/Documents/Pasteur/Projets/ChIPuana/chipuana_rtest/K4Me3_MOUSE_CuU_narrow_IDR0.05_Matrix_Optimal_Peak.mtx"
+file <- "K4Me3_MOUSE_CuU_Peak.mtx"
Conditions <- c("C","C","U","U")
-Replicates <- c("Rep1", "Rep2","Rep1", "Rep2")
+Replicates <- c("Rep1", "Rep2", "Rep1", "Rep2")
method <- "Limma" #Limma
normalisation <- "scale" #"scale"
spikes <- NULL
diff --git a/man/MAPlot.Rd b/man/MAPlot.Rd
index b29333ef3c8433558bc2cbc00ba1784cfc9e1b39..8fea48260fb38f538ad7567c93023703f81b54e5 100644
--- a/man/MAPlot.Rd
+++ b/man/MAPlot.Rd
@@ -4,7 +4,7 @@
\alias{MAPlot}
\title{MA-plots}
\usage{
-MAPlot(results, alpha = 0.05, outfile = FALSE)
+MAPlot(results, alpha = 0.05, outfile = FALSE, log2FClim = NULL)
}
\arguments{
\item{results}{A \code{list} of \code{data.frame} containing features results}
@@ -12,6 +12,8 @@ MAPlot(results, alpha = 0.05, outfile = FALSE)
\item{alpha}{cut-off to apply on each adjusted p-value}
\item{outfile}{TRUE to export the figure in a png file}
+
+\item{log2FClim}{numeric vector containing both upper and lower y-axis limits for all the MA-plots produced (NULL by default to set them automatically)}
}
\value{
One MA-plot per comparison
diff --git a/man/PCAPlot.Rd b/man/PCAPlot.Rd
index 790cc162e42892cf3014a68d0c1e60dedc791f7b..6afe3a855f52f81696a5a8e496cb7aadababb1c6 100644
--- a/man/PCAPlot.Rd
+++ b/man/PCAPlot.Rd
@@ -4,8 +4,14 @@
\alias{PCAPlot}
\title{PCA of samples (if use of DESeq2)}
\usage{
-PCAPlot(counts.trans, conditions, n = min(500, nrow(counts.trans)),
- col = c("#8d8888", "#97adb5", "#c9b6b6"), outfile = FALSE)
+PCAPlot(
+ counts.trans,
+ conditions,
+ n = min(500, nrow(counts.trans)),
+ col = c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482",
+ "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile = FALSE
+)
}
\arguments{
\item{counts.trans}{a matrix a transformed counts (VST- or rlog-counts)}
diff --git a/man/barplotTotal.Rd b/man/barplotTotal.Rd
index 241590e4217314e75b37ffca2352ce31c5968632..bdd37cf034769538c3d868547b83251ec6a07aaf 100644
--- a/man/barplotTotal.Rd
+++ b/man/barplotTotal.Rd
@@ -4,8 +4,13 @@
\alias{barplotTotal}
\title{Total number of reads per sample}
\usage{
-barplotTotal(counts, conditions, col = c("#8d8888", "#97adb5",
- "#c9b6b6"), outfile = FALSE)
+barplotTotal(
+ counts,
+ conditions,
+ col = c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482",
+ "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile = FALSE
+)
}
\arguments{
\item{counts}{\code{matrix} of counts}
diff --git a/man/countsBoxplots.Rd b/man/countsBoxplots.Rd
index 0ad800a48d4e3844ab7b0ba76a999edfd014d339..a64c94e05622fa987414d119145d741b75aae58d 100644
--- a/man/countsBoxplots.Rd
+++ b/man/countsBoxplots.Rd
@@ -4,8 +4,14 @@
\alias{countsBoxplots}
\title{Box-plots of (normalized) counts distribution per sample}
\usage{
-countsBoxplots(results, conditions, method, col = c("#8d8888", "#97adb5",
- "#c9b6b6"), outfile = FALSE)
+countsBoxplots(
+ results,
+ conditions,
+ method,
+ col = c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482",
+ "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile = FALSE
+)
}
\arguments{
\item{results}{results produced by DESeq2 or limma with ChIPuana}
diff --git a/man/densityPlot.Rd b/man/densityPlot.Rd
index 412857be51682b0331ea1d272660b5dbd4ddde8b..e2ddbd2b71dc7de83fa51c8e21e6b8d2c5ccbb67 100644
--- a/man/densityPlot.Rd
+++ b/man/densityPlot.Rd
@@ -4,8 +4,13 @@
\alias{densityPlot}
\title{Density plot of all samples}
\usage{
-densityPlot(counts, conditions, col = c("#8d8888", "#97adb5", "#c9b6b6"),
- outfile = TRUE)
+densityPlot(
+ counts,
+ conditions,
+ col = c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", "#c2b280", "#848482",
+ "#008856", "#e68fac", "#0067a5", "#f99379", "#604e97"),
+ outfile = TRUE
+)
}
\arguments{
\item{counts}{\code{matrix} of counts}
diff --git a/man/pairwiseScatterPlots.Rd b/man/pairwiseScatterPlots.Rd
index 1b7b31088606094be50e7e647b9301e4e7097a01..c7e7465453d7c9f7153e2487ed5e32326516ae8c 100644
--- a/man/pairwiseScatterPlots.Rd
+++ b/man/pairwiseScatterPlots.Rd
@@ -4,13 +4,11 @@
\alias{pairwiseScatterPlots}
\title{Scatter plots for pairwise comparaisons of log counts}
\usage{
-pairwiseScatterPlots(counts, conditions, outfile = FALSE)
+pairwiseScatterPlots(counts, outfile = FALSE)
}
\arguments{
\item{counts}{\code{matrix} of raw counts}
-\item{conditions}{factor vector of the condition from which each sample belongs}
-
\item{outfile}{TRUE to export the figure in a png file}
}
\value{
diff --git a/man/run.DESeq2.Rd b/man/run.DESeq2.Rd
index 664f2d0288834cf2b151bd1e33d18e75633ee4e2..6edef083a0c5c86997165a460b71af7f883a4608 100644
--- a/man/run.DESeq2.Rd
+++ b/man/run.DESeq2.Rd
@@ -4,8 +4,14 @@
\alias{run.DESeq2}
\title{Wrapper to run DESeq2}
\usage{
-run.DESeq2(counts, conditions, batch = NULL, pAdjustMethod = "BH",
- spikes = NULL, ...)
+run.DESeq2(
+ counts,
+ conditions,
+ batch = NULL,
+ pAdjustMethod = "BH",
+ spikes = NULL,
+ ...
+)
}
\arguments{
\item{counts}{\code{matrix} of raw counts}
diff --git a/man/run.Limma.Rd b/man/run.Limma.Rd
index 2e4bb35fe68dfeaca0180a03a2c360d63fd20a99..af5d6d2864cd0edac91518c49a8eb7af07928fb5 100644
--- a/man/run.Limma.Rd
+++ b/man/run.Limma.Rd
@@ -4,9 +4,16 @@
\alias{run.Limma}
\title{Wrapper to run Limma}
\usage{
-run.Limma(counts, conditions, normalize.method = "quantile",
- spikes = NULL, pAdjustMethod = "BH", batch = NULL,
- outfile = FALSE, ...)
+run.Limma(
+ counts,
+ conditions,
+ normalize.method = "quantile",
+ spikes = NULL,
+ pAdjustMethod = "BH",
+ batch = NULL,
+ outfile = FALSE,
+ ...
+)
}
\arguments{
\item{counts}{\code{matrix} of raw counts}