forked from PF2-pasteur-fr/SARTools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
template_script_edgeR.r
executable file
·82 lines (63 loc) · 4.58 KB
/
template_script_edgeR.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
################################################################################
### R script to compare several conditions with the SARTools and edgeR packages
### Hugo Varet
### March 20th, 2018
### designed to be executed with SARTools 1.6.3
################################################################################
################################################################################
### parameters: to be modified by the user ###
################################################################################
rm(list=ls()) # remove all the objects from the R session
workDir <- "C:/path/to/your/working/directory/" # working directory for the R session
projectName <- "projectName" # name of the project
author <- "Your name" # author of the statistical analysis/report
targetFile <- "target.txt" # path to the design/target file
rawDir <- "raw" # path to the directory containing raw counts files
featuresToRemove <- c("alignment_not_unique", # names of the features to be removed
"ambiguous", "no_feature", # (specific HTSeq-count information and rRNA for example)
"not_aligned", "too_low_aQual")# NULL if no feature to remove
varInt <- "group" # factor of interest
condRef <- "WT" # reference biological condition
batch <- NULL # blocking factor: NULL (default) or "batch" for example
alpha <- 0.05 # threshold of statistical significance
pAdjustMethod <- "BH" # p-value adjustment method: "BH" (default) or "BY"
cpmCutoff <- 1 # counts-per-million cut-off to filter low counts
gene.selection <- "pairwise" # selection of the features in MDSPlot
normalizationMethod <- "TMM" # normalization method: "TMM" (default), "RLE" (DESeq) or "upperquartile"
colors <- c("dodgerblue","firebrick1", # vector of colors of each biological condition on the plots
"MediumVioletRed","SpringGreen")
forceCairoGraph <- FALSE
################################################################################
### running script ###
################################################################################
setwd(workDir)
library(SARTools)
if (forceCairoGraph) options(bitmapType="cairo")
# checking parameters
checkParameters.edgeR(projectName=projectName,author=author,targetFile=targetFile,
rawDir=rawDir,featuresToRemove=featuresToRemove,varInt=varInt,
condRef=condRef,batch=batch,alpha=alpha,pAdjustMethod=pAdjustMethod,
cpmCutoff=cpmCutoff,gene.selection=gene.selection,
normalizationMethod=normalizationMethod,colors=colors)
# loading target file
target <- loadTargetFile(targetFile=targetFile, varInt=varInt, condRef=condRef, batch=batch)
# loading counts
counts <- loadCountData(target=target, rawDir=rawDir, featuresToRemove=featuresToRemove)
# description plots
majSequences <- descriptionPlots(counts=counts, group=target[,varInt], col=colors)
# edgeR analysis
out.edgeR <- run.edgeR(counts=counts, target=target, varInt=varInt, condRef=condRef,
batch=batch, cpmCutoff=cpmCutoff, normalizationMethod=normalizationMethod,
pAdjustMethod=pAdjustMethod)
# MDS + clustering
exploreCounts(object=out.edgeR$dge, group=target[,varInt], gene.selection=gene.selection, col=colors)
# summary of the analysis (boxplots, dispersions, export table, nDiffTotal, histograms, MA plot)
summaryResults <- summarizeResults.edgeR(out.edgeR, group=target[,varInt], counts=counts, alpha=alpha, col=colors)
# save image of the R session
save.image(file=paste0(projectName, ".RData"))
# generating HTML report
writeReport.edgeR(target=target, counts=counts, out.edgeR=out.edgeR, summaryResults=summaryResults,
majSequences=majSequences, workDir=workDir, projectName=projectName, author=author,
targetFile=targetFile, rawDir=rawDir, featuresToRemove=featuresToRemove, varInt=varInt,
condRef=condRef, batch=batch, alpha=alpha, pAdjustMethod=pAdjustMethod, cpmCutoff=cpmCutoff,
colors=colors, gene.selection=gene.selection, normalizationMethod=normalizationMethod)