Mercurial > repos > amawla > edger

comparison edgeR.pl @ 4:a8a56766694e draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author amawla
date Mon, 24 Aug 2015 18:50:49 -0400
parents
children
comparison
equal deleted inserted replaced
3:3fb55f96f065 4:a8a56766694e
1 #!/bin/perl
2
3 #EdgeR.pl Version 0.0.3
4 #Contributors: Monica Britton, Blythe Durbin-Johnson, Joseph Fass, Nikhil Joshi, Alex Mawla
5
6 use strict;
7 use warnings;
8 use Getopt::Std;
9 use File::Basename;
10 use File::Path qw(make_path remove_tree);
11
12 $| = 1;
13
14 my %OPTIONS = (a => "glm", d => "tag", f => "BH", r => 5, u => "movingave");
15
16 getopts('a:d:e:f:h:lmn:o:r:tu:', \%OPTIONS);
17
18
19 die qq(
20 Usage: edgeR.pl [OPTIONS] factor::factor1::levels [factor::factor2::levels ...] cp::cont_pred1::values [cp::cont_pred2::values ...] cnt::contrast1 [cnt::contrast2] matrix
21
22 OPTIONS: -a STR Type Of Analysis [glm, pw, limma] (default: $OPTIONS{a})
23 -d STR The dispersion estimate to use for GLM analysis [tag] (default: $OPTIONS{d})
24 -e STR Path to place additional output files
25 -f STR False discovery rate adjustment method [BH] (default: $OPTIONS{f})
26 -h STR Name of html file for additional files
27 -l Output the normalised digital gene expression matrix in log2 format (only applicable when using limma and -n is also specified)
28 -m Perform all pairwise comparisons
29 -n STR File name to output the normalised digital gene expression matrix (only applicable when usinf glm or limma model)
30 -o STR File name to output csv file with results
31 -r INT Common Dispersion Rowsum Filter, ony applicable when 1 factor analysis selected (default: $OPTIONS{r})
32 -t Estimate Tagwise Disp when performing 1 factor analysis
33 -u STR Method for allowing the prior distribution for the dispersion to be abundance- dependent ["movingave"] (default: $OPTIONS{u})
34
35 ) if(!@ARGV);
36
37 my $matrix = pop @ARGV;
38
39 make_path($OPTIONS{e});
40 open(Rcmd,">$OPTIONS{e}/r_script.R") or die "Cannot open $OPTIONS{e}/r_script.R\n\n";
41 print Rcmd "
42 zz <- file(\"$OPTIONS{e}/r_script.err\", open=\"wt\")
43 sink(zz)
44 sink(zz, type=\"message\")
45
46 library(edgeR)
47 library(limma)
48
49 toc <- read.table(\"$matrix\", sep=\"\\t\", comment=\"\", as.is=T)
50 groups <- sapply(toc[1, -1], strsplit, \":\")
51 for(i in 1:length(groups)) { g <- make.names(groups[[i]][2]); names(groups)[i] <- g; groups[[i]] <- groups[[i]][-2] }
52 colnames(toc) <- make.names(toc[2,])
53 toc[,1] <- gsub(\",\", \".\", toc[,1])
54 tagnames <- toc[-(1:2), 1]
55 rownames(toc) <- toc[,1]
56 toc <- toc[-(1:2), -1]
57 for(i in colnames(toc)) toc[, i] <- as.numeric(toc[,i])
58 norm_factors <- calcNormFactors(as.matrix(toc))
59
60 pw_tests <- list()
61 uniq_groups <- unique(names(groups))
62 for(i in 1:(length(uniq_groups)-1)) for(j in (i+1):length(uniq_groups)) pw_tests[[length(pw_tests)+1]] <- c(uniq_groups[i], uniq_groups[j])
63 DGE <- DGEList(toc, lib.size=norm_factors*colSums(toc), group=names(groups))
64 pdf(\"$OPTIONS{e}/MA_plots_normalisation.pdf\", width=14)
65 for(i in 1:length(pw_tests)) {
66 j <- c(which(names(groups) == pw_tests[[i]][1])[1], which(names(groups) == pw_tests[[i]][2])[1])
67 par(mfrow = c(1, 2))
68 maPlot(toc[, j[1]], toc[, j[2]], normalize = TRUE, pch = 19, cex = 0.2, ylim = c(-10, 10), main=paste(\"MA Plot\", colnames(toc)[j[1]], \"vs\", colnames(toc)[j[2]]))
69 grid(col = \"blue\")
70 abline(h = log2(norm_factors[j[2]]), col = \"red\", lwd = 4)
71 maPlot(DGE\$counts[, j[1]]/DGE\$samples\$lib.size[j[1]], DGE\$counts[, j[2]]/DGE\$samples\$lib.size[j[2]], normalize = FALSE, pch = 19, cex = 0.2, ylim = c(-8, 8), main=paste(\"MA Plot\", colnames(toc)[j[1]], \"vs\", colnames(toc)[j[2]], \"Normalised\"))
72 grid(col = \"blue\")
73 }
74 dev.off()
75 pdf(file=\"$OPTIONS{e}/MDSplot.pdf\")
76 plotMDS(DGE, main=\"MDS Plot\", col=as.numeric(factor(names(groups)))+1, xlim=c(-3,3))
77 dev.off()
78 tested <- list()
79 ";
80
81 my $all_cont;
82 my @add_cont;
83 my @fact;
84 my @fact_names;
85 my @cp;
86 my @cp_names;
87 if(@ARGV) {
88 foreach my $input (@ARGV) {
89 my @tmp = split "::", $input;
90 if($tmp[0] eq "factor") {
91 $tmp[1] =~ s/[ \?\(\)\[\]\/\\=+<>:;\"\',\*\^\|\&-]/./g;
92 push @fact_names, $tmp[1];
93 $tmp[2] =~ s/:/\", \"/g;
94 $tmp[2] = "\"".$tmp[2]."\"";
95 push @fact, $tmp[2];
96 } elsif($tmp[0] eq "cp") {
97 $tmp[1] =~ s/[ \?\(\)\[\]\/\\=+<>:;\"\',\*\^\|\&-]/./g;
98 push @cp_names, $tmp[1];
99 $tmp[2] =~ s/:/, /g;
100 push @cp, $tmp[2];
101 } elsif($tmp[0] eq "cnt") {
102 push @add_cont, $tmp[1];
103 } else {
104 die("Unknown Input: $input\n");
105 }
106 }
107 }
108
109 if($OPTIONS{a} eq "pw") {
110 print Rcmd "
111 disp <- estimateCommonDisp(DGE, rowsum.filter=$OPTIONS{r})
112 ";
113 if(defined $OPTIONS{t}) {
114 print Rcmd "
115 disp <- estimateTrendedDisp (disp)
116 disp <- estimateTagwiseDisp(disp, trend=\"$OPTIONS{u}\")
117 pdf(file=\"$OPTIONS{e}/Tagwise_Dispersion_vs_Abundance.pdf\")
118 plotBCV(disp, cex=0.4)
119 abline(h=disp\$common.dispersion, col=\"firebrick\", lwd=3)
120 dev.off()
121 ";
122 }
123 print Rcmd "
124 for(i in 1:length(pw_tests)) {
125 tested[[i]] <- exactTest(disp, pair=pw_tests[[i]])
126 names(tested)[i] <- paste(pw_tests[[i]][2], \"-\", pw_tests[[i]][1], sep=\"\")
127 }
128 pdf(file=\"$OPTIONS{e}/Smear_Plots.pdf\")
129 for(i in 1:length(pw_tests)) {
130 dt <- decideTestsDGE(tested[[i]], p.value=0.05, adjust.method=\"$OPTIONS{f}\")
131 if(sum(dt) > 0) {
132 de_tags <- rownames(disp)[which(dt != 0)]
133 ttl <- \"Diff. Exp. Genes With adj. Pvalue < 0.05\"
134 } else {
135 de_tags <- rownames(topTags(tested[[i]], n=100)\$table)
136 ttl <- \"Top 100 tags\"
137 }
138
139 if(length(dt) < 5000) {
140 pointcex = 0.5
141 } else {
142 pointcex = 0.2
143 }
144 plotSmear(disp, pair=pw_tests[[i]], de.tags = de_tags, main = paste(\"Smear Plot\", names(tested)[i]), cex=0.5)
145 abline(h = c(-1, 1), col = \"blue\")
146 legend(\"topright\", c(\"2 Fold Change\", ttl) , lty=c(1, NA), pch=c(NA, 19), pt.cex=0.5, col=c(\"blue\", \"red\"), bty=\"n\")
147 }
148 dev.off()
149 ";
150 }
151 elsif($OPTIONS{a} eq "glm") {
152 for(my $fct = 0; $fct <= $#fact_names; $fct++) {
153 print Rcmd "
154 $fact_names[$fct] <- c($fact[$fct])
155 ";
156 }
157 for(my $fct = 0; $fct <= $#cp_names; $fct++) {
158 print Rcmd "
159 $cp_names[$fct] <- c($cp[$fct])
160 ";
161 }
162 my $all_fact = "";
163 if(@fact_names) {
164 foreach (@fact_names) {
165 $all_fact .= " + factor($_)";
166 }
167 }
168 my $all_cp = "";
169 if(@cp_names) {
170 $all_cp = " + ".join(" + ", @cp_names);
171 }
172 print Rcmd "
173 group_fact <- factor(names(groups))
174 design <- model.matrix(~ -1 + group_fact${all_fact}${all_cp})
175 colnames(design) <- sub(\"group_fact\", \"\", colnames(design))
176 ";
177 foreach my $fct (@fact_names) {
178 print Rcmd "
179 colnames(design) <- make.names(sub(\"factor.$fct.\", \"\", colnames(design)))
180 ";
181 }
182 if($OPTIONS{d} eq "tag") {
183 print Rcmd "
184 disp <- estimateGLMCommonDisp(DGE, design)
185 disp <- estimateGLMTrendedDisp(disp, design)
186 disp <- estimateGLMTagwiseDisp(disp, design)
187 fit <- glmFit(disp, design)
188 pdf(file=\"$OPTIONS{e}/Tagwise_Dispersion_vs_Abundance.pdf\")
189 plotBCV(disp, cex=0.4)
190 dev.off()
191 ";
192 }
193 if(@add_cont) {
194 $all_cont = "\"".join("\", \"", @add_cont)."\"";
195 print Rcmd "
196 cont <- c(${all_cont})
197 for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \"([^0-9])\", sep=\"\"), paste(i, \"\\\\1\", sep=\"\"), cont)
198 for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \"\$\", sep=\"\"), i, cont)
199 ";
200 } else {
201 print Rcmd "
202 cont <- NULL
203 ";
204 }
205 if(defined $OPTIONS{m}) {
206 print Rcmd "
207 for(i in 1:length(pw_tests)) cont <- c(cont, paste(pw_tests[[i]][2], \"-\", pw_tests[[i]][1], sep=\"\"))
208 ";
209 }
210 if(!defined $OPTIONS{m} && !@add_cont){
211 die("No Contrasts have been specified, you must at least either select multiple pairwise comparisons or specify a custom contrast\n");
212 }
213 print Rcmd "
214 fit <- glmFit(disp, design)
215 cont <- makeContrasts(contrasts=cont, levels=design)
216 for(i in colnames(cont)) tested[[i]] <- glmLRT(fit, contrast=cont[,i])
217 pdf(file=\"$OPTIONS{e}/Smear_Plots.pdf\")
218 for(i in colnames(cont)) {
219 dt <- decideTestsDGE(tested[[i]], p.value=0.05, adjust.method=\"$OPTIONS{f}\")
220 if(sum(dt) > 0) {
221 de_tags <- rownames(disp)[which(dt != 0)]
222 ttl <- \"Diff. Exp. Genes With adj. Pvalue < 0.05\"
223 } else {
224 de_tags <- rownames(topTags(tested[[i]], n=100)\$table)
225 ttl <- \"Top 100 tags\"
226 }
227
228 if(length(dt) < 5000) {
229 pointcex = 0.5
230 } else {
231 pointcex = 0.2
232 }
233 plotSmear(disp, de.tags = de_tags, main = paste(\"Smear Plot\", i), cex=pointcex)
234 abline(h = c(-1, 1), col = \"blue\")
235 legend(\"topright\", c(\"2 Fold Change\", ttl) , lty=c(1, NA), pch=c(NA, 19), pt.cex=0.5, col=c(\"blue\", \"red\"), bty=\"n\")
236 }
237 dev.off()
238
239 ";
240 if(defined $OPTIONS{n}) {
241 print Rcmd "
242 tab <- data.frame(ID=rownames(fit\$fitted.values), fit\$fitted.values, stringsAsFactors=F)
243 write.table(tab, \"$OPTIONS{n}\", quote=F, sep=\"\\t\", row.names=F)
244 ";
245 }
246 } elsif($OPTIONS{a} eq "limma") {
247 for(my $fct = 0; $fct <= $#fact_names; $fct++) {
248 print Rcmd "
249 $fact_names[$fct] <- c($fact[$fct])
250 ";
251 }
252 for(my $fct = 0; $fct <= $#cp_names; $fct++) {
253 print Rcmd "
254 $cp_names[$fct] <- c($cp[$fct])
255 ";
256 }
257 my $all_fact = "";
258 if(@fact_names) {
259 foreach (@fact_names) {
260 $all_fact .= " + factor($_)";
261 }
262 }
263 my $all_cp = "";
264 if(@cp_names) {
265 $all_cp = " + ".join(" + ", @cp_names);
266 }
267 print Rcmd "
268
269 group_fact <- factor(names(groups))
270 design <- model.matrix(~ -1 + group_fact${all_fact}${all_cp})
271 colnames(design) <- sub(\"group_fact\", \"\", colnames(design))
272 ";
273 foreach my $fct (@fact_names) {
274 print Rcmd "
275 colnames(design) <- make.names(sub(\"factor.$fct.\", \"\", colnames(design)))
276 ";
277 }
278 print Rcmd "
279 isexpr <- rowSums(cpm(toc)>1) >= 1
280 toc <- toc[isexpr, ]
281 pdf(file=\"$OPTIONS{e}/LIMMA_voom.pdf\")
282 y <- voom(toc, design, plot=TRUE, lib.size=colSums(toc)*norm_factors)
283 dev.off()
284
285 pdf(file=\"$OPTIONS{e}/LIMMA_MDS_plot.pdf\")
286 plotMDS(y, labels=colnames(toc), col=as.numeric(factor(names(groups)))+1, gene.selection=\"common\")
287 dev.off()
288 fit <- lmFit(y, design)
289 ";
290 if(defined $OPTIONS{n}) {
291 if(defined $OPTIONS{l}) {
292 print Rcmd "
293 tab <- data.frame(ID=rownames(y\$E), y\$E, stringsAsFactors=F)
294 ";
295 } else {
296 print Rcmd "
297 tab <- data.frame(ID=rownames(y\$E), 2^y\$E, stringsAsFactors=F)
298 ";
299 }
300 print Rcmd "
301 write.table(tab, \"$OPTIONS{n}\", quote=F, sep=\"\\t\", row.names=F)
302 ";
303 }
304 if(@add_cont) {
305 $all_cont = "\"".join("\", \"", @add_cont)."\"";
306 print Rcmd "
307 cont <- c(${all_cont})
308 for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \"([^0-9])\", sep=\"\"), paste(i, \"\\\\1\", sep=\"\"), cont)
309 for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \"\$\", sep=\"\"), i, cont)
310 ";
311 } else {
312 print Rcmd "
313 cont <- NULL
314 ";
315 }
316 if(defined $OPTIONS{m}) {
317 print Rcmd "
318 for(i in 1:length(pw_tests)) cont <- c(cont, paste(pw_tests[[i]][2], \"-\", pw_tests[[i]][1], sep=\"\"))
319 ";
320 }
321 if(!defined $OPTIONS{m} && !@add_cont){
322 die("No Contrasts have been specified, you must at least either select multiple pairwise comparisons or specify a custom contrast\n");
323 }
324 print Rcmd "
325 cont <- makeContrasts(contrasts=cont, levels=design)
326 fit2 <- contrasts.fit(fit, cont)
327 fit2 <- eBayes(fit2)
328 ";
329 } else {
330 die("Anaysis type $OPTIONS{a} not found\n");
331
332 }
333 if($OPTIONS{a} ne "limma") {
334 print Rcmd "
335 options(digits = 6)
336 tab <- NULL
337 for(i in names(tested)) {
338 tab_tmp <- topTags(tested[[i]], n=Inf, adjust.method=\"$OPTIONS{f}\")[[1]]
339 colnames(tab_tmp) <- paste(i, colnames(tab_tmp), sep=\":\")
340 tab_tmp <- tab_tmp[tagnames,]
341 if(is.null(tab)) {
342 tab <- tab_tmp
343 } else tab <- cbind(tab, tab_tmp)
344 }
345 tab <- cbind(Feature=rownames(tab), tab)
346 ";
347 } else {
348 print Rcmd "
349 tab <- NULL
350 options(digits = 6)
351 for(i in colnames(fit2)) {
352 tab_tmp <- topTable(fit2, coef=i, n=Inf, sort.by=\"none\", adjust.method=\"$OPTIONS{f}\")
353 colnames(tab_tmp)[-1] <- paste(i, colnames(tab_tmp)[-1], sep=\":\")
354 if(is.null(tab)) {
355 tab <- tab_tmp
356 } else tab <- cbind(tab, tab_tmp)
357 }
358 tab <- cbind(Feature=rownames(tab), tab)
359 ";
360 }
361 print Rcmd "
362 write.table(tab, \"$OPTIONS{o}\", quote=F, sep=\"\\t\", row.names=F)
363 sink(type=\"message\")
364 sink()
365 ";
366 close(Rcmd);
367 system("R --no-restore --no-save --no-readline < $OPTIONS{e}/r_script.R > $OPTIONS{e}/r_script.out");
368
369 open(HTML, ">$OPTIONS{h}");
370 print HTML "<html><head><title>EdgeR: Empirical analysis of digital gene expression data</title></head><body><h3>EdgeR Additional Files:</h3><p><ul>\n";
371 print HTML "<li><a href=MA_plots_normalisation.pdf>MA_plots_normalisation.pdf</a></li>\n";
372 print HTML "<li><a href=MDSplot.pdf>MDSplot.pdf</a></li>\n";
373 if($OPTIONS{a} eq "pw") {
374 if(defined $OPTIONS{t}) {
375 print HTML "<li><a href=Tagwise_Dispersion_vs_Abundance.pdf>Tagwise_Dispersion_vs_Abundance.pdf</a></li>\n";
376 }
377 print HTML "<li><a href=Smear_Plots.pdf>Smear_Plots.pdf</a></li>\n";
378 } elsif($OPTIONS{a} eq "glm" && $OPTIONS{d} eq "tag") {
379 print HTML "<li><a href=Tagwise_Dispersion_vs_Abundance.pdf>Tagwise_Dispersion_vs_Abundance.pdf</a></li>\n";
380 print HTML "<li><a href=Smear_Plots.pdf>Smear_Plots.pdf</a></li>\n";
381 } elsif($OPTIONS{a} eq "limma") {
382 print HTML "<li><a href=LIMMA_MDS_plot.pdf>LIMMA_MDS_plot.pdf</a></li>\n";
383 print HTML "<li><a href=LIMMA_voom.pdf>LIMMA_voom.pdf</a></li>\n";
384 }
385 print HTML "<li><a href=r_script.R>r_script.R</a></li>\n";
386 print HTML "<li><a href=r_script.out>r_script.out</a></li>\n";
387 print HTML "<li><a href=r_script.err>r_script.err</a></li>\n";
388 print HTML "</ul></p>\n";
389 close(HTML);
390