annotate template_script_edgeR_CL.r @ 0:581d217c7337 draft

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