| Next changeset 1:4dd92ce61a5c (2017-12-08) |
|
Commit message:
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/ngsplot commit b'e9fcc157a7f2f2fa9d6ac9a58d425ff17c975f5c\n' |
|
added:
bin/install.db.tables.r bin/ngsplotdb.py bin/remove.db.tables.r bin/setTableDefaults.py database/bootstrap/dbfile.tbl database/bootstrap/default.tbl database/bootstrap/gn_list.txt database/dbfile.tbl database/default.tbl database/gn_list.txt lib/coverage.r lib/genedb.r lib/parse.args.r lib/plotlib.r ngs.plot.r ngsplot.xml replot.xml.notworking runNGSplot.pl runreplot.pl |
| b |
| diff -r 000000000000 -r 3ca58369469c bin/install.db.tables.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bin/install.db.tables.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| @@ -0,0 +1,73 @@ +#!/usr/bin/env Rscript +# Install the newly obtained DB tables into two ngs.plot system files: +# default.tbl and dbfile.tbl. + +args <- commandArgs(T) +# args <- c('/tmp/tmpLXAGXl', '/tmp/tmpUJduMe') +default.tbl <- args[1] +db.list <- args[2] + +# Genome-region default values. +anno.tbl <- read.table(default.tbl, header=TRUE, sep="\t", + blank.lines.skip=TRUE, stringsAsFactors=FALSE) +row.names(anno.tbl) <- paste(anno.tbl$Genome, anno.tbl$Region, sep=".") + +# Database table list. +anno.db.tbl <- read.table(db.list, + col.names=c("db.file", "Genome", "DB", "Region", + "FI.1", "FI.2", "FI.3"), + sep="\t", blank.lines.skip=TRUE, + stringsAsFactors=FALSE) +tss.tbl <- anno.db.tbl[anno.db.tbl$Region=="genebody", ] +tss.tbl$Region <- "tss" +tes.tbl <- anno.db.tbl[anno.db.tbl$Region=="genebody", ] +tes.tbl$Region <- "tes" +anno.db.tbl <- rbind(anno.db.tbl, tss.tbl, tes.tbl) +anno.db.tbl$URL <- NA # reserve for future use. + +# Extract further info columns. +default.fi <- anno.tbl[, c("DefaultFI1", "DefaultFI2", "DefaultFI3")] +db.fi <- anno.db.tbl[, c("FI.1", "FI.2", "FI.3")] +db.fi$tag <- paste(anno.db.tbl$Genome, anno.db.tbl$Region, sep=".") + +# DB scores. +getScore <- function(v.fi) { +# Calculate DB score based on intersection with default values. +# Args: +# v.fi: character vector of further infos. +# Returns: integer of DB score. + + default.vals <- default.fi[v.fi["tag"], ] + fi.vals <- v.fi[c("FI.1", "FI.2", "FI.3")] + 3 - length(intersect(default.vals, fi.vals)) +} +db.score <- apply(db.fi, 1, getScore) +anno.db.tbl$dbScore <- db.score + +# Save to text output. +progpath <- Sys.getenv('NGSPLOT') +if(progpath == "") { + stop("Set environment variable NGSPLOT before run the program. See README for details.\n") +} +def.f <- file.path(progpath, 'database', 'default.tbl') +db.f <- file.path(progpath, 'database', 'dbfile.tbl') +org.anno.tbl <- read.delim(def.f) +org.anno.db.tbl <- read.delim(db.f) +anno.tbl <- rbind(org.anno.tbl, anno.tbl) +anno.tbl <- unique(anno.tbl) +anno.db.tbl <- rbind(org.anno.db.tbl, anno.db.tbl) +anno.db.tbl <- unique(anno.db.tbl) +write.table(anno.tbl, file=def.f, row.names=F, sep="\t", quote=F) +write.table(anno.db.tbl, file=db.f, row.names=F, sep="\t", quote=F) + +# getScore <- function(i, db.fi, default.fi){ +# score <- 3 - length(intersect(as.matrix(default.fi[db.fi[i, "tag"], ]), as.matrix(db.fi[i, ]))) +# return(score) +# } + +# for (i in (1:dim(db.fi)[1])){ +# anno.db.tbl[i, "dbScore"] <- getScore(i, db.fi, default.fi) +# } + +# anno.version="0.01" +# save(anno.tbl, anno.db.tbl, anno.version, file="database.RData") |
| b |
| diff -r 000000000000 -r 3ca58369469c bin/ngsplotdb.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bin/ngsplotdb.py Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,670 @@\n+#!/usr/bin/env python\n+#\n+# ngsplotdb - PYTHON script to manipulate the ngs.plot database installation.\n+# Author: Li Shen, Mount Sinai School of Medicine\n+# Date created: Dec 28, 2012\n+# Last updated: May 28, 2013\n+#\n+# Functions implemented:\n+# \n+# list - List locally installed genomes.\n+# check - Check the integrity of local database.\n+# install - Install a genome from ngsplotdb.XXX.tar.gz.\n+# remove - Remove a locally installed genome.\n+\n+\n+# def dir_to_idspec(dirn):\n+# """Convert an FTP directory name to species ID and Name.\n+# This works with both UCSC and Ensembl FTP.\n+# """\n+# parts = dirn.split(\'_\')\n+# sid = parts[0][:3].lower() + parts[1][:3].lower()\n+# spec = str.title(parts[0] + \' \' + parts[1])\n+# return sid, spec\n+\n+# def listremote(database="refseq"):\n+# """List available genomes on a remote FTP site."""\n+\n+# import urllib2\n+\n+# url_ensembl = "ftp://ftp.ensembl.org/pub/current_gtf/"\n+# url_ucsc = "ftp://hgdownload.cse.ucsc.edu/goldenPath/currentGenomes/"\n+\n+# if database == "refseq":\n+# ftp = urllib2.urlopen(url_ucsc)\n+# elif database == "ensembl":\n+# ftp = urllib2.urlopen(url_ensembl)\n+# else:\n+# pass\n+\n+# # Read directory names from the remote FTP site.\n+# dir_names = ftp.read().splitlines()\n+# # Parse lines to obtain clean directory names only.\n+# if database == "refseq":\n+# dir_names = map(lambda x: x.split()[-3], dir_names)\n+# elif database == "ensembl":\n+# dir_names = map(lambda x: x.split()[-1], dir_names)\n+# else:\n+# pass\n+# # Filter directory names that do not correspond to a species.\n+# # On UCSC FTP, both "." and ".." will be retrieved from above and can \n+# # cause troubles.\n+\n+# id_spec_list = map(dir_to_idspec, dir_names)\n+# print "ID\\tName"\n+# for r in id_spec_list:\n+# print r[0], "\\t", r[1]\n+\n+def read_gnlist(root_path, mode):\n+ """Read installed genomes list.\n+\n+ Args:\n+ root_path: ngs.plot installation directory.\n+ mode: row reading mode - "vector" or "hash", correspond to each record\n+ read as a vector of hash table.\n+ Returns: (header split vector, hash table of installed genomes, \n+ vector of column widths)\n+ """\n+ from collections import defaultdict\n+\n+ gn_list = root_path + "/database/gn_list.txt"\n+ try:\n+ db_f = open(gn_list, "r")\n+ except IOError:\n+ print "Open {0} error: your ngs.plot database may be corrupted.".\\\n+ format(gn_list)\n+ sys.exit()\n+\n+ default_list = root_path + "/database/default.tbl"\n+ try:\n+ default_f = open(default_list, "r")\n+ except IOError:\n+ print "Open {0} error: your ngs.plot database may be corrupted.".\\\n+ format(default_list)\n+ sys.exit()\n+ default_l = []\n+ header = default_f.readline() # skip the header\n+ for rec in default_f:\n+ r_sp = rec.rstrip().split("\\t")\n+ default_l.append((r_sp[0], r_sp[2])) # tuple of genome and region, like ("dm3", "genebody")\n+ genome_region_dict = defaultdict(list)\n+ for genome,region in default_l:\n+ genome_region_dict[genome].append(region)\n+\n+ header = db_f.readline()\n+ h_sp = header.rstrip().split("\\t")\n+ h_sp.append("InstalledFeatures")\n+ v_cw = map(len, h_sp) # column widths initialize to header widths.\n+\n+ g_tbl = {}\n+ for rec in db_f:\n+ r_sp = rec.rstrip().split("\\t")\n+ r_sp.append(",".join(genome_region_dict[r_sp[0]]))\n+ r_cw = map(len, r_sp)\n+ v_cw = map(max, v_cw, r_cw)\n+\n+ r_tbl = dict(zip(h_sp, r_sp))\n+ if(mode == "vector"):\n+ g_tbl[r_tbl["ID"]] = r_sp\n+ elif(mode == "hash"):\n+ g_tbl[r_tbl["ID"]] = r_tbl\n+ else:\n+ pass\n+\n+ db_f.close()\n+\n+ return (h_sp, g_tbl, v_cw)\n+\n+\n+def update_gnlist(root_path, g_tbl, gn, assembly, species, ens_v, np_v):\n+ '..b'r)\n+ print "Done"\n+ else:\n+ print "Cannot find the genome in database. Nothing was done."\n+\n+\n+def rm_dbtbl(root_path, gn, sub_ftr=None):\n+ """Remove a genome from database meta tables.\n+\n+ Args:\n+ root_path: ngs.plot installation directory.\n+ gn: genome name to remove.\n+ Returns: None.\n+ """\n+\n+ import subprocess\n+\n+ if sub_ftr is None:\n+ subprocess.call(["remove.db.tables.r", gn])\n+ else:\n+ subprocess.call(["remove.db.tables.r", gn, sub_ftr])\n+\n+def chrnames(root_path, args):\n+ """List chromosome names for a given genome.\n+\n+ Args:\n+ root_path: ngs.plot installation directory.\n+ args.gn: genome name to remove.\n+ args.db: database(ensembl or refseq).\n+ Returns: None.\n+ """\n+\n+ import sys\n+\n+ db = args.db\n+ gn = args.gn\n+\n+ if db != "ensembl" and db != "refseq":\n+ print "Unrecognized database name: database must be ensembl or refseq."\n+ sys.exit()\n+\n+ chrnames_file = root_path + "/database/" + gn + "/.chrnames." + db\n+\n+ try:\n+ chrf = open(chrnames_file)\n+ except IOError:\n+ print "Open file: {0} error.".format(chrnames_file),\n+ print "Your database may be corrupted or you have an older version."\n+ sys.exit()\n+\n+ chr_list = chrf.read(1000000) # set 1MB limit to avoid nasty input.\n+ chrf.close()\n+\n+ chr_list = chr_list.strip()\n+ print chr_list\n+\n+\n+\n+###############################################\n+######## Main procedure. ######################\n+if __name__ == "__main__":\n+\n+ import argparse\n+ import os\n+ import sys\n+\n+ try:\n+ root_path = os.environ["NGSPLOT"]\n+ except KeyError:\n+ print "Environmental variable NGSPLOT is not set! Exit.\\n"\n+ sys.exit()\n+\n+ # Main parser.\n+ parser = argparse.ArgumentParser(description="Manipulate ngs.plot\'s \\\n+ annotation database",\n+ prog="ngsplotdb.py")\n+ parser.add_argument("-y", "--yes", help="Say yes to all prompted questions",\n+ action="store_true")\n+\n+ subparsers = parser.add_subparsers(title="Subcommands",\n+ help="additional help")\n+\n+ # ngsplotdb.py list parser.\n+ parser_list = subparsers.add_parser("list", help="List genomes installed \\\n+ in database")\n+ parser_list.set_defaults(func=listgn)\n+\n+ # ngsplotdb.py install parser.\n+ parser_install = subparsers.add_parser("install", \n+ help="Install genome from package \\\n+ file")\n+ parser_install.add_argument("pkg", help="Package file(.tar.gz) to install",\n+ type=str)\n+ parser_install.set_defaults(func=install)\n+\n+ # ngsplotdb.py remove parser.\n+ parser_remove = subparsers.add_parser("remove", \n+ help="Remove genome from database")\n+ parser_remove.add_argument("gn", help="Name of genome to be \\\n+ removed(e.g. hg18)", type=str)\n+ parser_remove.add_argument("--ftr", help="Remove cellline specific features \\\n+ (enhancer, dhs, etc)", type=str,\\\n+ default=None)\n+ parser_remove.set_defaults(func=remove)\n+\n+ # ngsplotdb.py chromosome names parser.\n+ parser_chrnames = subparsers.add_parser("chrnames", \n+ help="List chromosome names")\n+ parser_chrnames.add_argument("gn", help="Genome name(e.g. mm9)", type=str)\n+ parser_chrnames.add_argument("db", help="Database(ensembl or refseq)",\n+ type=str)\n+ parser_chrnames.set_defaults(func=chrnames)\n+\n+\n+ args = parser.parse_args()\n+ args.func(root_path, args)\n+\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c bin/remove.db.tables.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bin/remove.db.tables.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| @@ -0,0 +1,46 @@ +#!/usr/bin/env Rscript +# Remove the genome entries from two ngs.plot system files: +# default.tbl and dbfile.tbl. + +args <- commandArgs(T) +if(length(args)>1){ + gname.to.rm <- args[1] + feature.to.rm <- args[2] +}else{ + gname.to.rm <- args[1] +} + +# Save to text output. +progpath <- Sys.getenv('NGSPLOT') +if(progpath == "") { + stop("Set environment variable NGSPLOT before run the program. See README for details.\n") +} +default.file <- file.path(progpath, 'database', 'default.tbl') +dbfile.file <- file.path(progpath, 'database', 'dbfile.tbl') +default.tbl <- read.delim(default.file) +dbfile.tbl <- read.delim(dbfile.file) + +if(exists("feature.to.rm")){ + default.tbl <- subset(default.tbl, !(Genome==gname.to.rm & Region==feature.to.rm)) + dbfile.tbl <- subset(dbfile.tbl, !(Genome==gname.to.rm & Region==feature.to.rm)) +}else{ + default.tbl <- default.tbl[default.tbl$Genome != gname.to.rm, ] + dbfile.tbl <- dbfile.tbl[dbfile.tbl$Genome != gname.to.rm, ] +} + + +write.table(default.tbl, file=default.file, row.names=F, sep="\t", quote=F) +write.table(dbfile.tbl, file=dbfile.file, row.names=F, sep="\t", quote=F) + + +# getScore <- function(i, db.fi, default.fi){ +# score <- 3 - length(intersect(as.matrix(default.fi[db.fi[i, "tag"], ]), as.matrix(db.fi[i, ]))) +# return(score) +# } + +# for (i in (1:dim(db.fi)[1])){ +# anno.db.tbl[i, "dbScore"] <- getScore(i, db.fi, default.fi) +# } + +# anno.version="0.01" +# save(anno.tbl, anno.db.tbl, anno.version, file="database.RData") |
| b |
| diff -r 000000000000 -r 3ca58369469c bin/setTableDefaults.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bin/setTableDefaults.py Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| @@ -0,0 +1,91 @@ +#!/usr/bin/env python + +import sys +import string + +flankDict = { + "tss":2000, + "tes":2000, + "genebody":2000, + "exon":500, + "cgi":500, + "enhancer":1500, + "dhs":1000, +} + +pointLabDict = { + "tss":"TSS", + "tes":"TES", + "genebody":"TSS,TES", + "exon":"Acceptor,Donor", + "cgi":"Left,Right", + "enhancer":"Enhancer", + "dhs":"Left,Right", +} + +# FItype = { +# "tss":"endsite", +# "tes":"endsite", +# "genebody":"datatype", +# "exon":"subregion", +# "cgi":"subregion", +# "enhancer":"subregion", +# "dhs":"subregion" +# } + +in_f = file(sys.argv[1]) +out_f = file(sys.argv[2], "w") +out_f.write("\t".join(["Genome", "DefaultDB", "Region", "DefaultFI1", \ + "DefaultFI2", "DefaultFI3", "PointLab", "Flank"]) + \ + "\n") + +# while(True): +for line in in_f: + # line = in_f.readline() + # if len(line) == 0: + # break; + # line = string.strip(line) + + lineL = line.rstrip().split("\t") + genome = lineL[0] + defaultDB = lineL[1] + region = lineL[2] + + if region == "cgi": + fi_1 = "NA" + fi_2 = "ProximalPromoter" + fi_3 = "protein_coding" + elif region == "dhs": + fi_1 = "H1hesc" + fi_2 = "ProximalPromoter" + fi_3 = "protein_coding" + elif region == "exon": + fi_1 = "chipseq" + fi_2 = "canonical" + fi_3 = "protein_coding" + elif region == "genebody": + fi_1 = "chipseq" + fi_2 = "NA" + fi_3 = "protein_coding" + elif (region == "enhancer") and (genome == "hg19"): + fi_1 = "H1hesc" + fi_2 = "genebody" + fi_3 = "protein_coding" + elif (region == "enhancer") and (genome == "mm9"): + fi_1 = "mESC" + fi_2 = "genebody" + fi_3 = "protein_coding" + out_f.write("\t".join([genome, defaultDB, region, fi_1, fi_2, fi_3, \ + pointLabDict[region], str(flankDict[region])]) + "\n") + + # Extra: TSS and TES if region = genebody. + if region == "genebody": + region = "tss" + out_f.write("\t".join([genome, defaultDB, region, fi_1, fi_2, fi_3, \ + pointLabDict[region], str(flankDict[region])]) + "\n") + region = "tes" + out_f.write("\t".join([genome, defaultDB, region, fi_1, fi_2, fi_3, \ + pointLabDict[region], str(flankDict[region])]) + "\n") + +in_f.close() +out_f.close() |
| b |
| diff -r 000000000000 -r 3ca58369469c database/bootstrap/dbfile.tbl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/bootstrap/dbfile.tbl Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +db.file Genome DB Region FI.1 FI.2 FI.3 URL dbScore |
| b |
| diff -r 000000000000 -r 3ca58369469c database/bootstrap/default.tbl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/bootstrap/default.tbl Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +Genome DefaultDB Region DefaultFI1 DefaultFI2 DefaultFI3 PointLab Flank |
| b |
| diff -r 000000000000 -r 3ca58369469c database/bootstrap/gn_list.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/bootstrap/gn_list.txt Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +ID Assembly Species EnsVer NPVer |
| b |
| diff -r 000000000000 -r 3ca58369469c database/dbfile.tbl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/dbfile.tbl Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +db.file Genome DB Region FI.1 FI.2 FI.3 URL dbScore |
| b |
| diff -r 000000000000 -r 3ca58369469c database/default.tbl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/default.tbl Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +Genome DefaultDB Region DefaultFI1 DefaultFI2 DefaultFI3 PointLab Flank |
| b |
| diff -r 000000000000 -r 3ca58369469c database/gn_list.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/database/gn_list.txt Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,1 @@ +ID Assembly Species EnsVer NPVer |
| b |
| diff -r 000000000000 -r 3ca58369469c lib/coverage.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib/coverage.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,701 @@\n+# Function to check if the range exceeds coverage vector boundaries.\n+checkBound <- function(start, end, range, chrlen){\n+ if(end + range > chrlen || start - range < 1)\n+ return(FALSE) # out of boundary.\n+ else\n+ return(TRUE)\n+}\n+\n+Bin <- function(v, n) {\n+# Function to bin a coverage vector.\n+# Args:\n+# v: coverage vector.\n+# n: number of bins.\n+# Return: vector of binned values.\n+\n+ bin.bound <- seq(1, length(v), length.out=n + 1)\n+ sapply(1:n, function(i) {\n+ mean(v[bin.bound[i]:bin.bound[i + 1]])\n+ })\n+}\n+\n+SplineRev3Sec <- function(cov.list, v.fls, pts.list, v.strand, algo="spline") {\n+# For 3 sections of continuous coverage, spline each according to specified \n+# data points and return concatenated, interpolated curves.\n+# Args:\n+# cov.list: a list of coverage vectors. Each vector represents a gene.\n+# v.fls: vector of flanking region size.\n+# pts.list: a list of three integers for data points.\n+# v.strand: factor vector of strands.\n+# algo: algorithm used to normalize coverage vectors (spline, bin).\n+# Return: matrix of interpolated coverage: each row represents a gene; each \n+# column represents a data point. Coverage from exons are concatenated.\n+# Notes: the names of cov.list and pts.list must be the same for index purpose.\n+# Suggested: use \'l\', \'m\', and \'r\' for the names.\n+\n+ # Create an empty coveage matrix first.\n+ tot.pts <- sum(unlist(pts.list)) - 2\n+ cov.mat <- matrix(nrow=length(cov.list), ncol=tot.pts)\n+\n+ for(i in 1:length(cov.list)) {\n+ left.cov <- head(cov.list[[i]], n=v.fls[i])\n+ right.cov <- tail(cov.list[[i]], n=v.fls[i])\n+ mid.cov <- window(cov.list[[i]], start=v.fls[i] + 1, \n+ width=length(cov.list[[i]]) - 2*v.fls[i])\n+ if(algo == "spline") {\n+ left.cov <- spline(1:length(left.cov), left.cov, n=pts.list$l)$y\n+ right.cov <- spline(1:length(right.cov), right.cov, n=pts.list$r)$y\n+ mid.cov <- spline(1:length(mid.cov), mid.cov, n=pts.list$m)$y\n+ } else if(algo == "bin") {\n+ left.cov <- Bin(left.cov, n=pts.list$l)\n+ right.cov <- Bin(right.cov, n=pts.list$r)\n+ mid.cov <- Bin(mid.cov, n=pts.list$m)\n+ } else {\n+ # pass.\n+ }\n+ # Fuse the two points at the boundary and concatenate.\n+ f1 <- (tail(left.cov, n=1) + head(mid.cov, n=1)) / 2\n+ f2 <- (tail(mid.cov, n=1) + head(right.cov, n=1)) / 2\n+ con.cov <- c(left.cov[-length(left.cov)], f1, \n+ mid.cov[-c(1, length(mid.cov))], \n+ f2, right.cov[-1])\n+ # browser()\n+ if(v.strand[i] == \'+\') {\n+ cov.mat[i, ] <- con.cov\n+ } else {\n+ cov.mat[i, ] <- rev(con.cov)\n+ }\n+ }\n+\n+ cov.mat\n+}\n+\n+extrCov3Sec <- function(v.chrom, v.start, v.end, v.fls, v.strand, m.pts, f.pts, \n+ bufsize, algo, ...) {\n+# Extract and interpolate coverage vectors from genomic regions with 3 sections.\n+# Args:\n+# v.chrom: factor vector of chromosome names.\n+# v.start: integer vector of region start.\n+# v.end: integer vector of region end.\n+# v.fls: integer vector of flanking region size in bps.\n+# v.strand: factor vector of gene strands.\n+# m.pts: data points for middle interval.\n+# f.pts: data points for flanking region.\n+# bufsize: integer; buffers are added to both ends of each region.\n+# algo: algorithm used to normalize coverage vectors.\n+# Return: matrix of interpolated coverage: each row represents a gene; each \n+# column represents a data point.\n+\n+ interflank.gr <- GRanges(seqnames=v.chrom, ranges=IRanges(\n+ start=v.start - v.fls - bufsize, \n+ end=v.end + v.fls + bufsize))\n+ interflank.cov <- covBamExons(interflank.gr, v.strand, ...)\n+\n+ # Trim buffers from coverage vectors.\n+ interflank.cov <- lapply(interfla'..b' <- chk[1]:chk[2] # chunk: start -> end\n+ # If RNA-seq, retrieve exon ranges.\n+ if(rnaseq.gb) {\n+ exonranges.list <- unlist(exonmodel[coord[i, ]$tid])\n+ } else {\n+ exonranges.list <- NULL\n+ }\n+ doCov(coord[i, ], exonranges.list, ...)\n+ }\n+\n+ # Floor negative values which are caused by spline.\n+ result.matrix[result.matrix < 0] <- 0\n+ result.matrix / libsize * 1e6 # normalize to RPM.\n+\n+ } else {\n+ for(c in 1:length(chkidx.list)) {\n+ chk <- chkidx.list[[c]]\n+ i <- chk[1]:chk[2] # chunk: start -> end\n+ cat(".")\n+ # If RNA-seq, retrieve exon ranges.\n+ if(rnaseq.gb) {\n+ exonranges.list <- unlist(exonmodel[coord[i, ]$tid])\n+ } else {\n+ exonranges.list <- NULL\n+ }\n+ cov <- doCov(coord[i, ], exonranges.list, ...)\n+ if(c == 1) {\n+ result.matrix <- matrix(0, nrow=nrow(coord), ncol=ncol(cov))\n+ }\n+ result.matrix[i, ] <- cov\n+ }\n+ # Floor negative values which are caused by spline.\n+ # browser()\n+ result.matrix[result.matrix < 0] <- 0\n+ result.matrix / libsize * 1e6 # normalize to RPM.\n+\n+ }\n+}\n+\n+\n+doCov <- function(coord.mat, exonranges.list, chr.tag, pint, reg2plot, \n+ flanksize, flankfactor, m.pts, f.pts, ...) {\n+# Extract coverage from bam file into a matrix. According to the parameter\n+# values, call corresponding functions.\n+# Args:\n+# coord.mat: matrix of genomic coordinates to extract coverage.\n+# exonranges.list: list of IRanges objects, each represents a group of exons.\n+# pint: boolean of point interval.\n+# reg2plot: string of region to plot.\n+# flanksize: flanking region size.\n+# flankfactor: flanking region factor.\n+# m.pts: data points for middle interval.\n+# f.pts: data points for flanking region.\n+# Return: matrix of interpolated coverage: each row represents a gene; each \n+# column represents a data point. Coverage from exons are concatenated.\n+\n+ v.chrom <- coord.mat$chrom\n+ # if(!chr.tag) {\n+ # v.chrom <- sub(\'chr\', \'\', v.chrom)\n+ # }\n+ v.chrom <- as.factor(v.chrom)\n+ v.strand <- as.factor(coord.mat$strand)\n+\n+ # Figure out interval region sizes and calculate flanking region sizes.\n+ if(!pint) { # interval regions.\n+ if(!is.null(exonranges.list)) { # RNA-seq\n+ if(flankfactor > 0) {\n+ v.fls <- sapply(exonranges.list, function(t) {\n+ sum(end(t) - start(t) + 1) * flankfactor\n+ })\n+ } else {\n+ v.fls <- rep(flanksize, length=nrow(coord.mat))\n+ }\n+ extrCovExons(v.chrom, exonranges.list, v.fls, v.strand, \n+ m.pts, f.pts, ...)\n+ } else { # ChIP-seq with intervals.\n+ v.start <- coord.mat$start\n+ v.end <- coord.mat$end\n+ if(flankfactor > 0) {\n+ v.fls <- round((v.end - v.start + 1) * flankfactor)\n+ } else {\n+ v.fls <- rep(flanksize, length=nrow(coord.mat))\n+ }\n+ extrCov3Sec(v.chrom, v.start, v.end, v.fls, v.strand, \n+ m.pts, f.pts, ...)\n+ }\n+ } else { # point center with flanking regions.\n+ v.midp <- vector(\'integer\', length=nrow(coord.mat))\n+ for(r in 1:nrow(coord.mat)) {\n+ if(reg2plot == \'tss\' && v.strand[r] == \'+\' || \n+ reg2plot == \'tes\' && v.strand[r] == \'-\') {\n+ # the left site is center.\n+ v.midp[r] <- coord.mat$start[r]\n+ } else { # this also includes BED supplied point center.\n+ v.midp[r] <- coord.mat$end[r]\n+ }\n+ }\n+ # browser()\n+ extrCovMidp(v.chrom, v.midp, flanksize, v.strand, m.pts + f.pts*2, ...)\n+ }\n+}\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c lib/genedb.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib/genedb.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,258 @@\n+chromFormat <- function(crn, ...){\n+# Format chromosome names to our standard.\n+# Args:\n+# crn: character vector of chromosome names.\n+# Returns: character vector of formatted chromosome names.\n+\n+ crn <- sub(\'.fa$\', \'\', crn)\n+ nochr.i <- grep(\'^chr\', crn, invert=T)\n+ crn[nochr.i] <- paste(\'chr\', crn[nochr.i], sep=\'\')\n+ crn\n+}\n+\n+FIScoreIntersect <- function(db.info, v.finfo){\n+# Further info score based on intersection with database FI columns.\n+# Used to select database files.\n+# Args: \n+# db.info: one record of annotation database table.\n+# v.finfo: further info that is split into vector.\n+# Returns: a score for the intersection between database and further info.\n+\n+ length(intersect(as.vector(db.info[c("FI.1", "FI.2", "FI.3")]), v.finfo))\n+}\n+\n+SetupPlotCoord <- function(args.tbl, ctg.tbl, default.tbl, dbfile.tbl, \n+ progpath, genome, reg2plot, lgint, flanksize, \n+ samprate, galaxy) {\n+# Load genomic coordinates for plot based on the input arguments.\n+# Args:\n+# args.tbl: input argument table\n+# ctg.tbl: coverage-genelist-title table\n+# default.tbl: default settings of databases and plot setting\n+# dbfile.tbl: details of database files(.RData).\n+# progpath: program path derived from NGSPLOT\n+# genome: genome name, such as mm9, hg19.\n+# reg2plot: tss, tes, genebody, etc.\n+# lgint: boolean tag for large interval.\n+# flanksize: flanking region size.\n+# samprate: sampling rate\n+# galaxy: boolean tag for galaxy installation.\n+\n+ if(reg2plot!=\'bed\'){\n+ # Subset using genome-region combination.\n+ key <- default.tbl$Genome == genome & default.tbl$Region == reg2plot\n+ if(sum(key) == 0) {\n+ stop("The combination of genome and region does not exist. You \n+ may need to install the genome or the region does not exist yet.\\n")\n+ }\n+ anno.parameters <- default.tbl[key, ]\n+ db.match.mask <- dbfile.tbl$Genome == genome & \n+ dbfile.tbl$Region == reg2plot\n+ anno.db.candidates <- dbfile.tbl[db.match.mask, ]\n+ \n+ # Database flavor.\n+ if(\'-D\' %in% names(args.tbl)){ \n+ database <- as.character(args.tbl[\'-D\'])\n+ database.allowed <- unique(anno.db.candidates$DB)\n+ stopifnot(database %in% database.allowed)\n+ }else{\n+ database <- as.character(anno.parameters$DefaultDB)\n+ }\n+\n+ anno.db.candidates <- anno.db.candidates[anno.db.candidates$DB == database, ]\n+\n+ if (file.exists(file.path(progpath, \'database\', genome, reg2plot))){\n+ prefix <- file.path(progpath, \'database\', genome, reg2plot)\n+ }else{\n+ prefix <- file.path(progpath, \'database\', genome)\n+ }\n+\n+ # Further info to subset genomic regions.\n+ if(\'-F\' %in% names(args.tbl)){\n+ finfo <- as.character(args.tbl[\'-F\'])\n+ } else {\n+ finfo <- NULL\n+ }\n+\n+ Labs <- unlist(strsplit(as.character(anno.parameters$PointLab[1]), ","))\n+ if(length(Labs)==1){\n+ pint <- TRUE\n+ }else{\n+ pint <- FALSE\n+ }\n+\n+ if(!is.null(finfo)){ # use finfo to subset DB tables.\n+ v.finfo <- unlist(strsplit(finfo, \',\'))\n+ # Get the intersect of the v.finfo and FI columns of databases, \n+ # and choose the best fit.\n+ fi.score <- apply(anno.db.candidates, 1, FIScoreIntersect, \n+ v.finfo=v.finfo)\n+ anno.db.candidates <- anno.db.candidates[fi.score == max(fi.score), ]\n+ # RNA-seq tag.\n+ if(reg2plot == \'genebody\' && \'rnaseq\' %in% v.finfo) {\n+ rnaseq.gb <- TRUE\n+ }else{\n+ rnaseq.gb <- FALSE\n+ }\n+ } else {\n+ rnaseq.gb <- FALSE\n+ }\n+ anno.db.candidates <- anno.db.candidates[order(anno.db.candidates$dbScore), ]\n+ '..b' gid=NA, gname=\'N\', tid=\'N\', strand=\'+\', \n+ byname.uniq=T, bygid.uniq=NA)\n+ # Perform sanity check for bed file.\n+ if(!all(genome.coord$start <= genome.coord$end)) {\n+ stop(sprintf("Sanity check for bed file: %s failed. Bed files are \n+ 0-based and right-side open.\\n", bed.file))\n+ }\n+\n+ # Deal with columns 4-6 (Optional).\n+ if(ncol(bed.coord) >=4){\n+ genome.coord$gname <- bed.coord[, 4]\n+ }\n+ if(ncol(bed.coord) >=5){\n+ genome.coord$tid <- bed.coord[, 5]\n+ }\n+ if(ncol(bed.coord) >=6){\n+ genome.coord$strand <- bed.coord[, 6]\n+ }\n+\n+ genome.coord\n+ }\n+\n+ # Sample a vector but preserve its original sequential order.\n+ sampleInSequence <- function(x, samprate) {\n+ x[ifelse(runif(length(x)) <= samprate, T, F)]\n+ }\n+\n+ ValidateGeneSubset <- function(gid.match, tid.match, gname.match) {\n+ # To validate that there is no ambiguous hits in gene match.\n+ # Args:\n+ # XXX.match: positional hits with respect to the gene database. 0 means \n+ # no hit.\n+\n+ subset.matrix <- rbind(gid.match, tid.match, gname.match)\n+ g.valid <- sapply(subset.matrix, function(col.hits) {\n+ length(unique(col.hits[col.hits > 0])) <= 1\n+ })\n+\n+ all(g.valid)\n+ }\n+\n+ for(i in 1:length(uniq.reg)) {\n+ ur <- uniq.reg[i]\n+ \n+ if(reg2plot == "bed") {\n+ coord.list[[i]] <- ReadBedCoord(ur)\n+ } else if(ur == \'-1\') { # use whole genome.\n+ coord.list[[i]] <- genome.coord[genome.coord$byname.uniq, ]\n+ } else { # read gene list from text file.\n+ gene.list <- read.table(ur, as.is=T, comment.char=\'#\')$V1\n+ gid.match <- match(gene.list, genome.coord$gid, nomatch=0)\n+ gname.match <- match(gene.list, genome.coord$gname, nomatch=0)\n+ tid.match <- match(gene.list, genome.coord$tid, nomatch=0)\n+ if(!ValidateGeneSubset(gid.match, tid.match, gname.match)) {\n+ stop("\\nAmbiguous hits in gene database. This shall never \n+ happen! Contact ngs.plot maintainers.\\n")\n+ }\n+ subset.idx <- pmax(gid.match, tid.match, gname.match)\n+ # subset.idx <- gid.match + tid.match + gname.match\n+ subset.idx <- subset.idx[subset.idx != 0]\n+ if(length(subset.idx) == 0) {\n+ stop("\\nGene subset size becomes zero. Are you using the \n+ correct database?\\n")\n+ }\n+ coord.list[[i]] <- genome.coord[subset.idx, ]\n+ }\n+\n+ # Do sampling.\n+ if(samprate < 1) {\n+ samp.idx <- sampleInSequence(1:nrow(coord.list[[i]]), samprate)\n+ coord.list[[i]] <- coord.list[[i]][samp.idx, ]\n+ }\n+ }\n+\n+ # If bed file, set boolean tag for point interval, i.e. interval=1bp.\n+ if(reg2plot == "bed") {\n+ pint <- all(sapply(coord.list, function(cd) all(cd$start == cd$end)))\n+ if(pint) {\n+ Labs <- c("Center")\n+ } else {\n+ Labs <- c("5\'End", "3\'End")\n+ }\n+ }\n+\n+ # Set boolean tag for large interval display.\n+ if(!pint && is.na(lgint)) {\n+ if(median(sapply(coord.list, function(cd) {\n+ median(cd$end - cd$start + 1)\n+ })) > flanksize) {\n+ lgint <- 1\n+ } else {\n+ lgint <- 0\n+ }\n+ }\n+\n+ res <- list(coord.list=coord.list, rnaseq.gb=rnaseq.gb, lgint=lgint,\n+ reg.list=reg.list, uniq.reg=uniq.reg, pint=pint, Labs=Labs)\n+\n+ # Add exon models to the result if RNA-seq.\n+ if(rnaseq.gb) {\n+ em.load <- file.path(prefix, paste(genome, database, \'exonmodel.RData\', \n+ sep=\'.\'))\n+ load(em.load) # load exon models into variable: exonmodel\n+ res$exonmodel <- exonmodel\n+ }\n+\n+ res\n+}\n+\n+\n+\n+\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c lib/parse.args.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib/parse.args.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,521 @@\n+# Parse command line arguments and extract them into an associate array.\n+# Check if the required arguments are all satisfied.\n+parseArgs <- function(args, manditories) {\n+ if(length(args) %% 2 == 1 || length(args) == 0) {\n+ cat(\'Unpaired argument and value.\\n\')\n+ return(NULL)\n+ }\n+ n.i <- seq(1, length(args), by=2)\n+ v.i <- seq(2, length(args), by=2)\n+ args.name <- args[n.i]\n+ args.value <- args[v.i]\n+\n+ # Check if required argument values are supplied.\n+ miss_tag <- F\n+ man.bool <- manditories %in% args.name\n+ if(!all(man.bool)){\n+ cat(paste(\'Missing argument: \', paste(manditories[!man.bool], \n+ collapse=\',\'), \'.\', sep=\'\')\n+ )\n+ miss_tag <- T\n+ }\n+ if(miss_tag){\n+ res <- NULL\n+ }else{\n+ res <- args.value\n+ names(res) <- args.name\n+ }\n+ res\n+}\n+\n+ConfigTbl <- function(args.tbl, fraglen) {\n+# Create configuration table from "-C" argument.\n+# Args:\n+# args.tbl: named vector of program arguments.\n+# fraglen: fragment length.\n+# Returns: dataframe of configuration.\n+\n+ covfile <- args.tbl[\'-C\']\n+\n+ suppressWarnings(\n+ ctg.tbl <- tryCatch(\n+ read.table(covfile, colClasses=\'character\', comment.char=\'#\'), \n+ error=function(e) {\n+ if(\'-E\' %in% names(args.tbl)) {\n+ glist <- args.tbl[\'-E\']\n+ } else {\n+ glist <- \'-1\'\n+ }\n+ if(\'-T\' %in% names(args.tbl)) {\n+ title <- args.tbl[\'-T\']\n+ } else {\n+ title <- \'Noname\'\n+ }\n+ data.frame(cov=covfile, glist=glist, title=title, \n+ fraglen=as.character(fraglen), \n+ color=NA, stringsAsFactors=F)\n+ }\n+ )\n+ )\n+\n+ # Read a config file.\n+ if(ncol(ctg.tbl) < 3) {\n+ stop("Configuration file must contain at least 3 columns! \n+Insufficient information provided.\\n")\n+ }\n+ colnames(ctg.tbl)[1:3] <- c(\'cov\', \'glist\', \'title\')\n+ if(ncol(ctg.tbl) >= 4) {\n+ colnames(ctg.tbl)[4] <- \'fraglen\'\n+ fraglen.sp <- strsplit(ctg.tbl$fraglen, ":")\n+ if(!all(sapply(fraglen.sp, function(x) {\n+ length(x) == 1 || length(x) == 2}))) {\n+ stop("Fragment length format must be X or X:Y; X and Y are \n+integers.\\n")\n+ }\n+ if(!all(as.integer(unlist(fraglen.sp)) > 0)) {\n+ stop("Fragment length must be positive integers! Check your \n+configuration file.\\n")\n+ }\n+ } else {\n+ ctg.tbl <- data.frame(ctg.tbl, fraglen=as.character(fraglen),\n+ stringsAsFactors=F)\n+ }\n+ if(ncol(ctg.tbl) >= 5) {\n+ colnames(ctg.tbl)[5] <- \'color\'\n+ # Validate color specifications.\n+ col.validated <- col2rgb(ctg.tbl$color)\n+ } else {\n+ ctg.tbl <- data.frame(ctg.tbl, color=NA)\n+ }\n+ ctg.tbl\n+}\n+\n+CheckRegionAllowed <- function(reg2plot, anno.tbl) {\n+# Check if region to plot is an allowed value.\n+\n+ region.allowed <- c(as.vector(unique(anno.tbl$Region)), "bed")\n+ if(!reg2plot %in% region.allowed) {\n+ stop(paste(c("Unknown region specified. Must be one of:", \n+ region.allowed, "\\n"), collapse=" "))\n+ }\n+}\n+\n+CoverageVars <- function(args.tbl, reg2plot) {\n+# Setup variables from program arguments.\n+# Args:\n+# args.tbl: named vector of program arguments.\n+# reg2plot: string describing region to plot.\n+# Returns: list of variables.\n+\n+ vl <- list() # vl: configured variable list\n+\n+ #### Switch for debug ####\n+ if(\'-Debug\' %in% names(args.tbl)) {\n+ stopifnot(as.integer(args.tbl[\'-Debug\']) >= 0)\n+ vl$debug <- as.integer(args.tbl[\'-Debug\'])\n+ } else {\n+ vl$debug <- as.integer(0)\n+ }\n+\n+ #### Switch for Galaxy usage ####\n+ if(\'-Galaxy\' %in% names(args.tbl)) {\n'..b'window width to smooth avg. profiles, must be integer\\n")\n+ cat(" 1=no(default); 3=slightly; 5=somewhat; 9=quite; 13=super.\\n")\n+ cat(" -H Opacity of shaded area, suggested value:[0, 0.5]\\n")\n+ cat(" default=0, i.e. no shading, just lines\\n")\n+ cat(" -YAS Y-axis scale: auto(default) or min_val,max_val(custom scale)\\n")\n+ cat(" -LEG Draw legend? 1(default) or 0\\n")\n+ cat(" -BOX Draw box around plot? 1(default) or 0\\n")\n+ cat(" -VLN Draw vertical lines? 1(default) or 0\\n")\n+ cat(" -XYL Draw X- and Y-axis labels? 1(default) or 0\\n")\n+ cat(" -LWD Line width(default=3)\\n")\n+ cat("### Heatmap parameters:\\n")\n+ cat(" -GO Gene order algorithm used in heatmaps: total(default), hc, max,\\n")\n+ cat(" prod, diff, km and none(according to gene list supplied)\\n")\n+ cat(" -LOW Low count cutoff(default=10) in rank-based normalization\\n")\n+ cat(" -KNC K-means or HC number of clusters(default=5)\\n")\n+ cat(" -MIT Maximum number of iterations(default=20) for K-means\\n")\n+ cat(" -NRS Number of random starts(default=30) in K-means\\n")\n+ cat(" -RR Reduce ratio(default=30). The parameter controls the heatmap height\\n")\n+ cat(" The smaller the value, the taller the heatmap\\n")\n+ cat(" -SC Color scale used to map values to colors in a heatmap\\n")\n+ cat(" local(default): base on each individual heatmap\\n")\n+ cat(" region: base on all heatmaps belong to the same region\\n")\n+ cat(" global: base on all heatmaps together\\n")\n+ cat(" min_val,max_val: custom scale using a pair of numerics\\n")\n+ cat(" -FC Flooding fraction:[0, 1), default=0.02\\n")\n+ cat(" -CO Color for heatmap. For bam-pair, use color-tri(neg_color:[neu_color]:pos_color)\\n")\n+ cat(" Hint: must use R colors, such as darkgreen, yellow and blue2\\n")\n+ cat(" The neutral color is optional(default=black)\\n")\n+ cat(" -CD Color distribution for heatmap(default=0.6). Must be a positive number\\n")\n+ cat(" Hint: lower values give more widely spaced colors at the negative end\\n")\n+ cat(" In other words, they shift the neutral color to positive values\\n")\n+ cat(" If set to 1, the neutral color represents 0(i.e. no bias)\\n")\n+\n+}\n+\n+EchoCoverageArgs <- function() {\n+ cat("## Coverage-generation parameters:\\n")\n+ cat(" -F Further information provided to select database table or plottype:\\n")\n+ cat(" This is a string of description separated by comma.\\n")\n+ cat(" E.g. protein_coding,K562,rnaseq(order of descriptors does not matter)\\n")\n+ cat(" means coding genes in K562 cell line drawn in rnaseq mode.\\n")\n+ cat(" -D Gene database: ensembl(default), refseq\\n")\n+ cat(" -L Flanking region size(will override flanking factor)\\n")\n+ cat(" -N Flanking region factor\\n")\n+ cat(" -RB The fraction of extreme values to be trimmed on both ends\\n")\n+ cat(" default=0, 0.05 means 5% of extreme values will be trimmed\\n")\n+ cat(" -S Randomly sample the regions for plot, must be:(0, 1]\\n")\n+ cat(" -P #CPUs to use. Set 0(default) for auto detection\\n")\n+ cat(" -AL Algorithm used to normalize coverage vectors: spline(default), bin\\n")\n+ cat(" -CS Chunk size for loading genes in batch(default=100)\\n")\n+ cat(" -MQ Mapping quality cutoff to filter reads(default=20)\\n")\n+ cat(" -FL Fragment length used to calculate physical coverage(default=150)\\n")\n+ cat(" -SS Strand-specific coverage calculation: both(default), same, opposite\\n")\n+ cat(" -IN Shall interval be larger than flanking in plot?(0 or 1, default=automatic)\\n")\n+ cat(" -FI Forbid image output if set to 1(default=0)\\n")\n+}\n+############### End arguments configuration #####################\n+#################################################################\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c lib/plotlib.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib/plotlib.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b"@@ -0,0 +1,649 @@\n+#### plotlib.r ####\n+# This contains the library for plotting related functions.\n+#\n+# Authors: Li Shen, Ningyi Shao\n+# \n+# Created: Feb 19, 2013\n+# Last updated: May 21, 2013\n+#\n+\n+SetupHeatmapDevice <- function(reg.list, uniq.reg, ng.list, pts, font.size=12,\n+ unit.width=4, reduce.ratio=30) {\n+# Configure parameters for heatmap output device. The output is used by \n+# external procedures to setup pdf device ready for heatmap plotting.\n+# Args:\n+# reg.list: region list as in config file.\n+# uniq.reg: unique region list.\n+# ng.list: number of genes per heatmap in the order as config file.\n+# pts: data points (number of columns of heatmaps).\n+# font.size: font size.\n+# unit.width: image width per heatmap.\n+# reduce.ratio: how compressed are genes in comparison to data points? This \n+# controls image height.\n+\n+ # Number of plots per region.\n+ reg.np <- sapply(uniq.reg, function(r) sum(reg.list==r))\n+\n+ # Number of genes per region.\n+ reg.ng <- sapply(uniq.reg, function(r) {\n+ ri <- which(reg.list==r)[1]\n+ ng.list[ri]\n+ })\n+\n+ # Adjustment ratio.\n+ origin.fs <- 12 # default font size.\n+ fs.adj.ratio <- font.size / origin.fs\n+ # Margin size (in lines) adjusted by ratio.\n+ m.bot <- fs.adj.ratio * 2\n+ m.lef <- fs.adj.ratio * 1.5\n+ m.top <- fs.adj.ratio * 2\n+ m.rig <- fs.adj.ratio * 1.5 \n+ key.in <- fs.adj.ratio * 1.0 # colorkey in inches.\n+ m.lef.diff <- (fs.adj.ratio - 1) * 1.5\n+ m.rig.diff <- (fs.adj.ratio - 1) * 1.5\n+ # Setup image size.\n+ hm.width <- (unit.width + m.lef.diff + m.rig.diff) * max(reg.np)\n+ ipl <- .2 # inches per line. Obtained from par->'mai', 'mar'.\n+ # Convert #gene to image height.\n+ reg.hei <- sapply(reg.ng, function(r) {\n+ c(key.in, # colorkey + margin.\n+ r * unit.width / pts / reduce.ratio + \n+ m.bot * ipl + m.top * ipl) # heatmap + margin.\n+ })\n+ reg.hei <- c(reg.hei)\n+ hm.height <- sum(reg.hei)\n+\n+ # Setup layout of the heatmaps.\n+ lay.mat <- matrix(0, ncol=max(reg.np), nrow=length(reg.np) * 2)\n+ fig.n <- 1 # figure plotting number.\n+ for(i in 1:length(reg.np)) {\n+ row.upper <- i * 2 - 1\n+ row.lower <- i * 2\n+ for(j in 1:reg.np[i]) {\n+ lay.mat[row.upper, j] <- fig.n;\n+ fig.n <- fig.n + 1\n+ lay.mat[row.lower, j] <- fig.n;\n+ fig.n <- fig.n + 1\n+ }\n+ }\n+\n+ list(reg.hei=reg.hei, hm.width=hm.width, hm.height=hm.height, \n+ lay.mat=lay.mat, heatmap.mar=c(m.bot, m.lef, m.top, m.rig) * ipl)\n+}\n+\n+SetPtsSpline <- function(pint, lgint) {\n+# Set data points for spline function.\n+# Args:\n+# pint: tag for point interval.\n+# Return: list of data points, middle data points, flanking data points.\n+\n+ pts <- 100 # data points to plot: 0...pts\n+ if(pint){ # point interval.\n+ m.pts <- 1 # middle part points.\n+ f.pts <- pts / 2 # flanking part points.\n+ } else {\n+ if(lgint) {\n+ m.pts <- pts / 5 * 3 + 1\n+ f.pts <- pts / 5 + 1\n+ } else {\n+ m.pts <- pts / 5 + 1\n+ f.pts <- pts / 5 * 2 + 1\n+ }\n+ }\n+ list(pts=pts, m.pts=m.pts, f.pts=f.pts)\n+}\n+\n+CreatePlotMat <- function(pts, ctg.tbl) {\n+# Create matrix for avg. profiles.\n+# Args:\n+# pts: data points.\n+# ctg.tbl: configuration table.\n+# Return: avg. profile matrix initialized to zero.\n+\n+ regcovMat <- matrix(0, nrow=pts + 1, ncol=nrow(ctg.tbl))\n+ colnames(regcovMat) <- ctg.tbl$title\n+ regcovMat\n+}\n+\n+CreateConfiMat <- function(se, pts, ctg.tbl){\n+# Create matrix for standard errors.\n+# Args:\n+# se: tag for standard error plotting.\n+# pts: data points.\n+# ctg.tbl: configuration table.\n+# Return: standard error matrix initialized to zero or null.\n+\n+ if(se){\n+ confiMat <- matrix(0, nrow"..b'richCombined <- do.call(\'cbind\', enrichList[plist])\n+ enrichSelected <- enrichList[plist]\n+\n+ # If color scale is region, calculate breaks and quantile here.\n+ if(color.scale == \'region\') {\n+ flood.pts <- quantile(c(enrichSelected, recursive=T), \n+ c(flood.q, 1-flood.q))\n+ brk.use <- ColorBreaks(flood.pts[2], flood.pts[1], bam.pair, ncolor)\n+ }\n+\n+ # Order genes.\n+ if(is.matrix(enrichSelected[[1]]) && nrow(enrichSelected[[1]]) > 1) {\n+ if(bam.pair) {\n+ lowCutoffs <- 0\n+ } else {\n+ lowCutoffs <- v.low.cutoff[plist]\n+ }\n+ g.order.list <- OrderGenesHeatmap(enrichSelected, lowCutoffs, \n+ go.algo, go.paras)\n+ g.order <- g.order.list[[1]]\n+ g.cluster <- g.order.list[[2]]\n+ if(is.null(g.cluster)) {\n+ go.cluster[[ur]] <- NA\n+ } else{\n+ go.cluster[[ur]] <- rev(g.cluster[g.order])\n+ }\n+ go.list[[ur]] <- rev(rownames(enrichSelected[[1]][g.order, ]))\n+ } else {\n+ go.cluster[[ur]] <- NULL\n+ go.list[[ur]] <- NULL\n+ }\n+\n+ if(!do.plot) {\n+ next\n+ }\n+ \n+ # Go through each sample and do plot.\n+ for(pj in plist) {\n+ if(!is.null(g.order)) {\n+ enrichList[[pj]] <- enrichList[[pj]][g.order, ]\n+ }\n+\n+ # If color scale is local, calculate breaks and quantiles here.\n+ if(color.scale == \'local\') {\n+ flood.pts <- quantile(c(enrichList[[pj]], recursive=T), \n+ c(flood.q, 1-flood.q))\n+ brk.use <- ColorBreaks(flood.pts[2], flood.pts[1], bam.pair, \n+ ncolor)\n+ }\n+\n+ # Flooding extreme values.\n+ enrichList[[pj]][ enrichList[[pj]] < flood.pts[1] ] <- flood.pts[1]\n+ enrichList[[pj]][ enrichList[[pj]] > flood.pts[2] ] <- flood.pts[2]\n+\n+ # Draw colorkey.\n+ image(z=matrix(brk.use, ncol=1), col=enrich.palette(ncolor), \n+ breaks=brk.use, axes=F, useRaster=T, main=\'Colorkey\')\n+ axis(1, at=seq(0, 1, length.out=5), \n+ labels=format(brk.use[seq(1, ncolor + 1, length.out=5)], \n+ digits=1), \n+ lwd=1, lwd.ticks=1)\n+\n+ # Draw heatmap.\n+ image(z=t(enrichList[[pj]]), col=enrich.palette(ncolor), \n+ breaks=brk.use, axes=F, useRaster=T, main=title2plot[pj])\n+\n+ axis(1, at=xticks$pos, labels=xticks$lab, lwd=1, lwd.ticks=1)\n+ }\n+ }\n+ list(go.list,go.cluster)\n+}\n+\n+trim <- function(x, p){\n+# Trim a numeric vector on both ends.\n+# Args:\n+# x: numeric vector.\n+# p: percentage of data to trim.\n+# Return: trimmed vector.\n+\n+ low <- quantile(x, p)\n+ hig <- quantile(x, 1 - p)\n+ x[x > low & x < hig]\n+}\n+\n+CalcSem <- function(x, rb=.05){ \n+# Calculate standard error of mean for a numeric vector.\n+# Args:\n+# x: numeric vector\n+# rb: fraction of data to trim before calculating sem.\n+# Return: a scalar of the standard error of mean\n+\n+ if(rb > 0){\n+ x <- trim(x, rb)\n+ }\n+ sem <- sd(x) / sqrt(length(x))\n+ ifelse(is.na(sem), 0, sem)\n+ # NOTE: this should be improved to handle exception that "sd" calculation \n+ # emits errors.\n+}\n+\n+\n+\n+\n+## Leave for future reference.\n+#\n+# Set the antialiasing.\n+# type <- NULL\n+# if (capabilities()["aqua"]) {\n+# type <- "quartz"\n+# } else if (capabilities()["cairo"]) {\n+# type <- "cairo"\n+# } else if (capabilities()["X11"]) {\n+# type <- "Xlib"\n+# }\n+# Set the output type based on capabilities.\n+# if (is.null(type)){\n+# png(plot.name, width, height, pointsize=pointsize)\n+\n+# } else {\n+# png(plot.name, width, height, pointsize=pointsize, type=type)\n+# }\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c ngs.plot.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ngs.plot.r Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,341 @@\n+#!/usr/bin/env Rscript\n+#\n+# Program: ngs.plot.r\n+# Purpose: Plot sequencing coverages at different genomic regions.\n+# Allow overlaying various coverages with gene lists.\n+#\n+# -- by Li Shen, MSSM\n+# Created: Nov 2011.\n+# -- by Christophe Antoniewski\n+# recoded for Galaxy:\tDec 2017\n+\n+ngsplot.version <- \'2.63_artbio\'\n+# Program environment variable.\n+progpath <- Sys.getenv(\'NGSPLOT\')\n+if(progpath == "") {\n+ stop("Set environment variable NGSPLOT before run the program. See README \n+ for details.\\n")\n+}\n+\n+source(file.path(progpath, \'lib\', \'parse.args.r\'))\n+source(file.path(progpath, \'lib\', \'genedb.r\'))\n+source(file.path(progpath, \'lib\', \'plotlib.r\'))\n+source(file.path(progpath, \'lib\', \'coverage.r\'))\n+\n+# Deal with command line arguments.\n+cmd.help <- function(){\n+ cat("\\nVisit https://github.com/shenlab-sinai/ngsplot/wiki/ProgramArguments101 for details\\n")\n+ cat(paste("Version:", ngsplot.version, sep=" "))\n+ cat("\\nUsage: ngs.plot.r -G genome -R region -C [cov|config]file\\n")\n+ cat(" -O name [Options]\\n")\n+ cat("\\n## Mandatory parameters:\\n")\n+ cat(" -G Genome name. Use ngsplotdb.py list to show available genomes.\\n")\n+ cat(" -R Genomic regions to plot: tss, tes, genebody, exon, cgi, enhancer, dhs or bed\\n")\n+ cat(" -C Indexed bam file or a configuration file for multiplot\\n")\n+ cat(" -O Name for output: multiple files will be generated\\n")\n+ cat("## Optional parameters related to configuration file:\\n")\n+ cat(" -E Gene list to subset regions OR bed file for custom region\\n")\n+ cat(" -T Image title\\n")\n+ EchoCoverageArgs()\n+ EchoPlotArgs()\n+ cat("\\n")\n+}\n+\n+\n+###########################################################################\n+#################### Deal with program input arguments ####################\n+args <- commandArgs(T)\n+# args <- unlist(strsplit(\'-G mm10 -R tss -C K27M_no_stim_3_GATCAG_L006_R1_001Aligned.out.srt.rmdup.bam -O test -Debug 1\', \' \'))\n+\n+# Input argument parser.\n+args.tbl <- parseArgs(args, c(\'-G\', \'-C\', \'-R\', \'-O\')) # see lib/parse.args.r\n+if(is.null(args.tbl)){\n+ cmd.help()\n+ stop(\'Error in parsing command line arguments. Stop.\\n\')\n+}\n+genome <- args.tbl[\'-G\']\n+reg2plot <- args.tbl[\'-R\']\n+oname <- args.tbl[\'-O\']\n+warning(sprintf("%s", args.tbl))\n+\n+cat("Configuring variables...")\n+# Load tables of database: default.tbl, dbfile.tbl\n+default.tbl <- read.delim(file.path(progpath, \'database\', \'default.tbl\'))\n+dbfile.tbl <- read.delim(file.path(progpath, \'database\', \'dbfile.tbl\'))\n+CheckRegionAllowed(reg2plot, default.tbl)\n+\n+# Setup variables from arguments.\n+cov.args <- CoverageVars(args.tbl, reg2plot)\n+attach(cov.args) #\n+plot.args <- PlotVars(args.tbl)\n+attach(plot.args)\n+\n+# Configuration: coverage-genelist-title relationships.\n+ctg.tbl <- ConfigTbl(args.tbl, fraglen)\n+\n+# Setup plot-related coordinates and variables.\n+plotvar.list <- SetupPlotCoord(args.tbl, ctg.tbl, default.tbl, dbfile.tbl, \n+ progpath, genome, reg2plot, inttag, flanksize, \n+ samprate, galaxy)\n+attach(plotvar.list)\n+\n+# Setup data points for plot.\n+pts.list <- SetPtsSpline(pint, lgint)\n+pts <- pts.list$pts # data points for avg. profile and standard errors.\n+m.pts <- pts.list$m.pts # middle data points. For pint, m.pts=1.\n+f.pts <- pts.list$f.pts # flanking region data points.\n+\n+# Setup matrix for avg. profiles.\n+regcovMat <- CreatePlotMat(pts, ctg.tbl)\n+# Setup matrix for standard errors.\n+confiMat <- CreateConfiMat(se, pts, ctg.tbl)\n+\n+# Genomic enrichment for all profiles in the config. Use this for heatmaps.\n+enrichList <- vector(\'list\', nrow(ctg.tbl))\n+cat("Done\\n")\n+\n+# Load required libraries.\n+cat("Loading R libraries")\n+if(!suppressMessages(require(ShortRead, warn.conflicts=F))) {\n+ source("http://bioconductor.org/biocLite.R")\n+ biocLite(ShortRead)\n+ if(!suppressMessages(require(ShortRead, warn.conflicts=F))) {\n+ '..b'####\n+# Add row names to heatmap data.\n+for(i in 1:length(enrichList)) {\n+ reg <- ctg.tbl$glist[i] # gene list name.\n+ rownames(enrichList[[i]]) <- paste(coord.list[[reg]]$gname, \n+ coord.list[[reg]]$tid, sep=\':\')\n+}\n+# Some basic parameters.\n+xticks <- genXticks(reg2plot, pint, lgint, pts, flanksize, flankfactor, Labs)\n+unit.width <- 4\n+ng.list <- sapply(enrichList, nrow) # number of genes per heatmap.\n+\n+# Create image file and plot data into it.\n+if(!fi_tag){\n+ cat("Plotting figures...")\n+ #### Average profile plot. ####\n+ out.plot <- avgname\n+ pdf(out.plot, width=plot.width, height=plot.height, pointsize=font.size)\n+ plotmat(regcovMat, ctg.tbl$title, ctg.tbl$color, bam.pair, xticks, pts, \n+ m.pts, f.pts, pint, shade.alp, confiMat, mw, prof.misc)\n+ out.dev <- dev.off()\n+\n+ #### Heatmap. ####\n+ # Setup output device.\n+ hd <- SetupHeatmapDevice(reg.list, uniq.reg, ng.list, pts, font.size, \n+ unit.width, rr)\n+ reg.hei <- hd$reg.hei # list of image heights for unique regions.\n+ hm.width <- hd$hm.width # image width.\n+ hm.height <- hd$hm.height # image height.\n+ lay.mat <- hd$lay.mat # matrix for heatmap layout.\n+ heatmap.mar <- hd$heatmap.mar # heatmap margins in inches.\n+\n+ out.hm <- heatmapname\n+ pdf(out.hm, width=hm.width, height=hm.height, pointsize=font.size)\n+ par(mai=heatmap.mar)\n+ layout(lay.mat, heights=reg.hei)\n+\n+ # Do heatmap plotting.\n+ go.list <- plotheat(reg.list, uniq.reg, enrichList, v.low.cutoff, go.algo, \n+ go.paras, ctg.tbl$title, bam.pair, xticks, flood.frac, \n+ do.plot=T, hm.color=hm.color, color.distr=color.distr, \n+ color.scale=color.scale)\n+ out.dev <- dev.off()\n+ cat("Done\\n")\n+} else {\n+ go.list <- plotheat(reg.list, uniq.reg, enrichList, v.low.cutoff, go.algo, \n+ go.paras, ctg.tbl$title, bam.pair, xticks, flood.frac, \n+ do.plot=F, hm.color=hm.color, color.distr=color.distr, \n+ color.scale=color.scale)\n+}\n+\n+# Save plotting data.\n+oname1="data"\n+cat("Saving results...")\n+dir.create(oname1, showWarnings=F)\n+\n+# Average profiles.\n+out.prof <- file.path(oname1, \'avgprof.txt\')\n+write.table(regcovMat, file=out.prof, row.names=F, sep="\\t", quote=F)\n+\n+# Standard errors of mean.\n+if(!is.null(confiMat)){\n+ out.confi <- file.path(oname1, \'sem.txt\')\n+ write.table(confiMat, file=out.confi, row.names=F, sep="\\t", quote=F)\n+}\n+\n+# Heatmap density values.\n+for(i in 1:length(enrichList)) {\n+ reg <- ctg.tbl$glist[i] # gene list name.\n+ out.heat <- file.path(oname1, paste(\'hm\', i, \'.txt\', sep=\'\'))\n+ write.table(cbind(coord.list[[reg]][, c(\'gid\', \'gname\', \'tid\', \'strand\')], \n+ enrichList[[i]]), \n+ file=out.heat, row.names=F, sep="\\t", quote=F)\n+}\n+\n+# Avg. profile R data.\n+prof.dat <- file.path(oname1, \'avgprof.RData\')\n+save(plot.width, plot.height, regcovMat, ctg.tbl, bam.pair, xticks, pts, \n+ m.pts, f.pts, pint, shade.alp, confiMat, mw, prof.misc, se, v.lib.size, \n+ font.size,\n+ file=prof.dat)\n+\n+# Heatmap R data.\n+heat.dat <- file.path(oname1, \'heatmap.RData\')\n+save(reg.list, uniq.reg, ng.list, pts, enrichList, v.low.cutoff, go.algo, \n+ ctg.tbl, bam.pair, xticks, flood.frac, hm.color, unit.width, rr, \n+ go.list, color.scale, v.lib.size, font.size, go.paras, low.count,\n+ color.distr, \n+ file=heat.dat)\n+cat("Done\\n")\n+\n+# Wrap results up.\n+cat("Wrapping results up...")\n+cur.dir <- getwd()\n+out.dir <- dirname(oname1)\n+out.zip <- basename(oname1)\n+setwd(out.dir)\n+if(!zip(paste(out.zip, \'.zip\', sep=\'\'), out.zip, extras=\'-q\')) {\n+ if(unlink(oname, recursive=T)) {\n+ warning(sprintf("Unable to delete intermediate result folder: %s", \n+ oname))\n+ }\n+}\n+setwd(cur.dir)\n+cat("Done\\n")\n+cat("All done. Cheers!\\n")\n+\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c ngsplot.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ngsplot.xml Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,3185 @@\n+<tool id="ngsplot" name="ngsplot" version="0.1.0">\n+ <description>visualize NGS samples at functional genomic regions - beta</description>\n+ <requirements>\n+ <requirement type="package" version="1.17.1=r3.3.2_1">r-catools</requirement>\n+ <requirement type="package" version="1.3.4">r-domc</requirement>\n+ <requirement type="package" version="1.42.0">bioconductor-bsgenome</requirement>\n+<!-- <requirement type="package" version="1.26.1">bioconductor-rsamtools</requirement> -->\n+ <requirement type="package" version="1.32.0">bioconductor-shortread</requirement>\n+<!-- <requirement type="package" version="0.16.1=r3.2.2_0">bioconductor-biocgenerics</requirement> -->\n+ </requirements>\n+\n+ <command><![CDATA[\n+export NGSPLOT="$__tool_directory__" &&\n+#if $numsamples.numsamples2 == "1":\n+ perl "$__tool_directory__"/runNGSplot.pl \n+ $genome_name \n+ $genomic_region_source_type.genomic_region \n+ $genomic_region_source_type.further_information \n+ $genomic_region_source_type.interval_size \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_option \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_size \n+ $numsamples.numsamples2\n+ $numsamples.usepair.usepair1\n+\n+ $numsamples.usepair.bamfile1\n+ $numsamples.usepair.reffile1\n+ $numsamples.usepair.genelist1.usegenelist1\n+ $numsamples.usepair.genelist1.genelist1\n+ $numsamples.usepair.title1\n+ $numsamples.usepair.fraglen1\n+ $numsamples.usepair.linecol1\n+\n+ $gene_database\n+ $randomly_sample\n+ $GO.gene_order\n+ $GO.KNC\n+ $GO.MIT\n+ $GO.NRS\n+ $chunk_size\n+ $quality_requirement\n+ $standard_error\n+ $radius_size\n+ $flooding_fraction\n+ $smooth_method\n+ $shaded_area\n+ $out_zip $out_avg_png $out_hm_png\n+\n+#else if $numsamples.numsamples2 == "2":\n+ perl "$__tool_directory__"/runNGSplot.pl \n+ $genome_name \n+ $genomic_region_source_type.genomic_region \n+ $genomic_region_source_type.further_information \n+ $genomic_region_source_type.interval_size \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_option \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_size \n+ $numsamples.numsamples2\n+ $numsamples.usepair.usepair1\n+\n+ $numsamples.usepair.bamfile1\n+ $numsamples.usepair.reffile1\n+ $numsamples.usepair.genelist1.usegenelist1\n+ $numsamples.usepair.genelist1.genelist1\n+ $numsamples.usepair.title1\n+ $numsamples.usepair.fraglen1\n+ $numsamples.usepair.linecol1\n+\n+ $numsamples.usepair.bamfile2\n+ $numsamples.usepair.reffile2\n+ $numsamples.usepair.genelist2.usegenelist2\n+ $numsamples.usepair.genelist2.genelist2\n+ $numsamples.usepair.title2\n+ $numsamples.usepair.fraglen2\n+ $numsamples.usepair.linecol2\n+\n+ $gene_database\n+ $randomly_sample\n+ $GO.gene_order\n+ $GO.KNC\n+ $GO.MIT\n+ $GO.NRS\n+ $chunk_size\n+ $quality_requirement\n+ $standard_error\n+ $radius_size\n+ $flooding_fraction\n+ $smooth_method\n+ $shaded_area\n+ $out_zip $out_avg_png $out_hm_png\n+\n+#else if $numsamples.numsamples2 == "3":\n+ perl "$__tool_directory__"/runNGSplot.pl \n+ $genome_name \n+ $genomic_region_source_type.genomic_region \n+ $genomic_region_source_type.further_information \n+ $genomic_region_source_type.interval_size \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_option \n+ $genomic_region_source_type.flanking_region_option_source_type.flanking_region_size \n+ $numsamples.numsamples2\n+ $numsamples.usepair.usepair1\n+\n+ $numsamples.usepair.bamfile1\n+ $numsamples.usepair.reffile1\n+ $numsamples.usepair.genelist1.usegenelist1\n+ $numsamples.usepair.genelist1.genelist1\n+ $numsamples.usepair.title1\n+ $numsamples.usepair.fraglen1\n+ $numsamples.usepair.linecol1\n+\n+ $numsamples.usepair.bamfile2\n+ $numsamples.usepair.reffile2\n+ $numsamples.usepair.genelist2.usegenelist2\n+ $n'..b'eak value of the 1st profile</option>\n+ <option value="prod">Product of all profiles on the same region</option>\n+ <option value="diff">Difference between the 1st and 2nd profiles</option>\n+ <option value="km">K-means clustering</option>\n+ <option value="none">No ranking algorithm applied. Use order in gene list.</option>\n+ </param>\n+ <when value="km">\n+ <param type="text" name="KNC" value="5" label="K-means: Number of clusters" help=""/>\n+ <param type="text" name="MIT" value="20" label="K-means: Number of iterations" help=""/>\n+ <param type="text" name="NRS" value="30" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when>\n+ <when value="total">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ <when value="hc">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ <when value="max">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ <when value="prod">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ <when value="diff">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ <when value="none">\n+ <param type="hidden" name="KNC" value="NA" />\n+ <param type="hidden" name="MIT" value="NA" />\n+ <param type="hidden" name="NRS" value="NA" />\n+ </when>\n+ </conditional>\n+ <param name="chunk_size" size="10" type="text" value="100" label="Chunk size for loading genes in batch" help="This parameter controls the behavior of coverage calculation. A smaller value implies lower memory footprint but may increase processing time."/>\n+ <param name="quality_requirement" size="10" type="text" value="20" label="Mapping quality requirement" help="This is the Phred-scale mapping quality score. A score of 20 means an error rate of 1%."/>\n+ <param type="select" name="standard_error" label="Plot standard errors" help="Standard errors will be rendered as shaded area around each curve.">\n+ <option value="1">Yes</option>\n+ <option value="0">No</option>\n+ </param>\n+\n+ <param name="radius_size" size="10" type="text" value="0" label="Fraction of extreme values to be trimmed on both ends" help="The fraction of extreme values that will be trimmed on both ends. Eg. 0.05 will remove 5% of extreme values."/>\n+ <param name="flooding_fraction" size="10" type="text" value="0.02" label="Heatmap flooding fraction" help="Default of 0.02 means that the minimum value is truncated at 2% and the maximum value is truncated at 98%. A higher fraction results in plots that have higher brightness but are less dynamic."/>\n+ <param type="select" name="smooth_method" label="Moving window width to smooth avg. profiles">\n+ <option value="1">No</option>\n+ <option value="3">Slightly</option>\n+ <option value="5">Somewhat</option>\n+ <option value="9">Quite</option>\n+ <option value="13">Super</option>\n+ </param>\n+ <param name="shaded_area" size="10" type="text" value="0" label="Opacity of shaded area under curve" help="Suggested value between 0 and 0.5."/>\n+<!--\n+-->\n+</inputs>\n+\n+<help></help>\n+<outputs>\n+ <data format="pdf" name="out_avg_png" />\n+ <data format="pdf" name="out_hm_png" />\n+ <data format="zip" name="out_zip"/>\n+</outputs>\n+</tool>\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c replot.xml.notworking --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/replot.xml.notworking Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| b'@@ -0,0 +1,157 @@\n+<tool id="replot" name="replot">\n+<command interpreter="perl">\n+#if $outtype.outtype2 == "prof":\n+ runreplot.pl $input_zipfile $outtype.outtype2 $output_filename \n+ $outtype.WD\n+ $outtype.HG\n+ $outtype.SE\n+ $outtype.H\n+ $outtype.MW\n+ $outtype.Y.Y2\n+ $outtype.Y.Y3\n+ $outtype.LEG\n+ $outtype.BOX\n+ $outtype.VLN\n+ $outtype.XYL\n+ $outtype.LWD\n+ $outtype.FS\n+#else if $outtype.outtype2 == "heatmap":\n+ runreplot.pl $input_zipfile $outtype.outtype2 $output_filename \n+ $outtype.GO.GO2\n+ $outtype.GO.KNC\n+ $outtype.GO.MIT\n+ $outtype.GO.NRS\n+ $outtype.LOW\n+ $outtype.RR\n+ $outtype.FC\n+ $outtype.CO\n+ $outtype.CD\n+ $outtype.FS\n+ dummy\n+ dummy\n+ dummy\n+#end if:\n+\n+</command>\n+\n+<!-->\n+\n+\n+<-->\n+\n+<inputs>\n+ <param type="data" name="input_zipfile" multiple="false" label="Input zip file created by ngsplot"/>\n+\n+ <conditional name="outtype">\n+ <param type="select" name="outtype2" label="Select figure type to replot">\n+ <option value="prof">Coverage Profile</option>\n+ <option value="heatmap">Heatmap</option>\n+ </param>\n+\n+ <when value="prof">\n+ <param type="text" name="WD" value="8" label="Image width in inches" />\n+ <param type="text" name="HG" value="7" label="Image height in inches" />\n+ <param type="select" name="SE" label="Standard error shading" >\n+ <option value="1">Yes, include standard error shading</option>\n+ <option value="0">No, exclude standard error shading</option>\n+ </param>\n+ <param type="text" name="H" value="0" label="Opacity of shaded area under each curve" help="Use 0 for no shading. Suggested value of between 0 and 0.5 will add a semi-transparent shade under each curve."/>\n+ <param type="text" name="MW" value="1" label="Moving window width to smooth average profiles" help="The unit of the window size is a data point. In ngs.plot, the number of data points on the X-axis is 100. A moving window width of 10 means averaging over 10% of the entire X-axis."/>\n+ <conditional name="Y">\n+ <param type="select" name="Y2" label="Define Y-axis scale?" >\n+ <option value="no">No, let program decide</option>\n+ <option value="yes">Yes</option>\n+ </param>\n+ <when value="yes">\n+ <param type="text" name="Y3" value="0,100" label="Y-axis min_val,max_val" help="Enter minimum and maximum value separated by comma" />\n+ </when>\n+ <when value="no">\n+ <param type="hidden" name="Y3" value="na" label="test"/>\n+ </when>\n+ </conditional>\n+ <param type="select" name="LEG" label="Draw legend?" >\n+ <option value="1">Yes</option>\n+ <option value="0">No</option>\n+ </param>\n+ <param type="select" name="BOX" label="Draw box around plot?" >\n+ <option value="1">Yes</option>\n+ <option value="0">No</option>\n+ </param>\n+ <param type="select" name="VLN" label="Draw vertical lines?" >\n+ <option value="1">Yes</option>\n+ <option value="0">No</option>\n+ </param>\n+ <param type="select" name="XYL" label="Draw axis labels?" >\n+ <option value="1">Yes</option>\n+ <option value="0">No</option>\n+ </param>\n+ <param type="text" name="LWD" value="3" label="Line width" />\n+ <param type="text" name="FS" value="12" label="Font size" />\n+<!-->\n+<-->\n+\n+ </when>\n+ <when value="heatmap">\n+ <conditional name="GO">\n+ <param type="select" name="GO2" label="Gene order algorithm" >\n+ <option value="total">Overall enrichment of the 1st profile (default)</option>\n+ <option value="hc">Hierarchial clustering</option>\n+ <option value="max">Peak value of the 1st profile</option>\n+ <option value="prod">Product of all profiles on the same '..b'ptima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ <when value="hc">\n+ <param type="hidden" name="KNC" value="NA" label="K-means: Number of clusters" help=""/>\n+ <param type="hidden" name="MIT" value="NA" label="K-means: Number of iterations" help=""/>\n+ <param type="hidden" name="NRS" value="NA" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ <when value="max">\n+ <param type="hidden" name="KNC" value="NA" label="K-means: Number of clusters" help=""/>\n+ <param type="hidden" name="MIT" value="NA" label="K-means: Number of iterations" help=""/>\n+ <param type="hidden" name="NRS" value="NA" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ <when value="prod">\n+ <param type="hidden" name="KNC" value="NA" label="K-means: Number of clusters" help=""/>\n+ <param type="hidden" name="MIT" value="NA" label="K-means: Number of iterations" help=""/>\n+ <param type="hidden" name="NRS" value="NA" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ <when value="diff">\n+ <param type="hidden" name="KNC" value="NA" label="K-means: Number of clusters" help=""/>\n+ <param type="hidden" name="MIT" value="NA" label="K-means: Number of iterations" help=""/>\n+ <param type="hidden" name="NRS" value="NA" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ <when value="none">\n+ <param type="hidden" name="KNC" value="NA" label="K-means: Number of clusters" help=""/>\n+ <param type="hidden" name="MIT" value="NA" label="K-means: Number of iterations" help=""/>\n+ <param type="hidden" name="NRS" value="NA" label="K-means: Number of random starts" help="K-means is prone to local optima. Restarting it repeatedly may help to find a better solution."/>\n+ </when> \n+ </conditional>\n+ <param type="text" name="LOW" value="10" label="Low count cutoff in rank-based normalization" help="This number is in raw read counts."/>\n+ <param type="text" name="RR" value="30" label="Reduce ratio" help="This controls the heatmap height. The smaller the value, the taller the heatmap." />\n+ <param type="text" name="FC" value="0.02" label="Flooding fraction" help="The default value of 0.02 means that the minimum value is truncated at 2% and the maximum value is truncated at 98%. A higher fraction results in plots that have higher brightness but are less dynamic." />\n+ <param type="text" name="CO" value="default" label="Color for heatmap" help="For non bam-pair plots, input a single color. For bam-pair plots, use color-tri format (neg_color:[neu_color]:pos_color). Hint: must use R colors, such as darkgreen, yellow and blue2. The neutral color is optional(default=black)." />\n+ <param type="text" name="CD" value="0.6" label="Color distribution for heatmap" help="Must be a positive number. Hint: lower values give more widely spaced colors at the negative end. In other words, they shift the neutral color to positive values. If set to 1, the neutral color represents 0(i.e. no bias). If set to >1, the neutral color represents a negative value." />\n+ <param type="text" name="FS" value="12" label="Font size" />\n+ </when>\n+ </conditional>\n+<!--\n+-->\n+</inputs>\n+\n+<help></help>\n+<outputs>\n+ <data format="pdf" name="output_filename" />\n+</outputs>\n+</tool>\n' |
| b |
| diff -r 000000000000 -r 3ca58369469c runNGSplot.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/runNGSplot.pl Wed Dec 06 19:01:53 2017 -0500 |
| b |
| @@ -0,0 +1,129 @@ +#!/usr/bin/perl -w +use strict; +use File::Basename 'dirname'; +use File::Spec; +use Cwd 'abs_path'; + +my @inputs = @ARGV; +my @inputs2 = @inputs; + +my $genome_name = shift(@inputs); +my $genomic_region_source_type__genomic_region = shift(@inputs); +my $genomic_region_source_type__further_information = shift(@inputs); +my $genomic_region_source_type__interval_size = shift(@inputs); +my $genomic_region_source_type__flanking_region_option_source_type__flanking_region_option = shift(@inputs); +my $genomic_region_source_type__flanking_region_option_source_type__flanking_region_size = shift(@inputs); +my $numsamples = shift(@inputs); +my $usepairs = shift(@inputs); + +#print STDERR "inputs @inputs\n"; + +my $randfile = rand(100)."\.config\.txt"; +my $randfile2 = $randfile; +$randfile2 =~ s/config\.txt/logfile/gm; + +my $outfile = File::Spec->catfile(abs_path(dirname(__FILE__)),"$randfile"); +open(FILE,">$outfile"); + +for (my $i=1;$i<=$numsamples;$i++) { + my $bamfile=shift(@inputs); + my $reffile=shift(@inputs); + my $usegenelist=shift(@inputs); + my $genelist=shift(@inputs); + my $title=shift(@inputs); + my $fraglen=shift(@inputs); + my $color=shift(@inputs); + if ($usepairs eq 'yes') { + syswrite(FILE, "$bamfile\:$reffile\t$genelist\t$title\t$fraglen\t$color\n"); + }else { + syswrite(FILE, "$bamfile\t$genelist\t$title\t$fraglen\t$color\n"); + } +} +close(FILE); + +my $gene_database = shift(@inputs); +my $randomly_sample = shift(@inputs); +my $gene_order = shift(@inputs); +my $knc = shift(@inputs); +my $mit = shift(@inputs); +my $nrs = shift(@inputs); +my $chunk_size = shift(@inputs); +my $quality_requirement = shift(@inputs); +my $standard_error = shift(@inputs); +my $radius_size = shift(@inputs); +my $flooding_fraction = shift(@inputs); +my $smooth_method = shift(@inputs); +my $shaded_area = shift(@inputs); +my $out_name = shift(@inputs); +my $out_avg_name = shift(@inputs); +my $out_hm_name = shift(@inputs); + + +my $G = $genome_name; +my $R = $genomic_region_source_type__genomic_region; +my $C = $outfile; +my $O = $out_name; +my $O2 = $out_avg_name; +my $O3 = $out_hm_name; +my $F = $genomic_region_source_type__further_information; +my $D = $gene_database; +my $L = $genomic_region_source_type__flanking_region_option_source_type__flanking_region_size; +my $N = $genomic_region_source_type__flanking_region_option_source_type__flanking_region_size; +my $RB = $radius_size; +my $S = $randomly_sample; +my $CS = $chunk_size; +my $MQ = $quality_requirement; +my $IN = $genomic_region_source_type__interval_size; +my $SE = $standard_error; +my $MW = $smooth_method; +my $H = $shaded_area; +my $GO = $gene_order; +my $KNC = $knc; +my $MIT = $mit; +my $NRS = $nrs; +my $FC = $flooding_fraction; + +if ($GO eq 'km') { + $GO = "$GO -KNC $KNC -MIT $MIT -NRS $NRS"; +} + +my $logfile = File::Spec->catfile(abs_path(dirname(__FILE__)),"$randfile2"); +open(FILE2,">>$logfile"); +my $cmd5="pwd >> $logfile 2>&1"; +system($cmd5); +my $NGSPLOT = $ENV{"NGSPLOT"}; +my $cmd=''; +if (($R eq 'tss')||($R eq 'tes')) { + $cmd = "Rscript $NGSPLOT/ngs.plot.r -Galaxy 1 -P 0 -G $G -R $R -C $C -O $O -O2 $O2 -O3 $O3 -D $D -L $L -S $S -GO $GO -CS $CS -MQ $MQ -SE $SE -RB $RB -FC $FC -MW $MW -H $H >> $logfile 2>&1"; +}elsif ($genomic_region_source_type__flanking_region_option_source_type__flanking_region_option eq 'flanking_region_size') { + if ($IN eq 'automatic') { + $cmd = "$NGSPLOT/ngs.plot.r -Galaxy 1 -P 0 -G $G -R $R -C $C -O $O -O2 $O2 -O3 $O3 -D $D -L $L -S $S -GO $GO -CS $CS -MQ $MQ -SE $SE -RB $RB -FC $FC -MW $MW -H $H >> $logfile 2>&1"; + }else { + $cmd = "$NGSPLOT/ngs.plot.r -Galaxy 1 -P 0 -G $G -R $R -C $C -O $O -O2 $O2 -O3 $O3 -D $D -L $L -S $S -GO $GO -CS $CS -MQ $MQ -SE $SE -RB $RB -FC $FC -MW $MW -H $H -IN $IN >> $logfile 2>&1"; + } +}elsif ($genomic_region_source_type__flanking_region_option_source_type__flanking_region_option eq 'flanking_floating_size') { + if ($IN eq 'automatic') { + $cmd = "$NGSPLOT/ngs.plot.r -Galaxy 1 -P 0 -G $G -R $R -C $C -O $O -O2 $O2 -O3 $O3 -D $D -N $N -S $S -GO $GO -CS $CS -MQ $MQ -SE $SE -RB $RB -FC $FC -MW $MW -H $H >> $logfile 2>&1"; + }else { + $cmd = "$NGSPLOT/ngs.plot.r -Galaxy 1 -P 0 -G $G -R $R -C $C -O $O -O2 $O2 -O3 $O3 -D $D -N $N -S $S -GO $GO -CS $CS -MQ $MQ -SE $SE -RB $RB -FC $FC -MW $MW -H $H -IN $IN >> $logfile 2>&1"; + } +} + +print("$cmd\n"); +my $cmd2="cp data.zip $O"; +my $cmd3="rm $outfile"; +my $cmd4="rm $logfile"; +syswrite(FILE2, "\n$cmd\n"); +syswrite(FILE2, "\n$cmd2\n"); +syswrite(FILE2, "\n$cmd3\n"); +syswrite(FILE2, "\n$cmd4\n\n"); + +# print STDERR "cmd: $cmd\n"; +system($cmd); +system($cmd2); +#system($cmd3); +#system($cmd4); + +close(FILE2); + + |
| b |
| diff -r 000000000000 -r 3ca58369469c runreplot.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/runreplot.pl Wed Dec 06 19:01:53 2017 -0500 |
| [ |
| @@ -0,0 +1,63 @@ +#!/usr/bin/perl -w +use strict; +use File::Basename 'dirname'; +use File::Spec; +use Cwd 'abs_path'; + +my $input_file = $ARGV[0]; +my $output_type = $ARGV[1]; +my $output_filename = $ARGV[2]; + +my $cmd1 = "cp $input_file data.zip"; +my $cmd2 = ""; +my $cmd3 = "mv output.pdf $output_filename"; + +my $FS = ""; +if ($output_type eq 'prof') { + my $WD = "-WD $ARGV[3]"; + my $HG = "-HG $ARGV[4]"; + my $SE = "-SE $ARGV[5]"; + my $H = "-H $ARGV[6]"; + my $MW = "-MW $ARGV[7]"; + my $Y = "-Y $ARGV[8]"; + my $YAS = "-YAS $ARGV[9]"; + my $LEG = "-LEG $ARGV[10]"; + my $BOX = "-BOX $ARGV[11]"; + my $VLN = "-VLN $ARGV[12]"; + my $XYL = "-XYL $ARGV[13]"; + my $LWD = "-LWD $ARGV[14]"; + $FS = "-FS $ARGV[15]"; + + if ($Y =~ /no$/gm) { + $cmd2 = "replot.r $output_type -I data.zip -O output $WD $HG $SE $H $MW $LEG $BOX $VLN $XYL $LWD $FS"; + }else { + $cmd2 = "replot.r $output_type -I data.zip -O output $WD $HG $SE $H $MW $LEG $BOX $VLN $XYL $LWD $FS $YAS"; + } +}elsif ($output_type eq 'heatmap') { + my $GO = "-GO $ARGV[3]"; + my $KNC = "-KNC $ARGV[4]"; + my $MIT = "-MIT $ARGV[5]"; + my $NRS = "-NRS $ARGV[6]"; + my $LOW = "-LOW $ARGV[7]"; + my $RR = "-RR $ARGV[8]"; + my $FC = "-FC $ARGV[9]"; + my $CO = "-CO $ARGV[10]"; + my $CD = "-CD $ARGV[11]"; + $FS = "-FS $ARGV[12]"; + my $DUM1 = $ARGV[13]; + my $DUM2 = $ARGV[14]; + my $DUM3 = $ARGV[15]; + + if ($CO =~ /default/gm) { $CO = ""; } + + if ($GO =~ /km$/gm) { + $cmd2 = "replot.r $output_type -I data.zip -O output $LOW $RR $FC $CO $CD $FS $GO $KNC $MIT $NRS"; + }else { + $cmd2 = "replot.r $output_type -I data.zip -O output $LOW $RR $FC $CO $CD $FS $GO"; + } +} +system($cmd1); +system($cmd2); +system($cmd3); + + |