|
3
|
1 ################################################################################
|
|
|
2 ### R script to compare several conditions with the SARTools and edgeR packages
|
|
|
3 ### Hugo Varet
|
|
|
4 ### May 16th, 2018
|
|
|
5 ### designed to be executed with SARTools 1.6.3
|
|
|
6 ### run "Rscript template_script_edgeR_CL.r --help" to get some help
|
|
|
7 ################################################################################
|
|
|
8
|
|
|
9 rm(list=ls()) # remove all the objects from the R session
|
|
|
10 library(optparse) # to run the script in command lines
|
|
|
11
|
|
|
12 # options list with associated default value.
|
|
|
13 option_list <- list(
|
|
|
14 make_option(c("-P", "--projectName"),
|
|
|
15 default=basename(getwd()),
|
|
|
16 dest="projectName",
|
|
|
17 help="name of the project used for the report [default: name of the current directory]."),
|
|
|
18
|
|
|
19 make_option(c("-A", "--author"),
|
|
|
20 default=Sys.info()[7],
|
|
|
21 dest="author",
|
|
|
22 help="name of the report author [default: %default]."),
|
|
|
23
|
|
|
24 make_option(c("-t", "--targetFile"),
|
|
|
25 default="target.txt",
|
|
|
26 dest="targetFile",
|
|
|
27 help="path to the design/target file [default: %default]."),
|
|
|
28
|
|
|
29 make_option(c("-r", "--rawDir"),
|
|
|
30 default="raw",
|
|
|
31 dest="rawDir",
|
|
|
32 help="path to the directory containing the HTSeq files [default: %default]."),
|
|
|
33
|
|
|
34 make_option(c("-F", "--featuresToRemove"),
|
|
|
35 default="alignment_not_unique,ambiguous,no_feature,not_aligned,too_low_aQual",
|
|
|
36 dest="FTR",
|
|
|
37 help="names of the features to be removed, more than once can be specified [default: %default]"),
|
|
|
38
|
|
|
39 make_option(c("-v", "--varInt"),
|
|
|
40 default="group",
|
|
|
41 dest="varInt",
|
|
|
42 help="factor of interest [default: %default]"),
|
|
|
43
|
|
|
44 make_option(c("-c", "--condRef"),
|
|
|
45 default="WT",
|
|
|
46 dest="condRef",
|
|
|
47 help="reference biological condition [default: %default]"),
|
|
|
48
|
|
|
49 make_option(c("-b", "--batch"),
|
|
|
50 default=NULL,
|
|
|
51 dest="batch",
|
|
|
52 help="blocking factor [default: %default] or \"batch\" for example"),
|
|
|
53
|
|
|
54 make_option(c("-a", "--alpha"),
|
|
|
55 default=0.05,
|
|
|
56 dest="alpha",
|
|
|
57 help="threshold of statistical significance [default: %default]"),
|
|
|
58
|
|
|
59 make_option(c("-p", "--pAdjustMethod"),
|
|
|
60 default="BH",
|
|
|
61 dest="pAdjustMethod",
|
|
|
62 help="p-value adjustment method: \"BH\" or \"BY\" [default: %default]"),
|
|
|
63
|
|
|
64 make_option(c("-m", "--cpmCutoff"),
|
|
|
65 default=1,
|
|
|
66 dest="cpmCutoff",
|
|
|
67 help="counts-per-million cut-off to filter low counts"),
|
|
|
68
|
|
|
69 make_option(c("-g", "--gene.selection"),
|
|
|
70 default="pairwise",
|
|
|
71 dest="gene.selection",
|
|
|
72 help="selection of the features in MDSPlot [default: %default]"),
|
|
|
73
|
|
|
74 make_option(c("-n", "--normalizationMethod"),
|
|
|
75 default="TMM",
|
|
|
76 dest="normalizationMethod",
|
|
|
77 help="normalization method in calcNormFactors: \"TMM\", \"RLE\" or \"upperquartile\" [default: %default]"),
|
|
|
78
|
|
|
79 make_option(c("-C", "--colors"),
|
|
|
80 default="dodgerblue,firebrick1,MediumVioletRed,SpringGreen,chartreuse,cyan,darkorchid,darkorange",
|
|
|
81 dest="cols",
|
|
|
82 help="colors of each biological condition on the plots\n\t\t\"col1,col2,col3,col4\"\n\t\t[default: %default]"),
|
|
|
83
|
|
|
84 make_option(c("-f", "--forceCairoGraph"),
|
|
|
85 action="store_true",
|
|
|
86 default=FALSE,
|
|
|
87 dest="forceCairoGraph",
|
|
|
88 help="activate cairo type")
|
|
|
89 )
|
|
|
90
|
|
|
91 # now parse the command line to check which option is given and get associated values
|
|
|
92 parser <- OptionParser(usage="usage: %prog [options]",
|
|
|
93 option_list=option_list,
|
|
|
94 description="Compare two or more biological conditions in a RNA-Seq framework with edgeR.",
|
|
|
95 epilogue="For comments, bug reports etc... please contact Hugo Varet <hugo.varet@pasteur.fr>")
|
|
|
96 opt <- parse_args(parser, args=commandArgs(trailingOnly=TRUE), positional_arguments=0)$options
|
|
|
97
|
|
|
98 # get options and arguments
|
|
|
99 workDir <- getwd()
|
|
|
100 projectName <- opt$projectName # name of the project
|
|
|
101 author <- opt$author # author of the statistical analysis/report
|
|
|
102 targetFile <- opt$targetFile # path to the design/target file
|
|
|
103 rawDir <- opt$rawDir # path to the directory containing raw counts files
|
|
|
104 featuresToRemove <- unlist(strsplit(opt$FTR, ",")) # names of the features to be removed (specific HTSeq-count information and rRNA for example)
|
|
|
105 varInt <- opt$varInt # factor of interest
|
|
|
106 condRef <- opt$condRef # reference biological condition
|
|
|
107 batch <- opt$batch # blocking factor: NULL (default) or "batch" for example
|
|
|
108 alpha <- as.numeric(opt$alpha) # threshold of statistical significance
|
|
|
109 pAdjustMethod <- opt$pAdjustMethod # p-value adjustment method: "BH" (default) or "BY"
|
|
|
110 gene.selection <- opt$gene.selection # selection of the features in MDSPlot
|
|
|
111 normalizationMethod <- opt$normalizationMethod # normalization method in calcNormFactors
|
|
|
112 cpmCutoff <- opt$cpmCutoff # counts-per-million cut-off to filter low counts
|
|
|
113 colors <- unlist(strsplit(opt$cols, ",")) # vector of colors of each biologicial condition on the plots
|
|
|
114 forceCairoGraph <- opt$forceCairoGraph # force cairo as plotting device if enabled
|
|
|
115 # print(paste("workDir", workDir))
|
|
|
116 # print(paste("projectName", projectName))
|
|
|
117 # print(paste("author", author))
|
|
|
118 # print(paste("targetFile", targetFile))
|
|
|
119 # print(paste("rawDir", rawDir))
|
|
|
120 # print(paste("varInt", varInt))
|
|
|
121 # print(paste("condRef", condRef))
|
|
|
122 # print(paste("batch", batch))
|
|
|
123 # print(paste("alpha", alpha))
|
|
|
124 # print(paste("pAdjustMethod", pAdjustMethod))
|
|
|
125 # print(paste("featuresToRemove", featuresToRemove))
|
|
|
126 # print(paste("colors", colors))
|
|
|
127 # print(paste("gene.selection", gene.selection))
|
|
|
128 # print(paste("normalizationMethod", normalizationMethod))
|
|
|
129 # print(paste("cpmCutoff", cpmCutoff))
|
|
|
130
|
|
|
131 ################################################################################
|
|
|
132 ### running script ###
|
|
|
133 ################################################################################
|
|
|
134 # setwd(workDir)
|
|
|
135 library(SARTools)
|
|
|
136 if (forceCairoGraph) options(bitmapType="cairo")
|
|
|
137
|
|
|
138 # checking parameters
|
|
|
139 problem <- checkParameters.edgeR(projectName=projectName,author=author,targetFile=targetFile,
|
|
|
140 rawDir=rawDir,featuresToRemove=featuresToRemove,varInt=varInt,
|
|
|
141 condRef=condRef,batch=batch,alpha=alpha,pAdjustMethod=pAdjustMethod,
|
|
|
142 cpmCutoff=cpmCutoff,gene.selection=gene.selection,
|
|
|
143 normalizationMethod=normalizationMethod,colors=colors)
|
|
|
144 if (problem) quit(save="yes")
|
|
|
145
|
|
|
146 # loading target file
|
|
|
147 target <- loadTargetFile(targetFile=targetFile, varInt=varInt, condRef=condRef, batch=batch)
|
|
|
148
|
|
|
149 # loading counts
|
|
|
150 counts <- loadCountData(target=target, rawDir=rawDir, featuresToRemove=featuresToRemove)
|
|
|
151
|
|
|
152 # description plots
|
|
|
153 majSequences <- descriptionPlots(counts=counts, group=target[,varInt], col=colors)
|
|
|
154
|
|
|
155 # edgeR analysis
|
|
|
156 out.edgeR <- run.edgeR(counts=counts, target=target, varInt=varInt, condRef=condRef,
|
|
|
157 batch=batch, cpmCutoff=cpmCutoff, normalizationMethod=normalizationMethod,
|
|
|
158 pAdjustMethod=pAdjustMethod)
|
|
|
159
|
|
|
160 # MDS + clustering
|
|
|
161 exploreCounts(object=out.edgeR$dge, group=target[,varInt], gene.selection=gene.selection, col=colors)
|
|
|
162
|
|
|
163 # summary of the analysis (boxplots, dispersions, export table, nDiffTotal, histograms, MA plot)
|
|
|
164 summaryResults <- summarizeResults.edgeR(out.edgeR, group=target[,varInt], counts=counts, alpha=alpha, col=colors)
|
|
|
165
|
|
|
166 # save image of the R session
|
|
|
167 save.image(file=paste0(projectName, ".RData"))
|
|
|
168
|
|
|
169 # generating HTML report
|
|
|
170 writeReport.edgeR(target=target, counts=counts, out.edgeR=out.edgeR, summaryResults=summaryResults,
|
|
|
171 majSequences=majSequences, workDir=workDir, projectName=projectName, author=author,
|
|
|
172 targetFile=targetFile, rawDir=rawDir, featuresToRemove=featuresToRemove, varInt=varInt,
|
|
|
173 condRef=condRef, batch=batch, alpha=alpha, pAdjustMethod=pAdjustMethod, cpmCutoff=cpmCutoff,
|
|
|
174 colors=colors, gene.selection=gene.selection, normalizationMethod=normalizationMethod)
|
