| Previous changeset 3:3fb55f96f065 (2015-08-24) |
|
Commit message:
Uploaded |
|
added:
edgeR.pl |
| b |
| diff -r 3fb55f96f065 -r a8a56766694e edgeR.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/edgeR.pl Mon Aug 24 18:50:49 2015 -0400 |
| [ |
| b'@@ -0,0 +1,390 @@\n+#!/bin/perl\n+\n+#EdgeR.pl Version 0.0.3\n+#Contributors: Monica Britton, Blythe Durbin-Johnson, Joseph Fass, Nikhil Joshi, Alex Mawla\n+\n+use strict;\n+use warnings;\n+use Getopt::Std;\n+use File::Basename;\n+use File::Path qw(make_path remove_tree);\n+\n+$| = 1;\n+\n+my %OPTIONS = (a => "glm", d => "tag", f => "BH", r => 5, u => "movingave");\n+\n+getopts(\'a:d:e:f:h:lmn:o:r:tu:\', \\%OPTIONS);\n+\n+ \n+die qq(\n+Usage: edgeR.pl [OPTIONS] factor::factor1::levels [factor::factor2::levels ...] cp::cont_pred1::values [cp::cont_pred2::values ...] cnt::contrast1 [cnt::contrast2] matrix\n+\n+OPTIONS:\t-a\tSTR\tType Of Analysis [glm, pw, limma] (default: $OPTIONS{a})\n+\t\t\t-d\tSTR\tThe dispersion estimate to use for GLM analysis [tag] (default: $OPTIONS{d})\n+\t\t\t-e\tSTR\tPath to place additional output files\n+\t\t\t-f\tSTR\tFalse discovery rate adjustment method [BH] (default: $OPTIONS{f})\n+\t\t\t-h\tSTR\tName of html file for additional files\n+\t\t\t-l\t\tOutput the normalised digital gene expression matrix in log2 format (only applicable when using limma and -n is also specified)\n+\t\t\t-m\t\tPerform all pairwise comparisons\n+\t\t\t-n\tSTR\tFile name to output the normalised digital gene expression matrix (only applicable when usinf glm or limma model)\n+\t\t\t-o\tSTR\tFile name to output csv file with results\n+\t\t\t-r\tINT\tCommon Dispersion Rowsum Filter, ony applicable when 1 factor analysis selected (default: $OPTIONS{r})\n+\t\t\t-t\t\tEstimate Tagwise Disp when performing 1 factor analysis\n+\t\t\t-u\tSTR\tMethod for allowing the prior distribution for the dispersion to be abundance- dependent ["movingave"] (default: $OPTIONS{u})\n+ \n+) if(!@ARGV);\n+\n+my $matrix = pop @ARGV;\n+\n+make_path($OPTIONS{e});\n+open(Rcmd,">$OPTIONS{e}/r_script.R") or die "Cannot open $OPTIONS{e}/r_script.R\\n\\n";\n+print Rcmd "\n+zz <- file(\\"$OPTIONS{e}/r_script.err\\", open=\\"wt\\")\n+sink(zz)\n+sink(zz, type=\\"message\\")\n+\n+library(edgeR)\n+library(limma)\n+\n+toc <- read.table(\\"$matrix\\", sep=\\"\\\\t\\", comment=\\"\\", as.is=T)\n+groups <- sapply(toc[1, -1], strsplit, \\":\\")\n+for(i in 1:length(groups)) { g <- make.names(groups[[i]][2]); names(groups)[i] <- g; groups[[i]] <- groups[[i]][-2] }\n+colnames(toc) <- make.names(toc[2,])\n+toc[,1] <- gsub(\\",\\", \\".\\", toc[,1])\n+tagnames <- toc[-(1:2), 1]\n+rownames(toc) <- toc[,1]\n+toc <- toc[-(1:2), -1]\n+for(i in colnames(toc)) toc[, i] <- as.numeric(toc[,i])\n+norm_factors <- calcNormFactors(as.matrix(toc))\n+\n+pw_tests <- list()\n+uniq_groups <- unique(names(groups))\n+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])\n+DGE <- DGEList(toc, lib.size=norm_factors*colSums(toc), group=names(groups))\n+pdf(\\"$OPTIONS{e}/MA_plots_normalisation.pdf\\", width=14)\n+for(i in 1:length(pw_tests)) {\n+\tj <- c(which(names(groups) == pw_tests[[i]][1])[1], which(names(groups) == pw_tests[[i]][2])[1])\n+\tpar(mfrow = c(1, 2))\n+\tmaPlot(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]]))\n+\tgrid(col = \\"blue\\")\n+\tabline(h = log2(norm_factors[j[2]]), col = \\"red\\", lwd = 4)\n+\tmaPlot(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\\"))\n+\tgrid(col = \\"blue\\")\n+}\n+dev.off()\n+pdf(file=\\"$OPTIONS{e}/MDSplot.pdf\\")\n+plotMDS(DGE, main=\\"MDS Plot\\", col=as.numeric(factor(names(groups)))+1, xlim=c(-3,3))\n+dev.off()\n+tested <- list()\n+";\n+\n+my $all_cont;\n+my @add_cont;\n+my @fact;\n+my @fact_names;\n+my @cp;\n+my @cp_names;\n+if(@ARGV) {\n+\tforeach my $input (@ARGV) {\n+\t\tmy @tmp = split "::", $input;\n+\t\tif($tmp[0] eq "factor") {\n+\t\t\t$tmp[1] =~ s/[ \\?\\(\\)\\[\\]\\/\\\\=+<>:;\\"\\\',\\*\\^\\|\\&-]/./g;\n+\t\t\tpush @fact_names, $tmp[1];\n+\t\t\t$tmp[2] =~ s/:/\\", \\"/g;\n+\t\t\t$tmp[2] = "\\"".$tmp[2]."\\"";\n+\t\t\tpush @fact, $tmp[2];\n+\t\t} elsif($tmp[0] eq "cp") {\n+\t'..b' $fct (@fact_names) {\n+\t\tprint Rcmd "\n+colnames(design) <- make.names(sub(\\"factor.$fct.\\", \\"\\", colnames(design)))\n+";\n+\t}\n+\tprint Rcmd "\n+isexpr <- rowSums(cpm(toc)>1) >= 1\n+toc <- toc[isexpr, ]\n+pdf(file=\\"$OPTIONS{e}/LIMMA_voom.pdf\\")\n+y <- voom(toc, design, plot=TRUE, lib.size=colSums(toc)*norm_factors)\n+dev.off()\n+\n+pdf(file=\\"$OPTIONS{e}/LIMMA_MDS_plot.pdf\\")\n+plotMDS(y, labels=colnames(toc), col=as.numeric(factor(names(groups)))+1, gene.selection=\\"common\\")\n+dev.off()\n+fit <- lmFit(y, design)\n+";\n+\tif(defined $OPTIONS{n}) {\n+\t\tif(defined $OPTIONS{l}) {\n+\t\t\tprint Rcmd "\n+tab <- data.frame(ID=rownames(y\\$E), y\\$E, stringsAsFactors=F)\n+";\n+\t\t} else {\n+\t\t\tprint Rcmd "\n+tab <- data.frame(ID=rownames(y\\$E), 2^y\\$E, stringsAsFactors=F)\n+";\n+\t\t}\n+\t\tprint Rcmd "\n+write.table(tab, \\"$OPTIONS{n}\\", quote=F, sep=\\"\\\\t\\", row.names=F)\n+";\n+\t}\n+\tif(@add_cont) {\n+\t\t$all_cont = "\\"".join("\\", \\"", @add_cont)."\\"";\n+\t\tprint Rcmd "\n+cont <- c(${all_cont})\n+for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \\"([^0-9])\\", sep=\\"\\"), paste(i, \\"\\\\\\\\1\\", sep=\\"\\"), cont)\n+for(i in uniq_groups) cont <- gsub(paste(groups[[i]], \\"\\$\\", sep=\\"\\"), i, cont)\n+";\n+\t} else {\n+\t\tprint Rcmd "\n+cont <- NULL\n+";\n+\t}\n+\tif(defined $OPTIONS{m}) {\n+\t\tprint Rcmd "\n+for(i in 1:length(pw_tests)) cont <- c(cont, paste(pw_tests[[i]][2], \\"-\\", pw_tests[[i]][1], sep=\\"\\"))\n+";\n+\t}\n+\tif(!defined $OPTIONS{m} && !@add_cont){\n+\t\tdie("No Contrasts have been specified, you must at least either select multiple pairwise comparisons or specify a custom contrast\\n");\n+\t}\n+\tprint Rcmd "\n+cont <- makeContrasts(contrasts=cont, levels=design)\n+fit2 <- contrasts.fit(fit, cont)\n+fit2 <- eBayes(fit2)\n+";\n+} else {\n+\tdie("Anaysis type $OPTIONS{a} not found\\n");\n+\n+}\n+if($OPTIONS{a} ne "limma") {\n+\tprint Rcmd "\n+options(digits = 6)\n+tab <- NULL\n+for(i in names(tested)) {\n+\ttab_tmp <- topTags(tested[[i]], n=Inf, adjust.method=\\"$OPTIONS{f}\\")[[1]]\n+\tcolnames(tab_tmp) <- paste(i, colnames(tab_tmp), sep=\\":\\")\n+\ttab_tmp <- tab_tmp[tagnames,]\n+\tif(is.null(tab)) {\n+\t\ttab <- tab_tmp\n+\t} else tab <- cbind(tab, tab_tmp)\n+}\n+tab <- cbind(Feature=rownames(tab), tab)\n+";\n+} else {\n+\tprint Rcmd "\n+tab <- NULL\n+options(digits = 6)\n+for(i in colnames(fit2)) {\n+\ttab_tmp <- topTable(fit2, coef=i, n=Inf, sort.by=\\"none\\", adjust.method=\\"$OPTIONS{f}\\")\n+\tcolnames(tab_tmp)[-1] <- paste(i, colnames(tab_tmp)[-1], sep=\\":\\")\n+\tif(is.null(tab)) {\n+\t\ttab <- tab_tmp\n+\t} else tab <- cbind(tab, tab_tmp)\n+}\n+tab <- cbind(Feature=rownames(tab), tab)\n+";\n+}\n+print Rcmd "\n+write.table(tab, \\"$OPTIONS{o}\\", quote=F, sep=\\"\\\\t\\", row.names=F)\n+sink(type=\\"message\\")\n+sink()\n+";\n+close(Rcmd);\n+system("R --no-restore --no-save --no-readline < $OPTIONS{e}/r_script.R > $OPTIONS{e}/r_script.out");\n+\n+open(HTML, ">$OPTIONS{h}");\n+print HTML "<html><head><title>EdgeR: Empirical analysis of digital gene expression data</title></head><body><h3>EdgeR Additional Files:</h3><p><ul>\\n";\n+print HTML "<li><a href=MA_plots_normalisation.pdf>MA_plots_normalisation.pdf</a></li>\\n";\n+print HTML "<li><a href=MDSplot.pdf>MDSplot.pdf</a></li>\\n";\n+if($OPTIONS{a} eq "pw") {\n+\tif(defined $OPTIONS{t}) {\n+\t\tprint HTML "<li><a href=Tagwise_Dispersion_vs_Abundance.pdf>Tagwise_Dispersion_vs_Abundance.pdf</a></li>\\n";\n+\t}\n+\tprint HTML "<li><a href=Smear_Plots.pdf>Smear_Plots.pdf</a></li>\\n";\n+} elsif($OPTIONS{a} eq "glm" && $OPTIONS{d} eq "tag") {\n+\tprint HTML "<li><a href=Tagwise_Dispersion_vs_Abundance.pdf>Tagwise_Dispersion_vs_Abundance.pdf</a></li>\\n";\n+\tprint HTML "<li><a href=Smear_Plots.pdf>Smear_Plots.pdf</a></li>\\n";\n+} elsif($OPTIONS{a} eq "limma") {\n+\tprint HTML "<li><a href=LIMMA_MDS_plot.pdf>LIMMA_MDS_plot.pdf</a></li>\\n";\n+\tprint HTML "<li><a href=LIMMA_voom.pdf>LIMMA_voom.pdf</a></li>\\n";\n+}\n+print HTML "<li><a href=r_script.R>r_script.R</a></li>\\n";\n+print HTML "<li><a href=r_script.out>r_script.out</a></li>\\n";\n+print HTML "<li><a href=r_script.err>r_script.err</a></li>\\n";\n+print HTML "</ul></p>\\n";\n+close(HTML);\n+\n' |