| Previous changeset 2:2e3f0fbd745a (2015-08-24) Next changeset 4:a8a56766694e (2015-08-24) |
|
Commit message:
Deleted selected files |
|
removed:
edgeR.pl |
| b |
| diff -r 2e3f0fbd745a -r 3fb55f96f065 edgeR.pl --- a/edgeR.pl Mon Aug 24 18:50:22 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,390 +0,0 @@\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' |