script_edgeR_1factor.r

In biomics-pasteur-fr/RNADiff: R Tools for RNA-Seq Projects at Institut Pasteur (Biomics Platform)

###################################################
### edgeR_1factor parameters: to be modified by the user
###################################################
rm(list=ls())                                        # remove all the objects of the R session

workspace <- "."           # workspace for the R session

projectName <- "project"                             # name of the project (cannot contain any ".")
analysisVersion <- "vN"                              # name of the analysis version (cannot contain any ".")

author <- "(Biomics Platform - Institut Pasteur)"      # author of the statistical report
researcher <- "Charles Friedel"                      # name of the researcher
chief <- "Louis Pasteur"                             # name of the head of unit

targetFile <- "target.txt"                           # path to the design/target file
infoFile <- NULL                                     # path to the annotation file (needed if 0 counts not in counts files)
rawDir <- "raw"                                      # path to the directory containing raw counts files
varInt <- "group"                                    # factor of interest
condRef <- "condRef"                                 # reference biological condition
batch <- NULL                                        # factor on which to adjust the statistical model: NULL (default) or "batch" for example

outfile <- TRUE                                      # TRUE to export figures, FALSE to display them in R
colors <- c("#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", # vector of colors of each group on the plots
            "#c2b280", "#848482", "#008856", "#e68fac", "#0067a5", 
            "#f99379", "#604e97", "#f6a600", "#b3446c", "#dcd300", 
            "#882d17", "#8db600", "#654522", "#e25822", "#2b3d26")

allComp <- TRUE                                      # make all the possible comparisons or only those to the reference level?
alpha <- 0.05                                        # threshold of statistical significance
adjMethod <- "BH"                                    # p-value adjustment method: "BH" (default) or "BY"

gene.selection <- "pairwise"                         # type de selection dans MDSplot

geneLengthFile <- NULL                               # path to the genes lenghts file (default is NULL)
featuresToRemove <- c("alignment_not_unique",        # names of the features to be removed (default is the HTSeq-count specific lines)
                      "ambiguous", "no_feature",
                      "not_aligned", "too_low_aQual") 

###################################################
### code chunk number 1: autres parametres et divers chargements
###################################################
setwd(workspace)
library(RNADiff)
library(knitr)

versionName <- paste(projectName, analysisVersion, sep="-")
ncol <- NULL                                         # largeur des tableaux dans le rapport

cat("Creation des dossiers d'exports\n") 
dir.create("figures", showWarnings=FALSE)
dir.create("tables", showWarnings=FALSE)

###################################################
### code chunk number 2: loadData
###################################################
cat("Chargement des annotations et longueurs des genes si besoin\n")
if (!is.null(infoFile)) print(head(info <- read.delim(infoFile, sep="\t", header=TRUE, stringsAsFactors=FALSE))) else info <- NULL
if (!is.null(geneLengthFile)) print(head(glength <- read.table(geneLengthFile, sep="\t", header=TRUE, stringsAsFactors=FALSE))) else glength <- NULL

cat("Chargement du target file\n")
print(target <- loadTargetFile(targetFile, varInt=varInt, condRef=condRef))
conds <- levels(target[,varInt])
group <- data.frame(group=factor(target[,varInt]))

cat("Chargement des donnees\n")
counts <- loadCountData(target, rawDir=rawDir, versionName=versionName, featuresToRemove=featuresToRemove)

cat("Verifier que les echantillons de counts sont dans le meme ordre que le target\n")
print(cbind(target=as.character(target[,1]),counts=colnames(counts)))

cat("Verifier que les identifiants dans info et glength sont les memes que dans les comptages\n")
checkInfoGlength(counts=counts, info=info, glength=glength)

####################################################
### code chunk number 3: description of raw data
####################################################
cat("\nFigure : nombre de reads par echantillon\n")
barplotTC(counts=counts, group=group, col=colors, out=outfile, versionName=versionName)

cat("Figure : nombre de comptages nuls par echantillon\n")
barplotNul(counts=counts, group=group, col=colors, out=outfile, versionName=versionName)
N <- nrow(counts) - nrow(removeNul(counts))
cat("\nNombre de genes avec que des comptages nuls :", N,"\n")

cat("\nFigure : estimation de la densite des comptages de chaque echantillon\n")
densityPlot(counts=counts, group=group, col=colors, out=outfile, versionName=versionName)

cat("\nFigure + tableau : sequences majoritaires pour chaque echantillon\n")
majSequences <- majSequences(counts=counts, group=group, versionName=versionName, col=colors, out=outfile)

cat("\nCalcul des SERE\n")
print(sere <- pairwiseSERE(counts, versionName=versionName))

cat("\nFigure : pairwise scatterplots of samples\n")
pairwiseScatterPlots(counts=counts, group=group, out=outfile, versionName=versionName)

####################################################
### code chunk number 4: creating DESeqDataSet object, normalization and estimateDispersion
####################################################
design <- formula(paste("~", ifelse(!is.null(batch), paste(batch,"+"), ""), varInt))
dge <- DGEList(counts=counts, remove.zeros=TRUE)
dge$design <- model.matrix(design, data=target)
print(design)

cat("Estimation des Effective Library Sizes\n")
dge <- calcNormFactors(dge)
tmm <- dge$samples$norm.factors
cat("TMM coefficients:\n")
print(tmm)

cat("\nCalcul des dispersions\n")
dge <- estimateGLMCommonDisp(dge, dge$design)
dge <- estimateGLMTrendedDisp(dge, dge$design)
dge <- estimateGLMTagwiseDisp(dge, dge$design)

cat("\nFigure : relation mean-dispersion\n")
BCVPlot(dge=dge, out=outfile, versionName=versionName)

####################################################
### code chunk number 5: Boxplot avant et apres normalisation
####################################################
cat("Figure : boxplots sur comptages bruts et normalises\n")
boxplotCounts(counts=dge$counts, group=group, col=colors, out=outfile, versionName=versionName)
boxplotCounts(counts=normCounts.edgeR(dge), group=group, col=colors, type="norm", out=outfile, versionName=versionName)

###################################################
### code chunk number 6: clustering + PCA of samples
###################################################
cat("Figure : dendrogramme de la classification sur comptages transformes\n")
clusterPlot(counts=cpm(dge, prior.count=2, log=TRUE), out=outfile, versionName=versionName)

cat("Figure : premier plan de Multi Dimensional Scaling sur les comptages transformes\n")
MDSPlot(dge=dge, group=group, gene.selection=gene.selection, col=colors, out=outfile, versionName=versionName)

####################################################
### code chunk number 7: analyse differentielle
####################################################
fit <- glmFit(dge, dge$design)

colsToTest <- grep(varInt, colnames(fit$design))
namesToTest <- paste0(gsub(varInt,"",colnames(fit$design)[colsToTest]),"_vs_",condRef)

res <- vector("list",length(colsToTest)); names(res) <- namesToTest;
for (i in 1:length(colsToTest)){
  cat(paste0("Comparison ", gsub("_", " ", namesToTest[i]),": testing coefficient ", colnames(fit$design)[colsToTest[i]]),"\n")
  lrt <- glmLRT(fit, coef=colsToTest[i])
  res[[namesToTest[i]]] <- topTags(lrt, n=nrow(dge$counts), adjust.method=adjMethod, sort.by="none")$table
}

# defining contrasts for the other comparisons (if applicable and wanted)
if (length(colsToTest)>=2 & allComp){
  print(colnames <- gsub(varInt,"",colnames(fit$design)))
  for (comp in combn(length(colsToTest),2,simplify=FALSE)){ 
    contrast <- numeric(ncol(dge$design))
    contrast[colsToTest[comp[1:2]]] <- c(-1,1)
    compname <- paste0(colnames[colsToTest[comp[2]]],"_vs_",colnames[colsToTest[comp[1]]])
    cat(paste0("Comparison ", gsub("_", "", compname),": testing contrast (", paste(contrast, collapse=", "),")"),"\n")
    lrt <- glmLRT(fit, contrast=contrast)
    res[[compname]] <- topTags(lrt, n=nrow(dge$counts), adjust.method=adjMethod, sort.by="none")$table
  }
}

###################################################
### code chunk number 8: export tables
###################################################
cat("Export des resultats\n")
complete <- exportComplete.edgeR(dge=dge, res=res, alpha=alpha, group=group[,1], conds=conds, versionName=versionName, info=info, export=TRUE)

cat("# genes up, down et total par comparaison\n")
print(nDiffTotal <- nDiffTotal(complete, alpha=alpha, versionName=versionName), quote=FALSE)

cat("Figure : nb de genes DE selon seuil FDR\n")
nbDiffSeuil(complete=complete, out=outfile, versionName=versionName)

if (!is.null(geneLengthFile)){
  cat("Export : comptages normalises par la longueur des genes\n")
  normGeneLength(counts=normCounts.edgeR(dge), glength=glength, versionName=versionName)
  geneLengthEffect(counts, complete, glength, out=outfile, versionName=versionName)
}

###################################################
### code chunk number 9: distribution of raw p-values and MA-plot
###################################################
cat("Figure : distribution des log2(Fold-Changes)\n")
diagLogFC(complete=complete, out=outfile, versionName=versionName)

cat("Figure : histogramme des p-valeurs brutes\n")
histoRawp(complete=complete, out=outfile, versionName=versionName)

cat("\nFigure : MA-plot\n")
MAplotDE(complete=complete, pvalCutoff=alpha, out=outfile, versionName=versionName)

cat("\nFigure : volcano-plot\n")
volcanoPlotDE(complete=complete, pvalCutoff=alpha, out=outfile, versionName=versionName)

cat("\nFigure : Venn diagram\n")
vennDiagramDE(complete=complete, alpha=alpha, out=outfile, versionName=versionName)

cat("\nFigure : heatmap\n")
heatmapDE(counts.trans=cpm(dge, prior.count=2, log=TRUE), complete=complete, 
          alpha=alpha, out=outfile, key.xlab="logCPM-centered data", versionName=versionName)

###################################################
### code chunk number 10: sessionInfo and saving
###################################################
cat("Sauvegarde des resultats\n")
sessionInfo <- sessionInfo()
pckVersionRNADiff <- packageVersion("RNADiff")
pckVersionedgeR <- packageVersion("edgeR")
save.image(file=paste0(versionName, ".RData"))