Mercurial > repos > fubar > rglasso_1_9_8
changeset 3:e0e11c2cae3f draft
Uploaded
| author | fubar |
|---|---|
| date | Fri, 19 Dec 2014 05:57:40 -0500 |
| parents | 359c66be8677 |
| children | 5f2db639f8eb |
| files | readme.rst rgToolFactory.py rg_nri.xml rglasso_cox.xml test-data/cox_coxlassotest_glmnet_cvdeviance.pdf test-data/cox_coxlassotest_glmnetdev.pdf test-data/cox_test.xls test-data/coxlassotest.html test-data/coxlassotest_modelres.xls tool_dependencies.xml |
| diffstat | 10 files changed, 666 insertions(+), 1880 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/readme.rst Fri Dec 19 05:57:40 2014 -0500 @@ -0,0 +1,39 @@ +glmnet wrappers +=============== + +This is a self installing Galaxy tool exposing the glmnet_ R package which has excellent documentation at +glmnet_ Minimal details are provided in this wrapper - please RTM to get the best out of it. + +The tool exposes the entire range of penalised maximum likelihood +GLM models ranging from pure lasso (set alpha to 1) to pure ridge-regression (set alpha to 0). + +These models can be k-fold internally cross validated to help select an "optimal" predictive or classification +algorithm. Predictive coefficients for each included independent variable are output for each model. + +Predictors can be forced into models to adjust for known confounders or explanatory factors. + +The glmnet_ implementation of the coordinate descent algorithm is fast and efficient even on relatively large problems +with tens of thousands of predictors and thousands of samples - such as normalised microarray intensities and anthropometry +on a very large sample of obese patients. + +The user supplies a tabular file with rows as samples and columns containing observations, then chooses +as many predictors as required. A separate model will be output for each of potentially multiple dependent +variables. Models are reported as the coefficients for terms in an 'optimal' model. +These optimal predictors are selected by repeatedly setting +aside a random subsample, building a model in the remainder and estimating AUC or deviance +using k (default 10) fold internal cross validation. For each of these steps, a random 1/k +of the samples are set aside and used to estiamte performance of an optimal model estimated +from the remaining samples. Plots are provided showing the range of these (eg 10) internal validation +estimates and mean model AUC (binomial) or residual deviance plots at each penalty increment step. + +A full range of link functions are available including Gaussian, Poisson, Binomial and +Cox proportional hazard time to failure for censored data in this wrapper. + +Note that multinomial and multiresponse gaussian models are NOT yet implemented since I have not yet +had use for them - send code! + +.. _glmnet: http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html + +Wrapper author: Ross Lazarus +19 october 2014 +
--- a/rgToolFactory.py Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,658 +0,0 @@ -# rgToolFactory.py -# see https://bitbucket.org/fubar/galaxytoolfactory/wiki/Home -# -# copyright ross lazarus (ross stop lazarus at gmail stop com) May 2012 -# -# all rights reserved -# Licensed under the LGPL -# suggestions for improvement and bug fixes welcome at https://bitbucket.org/fubar/galaxytoolfactory/wiki/Home -# -# july 2014 -# added buffered read of sterror after run -# -# august 2013 -# found a problem with GS if $TMP or $TEMP missing - now inject /tmp and warn -# -# july 2013 -# added ability to combine images and individual log files into html output -# just make sure there's a log file foo.log and it will be output -# together with all images named like "foo_*.pdf -# otherwise old format for html -# -# January 2013 -# problem pointed out by Carlos Borroto -# added escaping for <>$ - thought I did that ages ago... -# -# August 11 2012 -# changed to use shell=False and cl as a sequence - -# This is a Galaxy tool factory for simple scripts in python, R or whatever ails ye. -# It also serves as the wrapper for the new tool. -# -# you paste and run your script -# Only works for simple scripts that read one input from the history. -# Optionally can write one new history dataset, -# and optionally collect any number of outputs into links on an autogenerated HTML page. - -# DO NOT install on a public or important site - please. - -# installed generated tools are fine if the script is safe. -# They just run normally and their user cannot do anything unusually insecure -# but please, practice safe toolshed. -# Read the fucking code before you install any tool -# especially this one - -# After you get the script working on some test data, you can -# optionally generate a toolshed compatible gzip file -# containing your script safely wrapped as an ordinary Galaxy script in your local toolshed for -# safe and largely automated installation in a production Galaxy. - -# If you opt for an HTML output, you get all the script outputs arranged -# as a single Html history item - all output files are linked, thumbnails for all the pdfs. -# Ugly but really inexpensive. -# -# Patches appreciated please. -# -# -# long route to June 2012 product -# Behold the awesome power of Galaxy and the toolshed with the tool factory to bind them -# derived from an integrated script model -# called rgBaseScriptWrapper.py -# Note to the unwary: -# This tool allows arbitrary scripting on your Galaxy as the Galaxy user -# There is nothing stopping a malicious user doing whatever they choose -# Extremely dangerous!! -# Totally insecure. So, trusted users only -# -# preferred model is a developer using their throw away workstation instance - ie a private site. -# no real risk. The universe_wsgi.ini admin_users string is checked - only admin users are permitted to run this tool. -# - -import sys -import shutil -import subprocess -import os -import time -import tempfile -import optparse -import tarfile -import re -import shutil -import math - -progname = os.path.split(sys.argv[0])[1] -myversion = 'V000.2 June 2012' -verbose = False -debug = False -toolFactoryURL = 'https://bitbucket.org/fubar/galaxytoolfactory' -buffsize = 1048576 - - -def timenow(): - """return current time as a string - """ - return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time())) - -html_escape_table = { - "&": "&", - ">": ">", - "<": "<", - "$": "\$" - } - -def html_escape(text): - """Produce entities within text.""" - return "".join(html_escape_table.get(c,c) for c in text) - -def cmd_exists(cmd): - return subprocess.call("type " + cmd, shell=True, - stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0 - - -class ScriptRunner: - """class is a wrapper for an arbitrary script - """ - - def __init__(self,opts=None,treatbashSpecial=True): - """ - cleanup inputs, setup some outputs - - """ - self.useGM = cmd_exists('gm') - self.useIM = cmd_exists('convert') - self.useGS = cmd_exists('gs') - self.temp_warned = False # we want only one warning if $TMP not set - self.treatbashSpecial = treatbashSpecial - if opts.output_dir: # simplify for the tool tarball - os.chdir(opts.output_dir) - self.thumbformat = 'png' - self.opts = opts - self.toolname = re.sub('[^a-zA-Z0-9_]+', '', opts.tool_name) # a sanitizer now does this but.. - self.toolid = self.toolname - self.myname = sys.argv[0] # get our name because we write ourselves out as a tool later - self.pyfile = self.myname # crude but efficient - the cruft won't hurt much - self.xmlfile = '%s.xml' % self.toolname - s = open(self.opts.script_path,'r').readlines() - s = [x.rstrip() for x in s] # remove pesky dos line endings if needed - self.script = '\n'.join(s) - fhandle,self.sfile = tempfile.mkstemp(prefix=self.toolname,suffix=".%s" % (opts.interpreter)) - tscript = open(self.sfile,'w') # use self.sfile as script source for Popen - tscript.write(self.script) - tscript.close() - self.indentedScript = '\n'.join([' %s' % x for x in s]) # for restructured text in help - self.escapedScript = '\n'.join([html_escape(x) for x in s]) - self.elog = os.path.join(self.opts.output_dir,"%s_error.log" % self.toolname) - self.tlog = os.path.join(self.opts.output_dir,"%s_runner.log" % self.toolname) - if opts.output_dir: # may not want these complexities - art = '%s.%s' % (self.toolname,opts.interpreter) - artpath = os.path.join(self.opts.output_dir,art) # need full path - artifact = open(artpath,'w') # use self.sfile as script source for Popen - artifact.write(self.script) - artifact.close() - self.cl = [] - self.html = [] - a = self.cl.append - a(opts.interpreter) - if self.treatbashSpecial and opts.interpreter in ['bash','sh']: - a(self.sfile) - else: - a('-') # stdin - a(opts.input_tab) - a(opts.output_tab) - self.outFormats = 'tabular' # TODO make this an option at tool generation time - self.inputFormats = 'tabular' # TODO make this an option at tool generation time - self.test1Input = '%s_test1_input.xls' % self.toolname - self.test1Output = '%s_test1_output.xls' % self.toolname - self.test1HTML = '%s_test1_output.html' % self.toolname - - def makeXML(self): - """ - Create a Galaxy xml tool wrapper for the new script as a string to write out - fixme - use templating or something less fugly than this example of what we produce - - <tool id="reverse" name="reverse" version="0.01"> - <description>a tabular file</description> - <command interpreter="python"> - reverse.py --script_path "$runMe" --interpreter "python" - --tool_name "reverse" --input_tab "$input1" --output_tab "$tab_file" - </command> - <inputs> - <param name="input1" type="data" format="tabular" label="Select a suitable input file from your history"/><param name="job_name" type="text" label="Supply a name for the outputs to remind you what they contain" value="reverse"/> - - </inputs> - <outputs> - <data format="tabular" name="tab_file" label="${job_name}"/> - - </outputs> - <help> - -**What it Does** - -Reverse the columns in a tabular file - - </help> - <configfiles> - <configfile name="runMe"> - -# reverse order of columns in a tabular file -import sys -inp = sys.argv[1] -outp = sys.argv[2] -i = open(inp,'r') -o = open(outp,'w') -for row in i: - rs = row.rstrip().split('\t') - rs.reverse() - o.write('\t'.join(rs)) - o.write('\n') -i.close() -o.close() - - - </configfile> - </configfiles> - </tool> - - """ - newXML="""<tool id="%(toolid)s" name="%(toolname)s" version="%(tool_version)s"> - %(tooldesc)s - %(command)s - <inputs> - %(inputs)s - </inputs> - <outputs> - %(outputs)s - </outputs> - <configfiles> - <configfile name="runMe"> - %(script)s - </configfile> - </configfiles> - %(tooltests)s - <help> - %(help)s - </help> - </tool>""" # needs a dict with toolname, toolid, interpreter, scriptname, command, inputs as a multi line string ready to write, outputs ditto, help ditto - - newCommand="""<command interpreter="python"> - %(toolname)s.py --script_path "$runMe" --interpreter "%(interpreter)s" - --tool_name "%(toolname)s" %(command_inputs)s %(command_outputs)s - </command>""" # may NOT be an input or htmlout - tooltestsTabOnly = """<tests><test> - <param name="input1" value="%(test1Input)s" ftype="tabular"/> - <param name="job_name" value="test1"/> - <param name="runMe" value="$runMe"/> - <output name="tab_file" file="%(test1Output)s" ftype="tabular"/> - </test></tests>""" - tooltestsHTMLOnly = """<tests><test> - <param name="input1" value="%(test1Input)s" ftype="tabular"/> - <param name="job_name" value="test1"/> - <param name="runMe" value="$runMe"/> - <output name="html_file" file="%(test1HTML)s" ftype="html" lines_diff="5"/> - </test></tests>""" - tooltestsBoth = """<tests><test> - <param name="input1" value="%(test1Input)s" ftype="tabular"/> - <param name="job_name" value="test1"/> - <param name="runMe" value="$runMe"/> - <output name="tab_file" file="%(test1Output)s" ftype="tabular" /> - <output name="html_file" file="%(test1HTML)s" ftype="html" lines_diff="10"/> - </test></tests>""" - xdict = {} - xdict['tool_version'] = self.opts.tool_version - xdict['test1Input'] = self.test1Input - xdict['test1HTML'] = self.test1HTML - xdict['test1Output'] = self.test1Output - if self.opts.make_HTML and self.opts.output_tab <> 'None': - xdict['tooltests'] = tooltestsBoth % xdict - elif self.opts.make_HTML: - xdict['tooltests'] = tooltestsHTMLOnly % xdict - else: - xdict['tooltests'] = tooltestsTabOnly % xdict - xdict['script'] = self.escapedScript - # configfile is least painful way to embed script to avoid external dependencies - # but requires escaping of <, > and $ to avoid Mako parsing - if self.opts.help_text: - xdict['help'] = open(self.opts.help_text,'r').read() - else: - xdict['help'] = 'Please ask the tool author for help as none was supplied at tool generation' - coda = ['**Script**','Pressing execute will run the following code over your input file and generate some outputs in your history::'] - coda.append(self.indentedScript) - coda.append('**Attribution** This Galaxy tool was created by %s at %s\nusing the Galaxy Tool Factory.' % (self.opts.user_email,timenow())) - coda.append('See %s for details of that project' % (toolFactoryURL)) - coda.append('Please cite: Creating re-usable tools from scripts: The Galaxy Tool Factory. Ross Lazarus; Antony Kaspi; Mark Ziemann; The Galaxy Team. ') - coda.append('Bioinformatics 2012; doi: 10.1093/bioinformatics/bts573') - xdict['help'] = '%s\n%s' % (xdict['help'],'\n'.join(coda)) - if self.opts.tool_desc: - xdict['tooldesc'] = '<description>%s</description>' % self.opts.tool_desc - else: - xdict['tooldesc'] = '' - xdict['command_outputs'] = '' - xdict['outputs'] = '' - if self.opts.input_tab <> 'None': - xdict['command_inputs'] = '--input_tab "$input1" ' # the space may matter a lot if we append something - xdict['inputs'] = '<param name="input1" type="data" format="%s" label="Select a suitable input file from your history"/> \n' % self.inputFormats - else: - xdict['command_inputs'] = '' # assume no input - eg a random data generator - xdict['inputs'] = '' - xdict['inputs'] += '<param name="job_name" type="text" label="Supply a name for the outputs to remind you what they contain" value="%s"/> \n' % self.toolname - xdict['toolname'] = self.toolname - xdict['toolid'] = self.toolid - xdict['interpreter'] = self.opts.interpreter - xdict['scriptname'] = self.sfile - if self.opts.make_HTML: - xdict['command_outputs'] += ' --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" ' - xdict['outputs'] += ' <data format="html" name="html_file" label="${job_name}.html"/>\n' - if self.opts.output_tab <> 'None': - xdict['command_outputs'] += ' --output_tab "$tab_file"' - xdict['outputs'] += ' <data format="%s" name="tab_file" label="${job_name}"/>\n' % self.outFormats - xdict['command'] = newCommand % xdict - xmls = newXML % xdict - xf = open(self.xmlfile,'w') - xf.write(xmls) - xf.write('\n') - xf.close() - # ready for the tarball - - - def makeTooltar(self): - """ - a tool is a gz tarball with eg - /toolname/tool.xml /toolname/tool.py /toolname/test-data/test1_in.foo ... - """ - retval = self.run() - if retval: - print >> sys.stderr,'## Run failed. Cannot build yet. Please fix and retry' - sys.exit(1) - self.makeXML() - tdir = self.toolname - os.mkdir(tdir) - if self.opts.input_tab <> 'None': # no reproducible test otherwise? TODO: maybe.. - testdir = os.path.join(tdir,'test-data') - os.mkdir(testdir) # make tests directory - shutil.copyfile(self.opts.input_tab,os.path.join(testdir,self.test1Input)) - if self.opts.output_tab <> 'None': - shutil.copyfile(self.opts.output_tab,os.path.join(testdir,self.test1Output)) - if self.opts.make_HTML: - shutil.copyfile(self.opts.output_html,os.path.join(testdir,self.test1HTML)) - if self.opts.output_dir: - shutil.copyfile(self.tlog,os.path.join(testdir,'test1_out.log')) - op = '%s.py' % self.toolname # new name - outpiname = os.path.join(tdir,op) # path for the tool tarball - pyin = os.path.basename(self.pyfile) # our name - we rewrite ourselves (TM) - notes = ['# %s - a self annotated version of %s generated by running %s\n' % (op,pyin,pyin),] - notes.append('# to make a new Galaxy tool called %s\n' % self.toolname) - notes.append('# User %s at %s\n' % (self.opts.user_email,timenow())) - pi = open(self.pyfile,'r').readlines() # our code becomes new tool wrapper (!) - first Galaxy worm - notes += pi - outpi = open(outpiname,'w') - outpi.write(''.join(notes)) - outpi.write('\n') - outpi.close() - stname = os.path.join(tdir,self.sfile) - if not os.path.exists(stname): - shutil.copyfile(self.sfile, stname) - xtname = os.path.join(tdir,self.xmlfile) - if not os.path.exists(xtname): - shutil.copyfile(self.xmlfile,xtname) - tarpath = "%s.gz" % self.toolname - tar = tarfile.open(tarpath, "w:gz") - tar.add(tdir,arcname=self.toolname) - tar.close() - shutil.copyfile(tarpath,self.opts.new_tool) - shutil.rmtree(tdir) - ## TODO: replace with optional direct upload to local toolshed? - return retval - - - def compressPDF(self,inpdf=None,thumbformat='png'): - """need absolute path to pdf - note that GS gets confoozled if no $TMP or $TEMP - so we set it - """ - assert os.path.isfile(inpdf), "## Input %s supplied to %s compressPDF not found" % (inpdf,self.myName) - our_env = os.environ.copy() - if not (our_env.get('TMP',None) or our_env.get('TEMP',None)): - our_env['TMP'] = '/tmp' - if not self.temp_warned: - print >> sys.stdout,'## WARNING - no $TMP or $TEMP!!! Please fix - using /tmp temporarily' - self.temp_warned = True - hlog = os.path.join(self.opts.output_dir,"compress_%s.txt" % os.path.basename(inpdf)) - sto = open(hlog,'w') - outpdf = '%s_compressed' % inpdf - cl = ["gs", "-sDEVICE=pdfwrite", "-dNOPAUSE", "-dUseCIEColor", "-dBATCH","-dPDFSETTINGS=/printer", "-sOutputFile=%s" % outpdf,inpdf] - x = subprocess.Popen(cl,stdout=sto,stderr=sto,cwd=self.opts.output_dir,env=our_env) - retval1 = x.wait() - sto.close() - if retval1 == 0: - os.unlink(inpdf) - shutil.move(outpdf,inpdf) - os.unlink(hlog) - else: - x = open(hlog,'r').readlines() - print >> sys.stdout,x - hlog = os.path.join(self.opts.output_dir,"thumbnail_%s.txt" % os.path.basename(inpdf)) - sto = open(hlog,'w') - outpng = '%s.%s' % (os.path.splitext(inpdf)[0],thumbformat) - if self.useGM: - cl2 = ['gm', 'convert', inpdf, outpng] - else: # assume imagemagick - cl2 = ['convert', inpdf, outpng] - x = subprocess.Popen(cl2,stdout=sto,stderr=sto,cwd=self.opts.output_dir,env=our_env) - retval2 = x.wait() - sto.close() - if retval2 <> 0: - x = open(hlog,'r').readlines() - print >> sys.stdout,x - else: - os.unlink(hlog) - retval = retval1 or retval2 - return retval - - - def getfSize(self,fpath,outpath): - """ - format a nice file size string - """ - size = '' - fp = os.path.join(outpath,fpath) - if os.path.isfile(fp): - size = '0 B' - n = float(os.path.getsize(fp)) - if n > 2**20: - size = '%1.1f MB' % (n/2**20) - elif n > 2**10: - size = '%1.1f KB' % (n/2**10) - elif n > 0: - size = '%d B' % (int(n)) - return size - - def makeHtml(self): - """ Create an HTML file content to list all the artifacts found in the output_dir - """ - - galhtmlprefix = """<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> - <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"> - <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> - <meta name="generator" content="Galaxy %s tool output - see http://g2.trac.bx.psu.edu/" /> - <title></title> - <link rel="stylesheet" href="/static/style/base.css" type="text/css" /> - </head> - <body> - <div class="toolFormBody"> - """ - galhtmlattr = """<hr/><div class="infomessage">This tool (%s) was generated by the <a href="https://bitbucket.org/fubar/galaxytoolfactory/overview">Galaxy Tool Factory</a></div><br/>""" - galhtmlpostfix = """</div></body></html>\n""" - - flist = os.listdir(self.opts.output_dir) - flist = [x for x in flist if x <> 'Rplots.pdf'] - flist.sort() - html = [] - html.append(galhtmlprefix % progname) - html.append('<div class="infomessage">Galaxy Tool "%s" run at %s</div><br/>' % (self.toolname,timenow())) - fhtml = [] - if len(flist) > 0: - logfiles = [x for x in flist if x.lower().endswith('.log')] # log file names determine sections - logfiles.sort() - logfiles = [x for x in logfiles if os.path.abspath(x) <> os.path.abspath(self.tlog)] - logfiles.append(os.path.abspath(self.tlog)) # make it the last one - pdflist = [] - npdf = len([x for x in flist if os.path.splitext(x)[-1].lower() == '.pdf']) - for rownum,fname in enumerate(flist): - dname,e = os.path.splitext(fname) - sfsize = self.getfSize(fname,self.opts.output_dir) - if e.lower() == '.pdf' : # compress and make a thumbnail - thumb = '%s.%s' % (dname,self.thumbformat) - pdff = os.path.join(self.opts.output_dir,fname) - retval = self.compressPDF(inpdf=pdff,thumbformat=self.thumbformat) - if retval == 0: - pdflist.append((fname,thumb)) - else: - pdflist.append((fname,fname)) - if (rownum+1) % 2 == 0: - fhtml.append('<tr class="odd_row"><td><a href="%s">%s</a></td><td>%s</td></tr>' % (fname,fname,sfsize)) - else: - fhtml.append('<tr><td><a href="%s">%s</a></td><td>%s</td></tr>' % (fname,fname,sfsize)) - for logfname in logfiles: # expect at least tlog - if more - if os.path.abspath(logfname) == os.path.abspath(self.tlog): # handled later - sectionname = 'All tool run' - if (len(logfiles) > 1): - sectionname = 'Other' - ourpdfs = pdflist - else: - realname = os.path.basename(logfname) - sectionname = os.path.splitext(realname)[0].split('_')[0] # break in case _ added to log - ourpdfs = [x for x in pdflist if os.path.basename(x[0]).split('_')[0] == sectionname] - pdflist = [x for x in pdflist if os.path.basename(x[0]).split('_')[0] <> sectionname] # remove - nacross = 1 - npdf = len(ourpdfs) - - if npdf > 0: - nacross = math.sqrt(npdf) ## int(round(math.log(npdf,2))) - if int(nacross)**2 != npdf: - nacross += 1 - nacross = int(nacross) - width = min(400,int(1200/nacross)) - html.append('<div class="toolFormTitle">%s images and outputs</div>' % sectionname) - html.append('(Click on a thumbnail image to download the corresponding original PDF image)<br/>') - ntogo = nacross # counter for table row padding with empty cells - html.append('<div><table class="simple" cellpadding="2" cellspacing="2">\n<tr>') - for i,paths in enumerate(ourpdfs): - fname,thumb = paths - s= """<td><a href="%s"><img src="%s" title="Click to download a PDF of %s" hspace="5" width="%d" - alt="Image called %s"/></a></td>\n""" % (fname,thumb,fname,width,fname) - if ((i+1) % nacross == 0): - s += '</tr>\n' - ntogo = 0 - if i < (npdf - 1): # more to come - s += '<tr>' - ntogo = nacross - else: - ntogo -= 1 - html.append(s) - if html[-1].strip().endswith('</tr>'): - html.append('</table></div>\n') - else: - if ntogo > 0: # pad - html.append('<td> </td>'*ntogo) - html.append('</tr></table></div>\n') - logt = open(logfname,'r').readlines() - logtext = [x for x in logt if x.strip() > ''] - html.append('<div class="toolFormTitle">%s log output</div>' % sectionname) - if len(logtext) > 1: - html.append('\n<pre>\n') - html += logtext - html.append('\n</pre>\n') - else: - html.append('%s is empty<br/>' % logfname) - if len(fhtml) > 0: - fhtml.insert(0,'<div><table class="colored" cellpadding="3" cellspacing="3"><tr><th>Output File Name (click to view)</th><th>Size</th></tr>\n') - fhtml.append('</table></div><br/>') - html.append('<div class="toolFormTitle">All output files available for downloading</div>\n') - html += fhtml # add all non-pdf files to the end of the display - else: - html.append('<div class="warningmessagelarge">### Error - %s returned no files - please confirm that parameters are sane</div>' % self.opts.interpreter) - html.append(galhtmlpostfix) - htmlf = file(self.opts.output_html,'w') - htmlf.write('\n'.join(html)) - htmlf.write('\n') - htmlf.close() - self.html = html - - - def run(self): - """ - scripts must be small enough not to fill the pipe! - """ - my_env = os.environ.copy() - if self.treatbashSpecial and self.opts.interpreter in ['bash','sh']: - retval = self.runBash(pth) - else: - if self.opts.output_dir: - ste = open(self.elog,'w') - sto = open(self.tlog,'w') - sto.write('## Toolfactory running %s as %s script\n' % (self.toolname,self.opts.interpreter)) - sto.flush() - p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=ste,stdin=subprocess.PIPE,cwd=self.opts.output_dir,env=my_env) - else: - p = subprocess.Popen(self.cl,shell=False,stdin=subprocess.PIPE,env=my_env) - p.stdin.write(self.script) - p.stdin.close() - retval = p.wait() - if self.opts.output_dir: - sto.close() - ste.close() - # get stderr, allowing for case where it's very large - tmp_stderr = open( self.elog, 'rb' ) - stderr = '' - try: - while True: - stderr += tmp_stderr.read( buffsize ) - if not stderr or len( stderr ) % buffsize != 0: - break - except OverflowError: - pass - tmp_stderr.close() - if self.opts.make_HTML: - self.makeHtml() - return retval - - def runBash(self): - """ - cannot use - for bash so use self.sfile - """ - if self.opts.output_dir: - s = '## Toolfactory generated command line = %s\n' % ' '.join(self.cl) - ste = open(self.elog,'w') - sto = open(self.tlog,'w') - sto.write(s) - sto.flush() - p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=ste,cwd=self.opts.output_dir) - else: - p = subprocess.Popen(self.cl,shell=False) - retval = p.wait() - if self.opts.output_dir: - sto.close() - ste.close() - # get stderr, allowing for case where it's very large - tmp_stderr = open(self.elog, 'rb' ) - stderr = '' - try: - while True: - stderr += tmp_stderr.read( buffsize ) - if not stderr or len( stderr ) % buffsize != 0: - break - except OverflowError: - pass - tmp_stderr.close() - if self.opts.make_HTML: - self.makeHtml() - return retval - - -def main(): - u = """ - This is a Galaxy wrapper. It expects to be called by a special purpose tool.xml as: - <command interpreter="python">rgBaseScriptWrapper.py --script_path "$scriptPath" --tool_name "foo" --interpreter "Rscript" - </command> - """ - op = optparse.OptionParser() - a = op.add_option - a('--script_path',default=None) - a('--tool_name',default=None) - a('--interpreter',default=None) - a('--output_dir',default=None) - a('--output_html',default=None) - a('--input_tab',default="None") - a('--output_tab',default="None") - a('--user_email',default='Unknown') - a('--bad_user',default=None) - a('--make_Tool',default=None) - a('--make_HTML',default=None) - a('--help_text',default=None) - a('--tool_desc',default=None) - a('--new_tool',default=None) - a('--tool_version',default=None) - opts, args = op.parse_args() - assert not opts.bad_user,'UNAUTHORISED: %s is NOT authorized to use this tool until Galaxy admin adds %s to admin_users in universe_wsgi.ini' % (opts.bad_user,opts.bad_user) - assert opts.tool_name,'## Tool Factory expects a tool name - eg --tool_name=DESeq' - assert opts.interpreter,'## Tool Factory wrapper expects an interpreter - eg --interpreter=Rscript' - assert os.path.isfile(opts.script_path),'## Tool Factory wrapper expects a script path - eg --script_path=foo.R' - if opts.output_dir: - try: - os.makedirs(opts.output_dir) - except: - pass - r = ScriptRunner(opts) - if opts.make_Tool: - retcode = r.makeTooltar() - else: - retcode = r.run() - os.unlink(r.sfile) - if retcode: - sys.exit(retcode) # indicate failure to job runner - - -if __name__ == "__main__": - main() - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rg_nri.xml Fri Dec 19 05:57:40 2014 -0500 @@ -0,0 +1,627 @@ +<tool id="rg_nri" name="NRI" version="0.03"> + <description>and other model improvement measures</description> + <requirements> + <requirement type="package" version="3.1.1">R_3_1_1</requirement> + <requirement type="package" version="1.3.18">graphicsmagick</requirement> + <requirement type="package" version="9.10">ghostscript</requirement> + <requirement type="package" version="2.14">glmnet_lars_2_14</requirement> + </requirements> + <command interpreter="python"> + rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rg_NRI" + --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" + </command> + <inputs> + <param name="title" type="text" value="NRI test" size="80" label="Plot Title" help="Will appear as the title for the comparison plot"/> + <param name="input1" type="data" format="tabular" label="Select a tabular file from the baseline model with predicted and observed outcome column for each subject" + multiple='False' help="Observed and predicted status columns must be selected from this file below - NOTE both models must be in same order with exact matches in all observed outcomes" optional="False"/> + <param name="input1_observed" label="Select column containing observed outcome (0 for no event, 1 for an event)" type="data_column" data_ref="input1" numerical="True" + multiple="False" use_header_names="True" optional="False" help = "Observed outcomes are compared in the two files to check that the datasets are from the same data"/> + <param name="input1_predicted" label="Select column containing predicted event probabilies from baseline model" type="data_column" data_ref="input1" numerical="True" + multiple="False" use_header_names="True" optional="False" help="Must be in range 0-1"/> + <param name="input1_id" label="Select column containing subject ID from baseline model" type="data_column" data_ref="input1" numerical="True" + multiple="False" use_header_names="True" optional="False" help="Subect IDs are needed to match subjects to compare predictions in the two inputs"/> + <param name="input2" type="data" format="tabular" label="Select a tabular file from the new model with predicted and observed outcome columns for each subject" + multiple='False' help="Observed and predicted status columns must be selected from this file below" /> + <param name="input2_observed" label="Select column containing observed outcome (0 for no event, 1 for an event)" type="data_column" data_ref="input2" numerical="True" + multiple="False" use_header_names="True" optional="False" help = "Observed outcomes are compared in the two files to check that the datasets are from the same data"/> + <param name="input2_predicted" label="Select column containing predicted event probabilities from the new model" type="data_column" data_ref="input2" numerical="True" + multiple="False" use_header_names="True" optional="False" help="Must be in range 0-1"/> + <param name="input2_id" label="Select column containing subject ID from the new model" type="data_column" data_ref="input2" numerical="True" + multiple="False" use_header_names="True" optional="False" help="Subect IDs are needed to match subjects to compare predictions in the two inputs"/> + <conditional name="CImeth"> + <param name="cis" type="select" label="CI calculation method" + help="Bootstrap will take time - a long time for thousands - asymptotic is quick and informative"> + <option value="asymptotic" selected="true">Asymptotic estimate</option> + <option value="boot">Bootstrap for empirical CIs</option> + </param> + <when value="boot"> + <param name="nboot" type="integer" value="1000" label="Number of bootstrap replicates"/> + </when> + <when value="asymptotic"> + <param name="nboot" type="hidden" value="1000"/> + </when> + </conditional> + </inputs> + <outputs> + <data format="html" name="html_file" label="${title}.html"/> + <data format="tabular" name="nri_file" label="${title}_nrires.xls"/> + </outputs> + <tests> + <test> + <param name='input1' value='nri_test1.xls' ftype='tabular' /> + <param name='input2' value='nri_test2.xls' ftype='tabular' /> + <output name='html_file' file='nri_test1_out.html' compare='diff' lines_diff='10' /> + <output name='nri_file' file='nri_test1_out.xls' /> + </test> +</tests> +<help> + +**Before you start** + +This is a simple tool to calculate various measures of improvement in prediction between two models described in pickering_paper_ +It is based on an R script pickering_code_ written by Dr John W Pickering and Dr David Cairns from sunny Otago University which +has been debugged and slightly adjusted to fit a Galaxy tool wrapper. + + +**What it does** + +Copied from the documentation in pickering_code_ :: + + + Functions to create risk assessment plots and associated summary statistics + + + (c) 2012 Dr John W Pickering, john.pickering@otago.ac.nz, and Dr David Cairns + Last modified August 2014 + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + + FUNCTIONS + raplot + Produces a Risk Assessment Plot and outputs the coordinates of the four curves + Based on: Pickering, J. W. and Endre, Z. H. (2012). New Metrics for Assessing Diagnostic Potential of + Candidate Biomarkers. Clinical Journal of the American Society of Nephrology, 7, 1355–1364. doi:10.2215/CJN.09590911 + + statistics.raplot + Produces the NRIs, IDIs, IS, IP, AUCs. + Based on: Pencina, M. J., D'Agostino, R. B. and Steyerberg, E. W. (2011). Extensions of net reclassification improvement calculations to + measure usefulness of new biomarkers. Statistics in Medicine, 30(1), 11–21. doi:10.1002/sim.4085 + Pencina, M. J., D'Agostino, R. B. and Vasan, R. S. (2008). Evaluating the added predictive ability of a new marker: From area under the + ROC curve to reclassification and beyond. + Statistics in Medicine, 27(2), 157–172. doi:10.1002/sim.2929 + DeLong, E., DeLong, D. and Clarke-Pearson, D. (1988). Comparing the areas under 2 or more correlated receiver operating characteristic curves - a nonparametric approach. + Biometrics, 44(3), 837–845. + + summary.raplot + Produces the NRIs, IDIs, IS, IP, AUCs with confidence intervals using a bootstrap or asymptotic procedure. (I prefer bootstrap which is chosed by cis=c("boot")) + + + Required arguments for all functions: + x1 is calculated risk (eg from a glm) for the null model, i.e. predict(,type="response") on a glm object + x2 is calculated risk (eg from a glm) for the alternative model + y is the case-control indicator (0 for controls, 1 for cases) + Optional argument + t are the boundaries of the risks for each group (ie 0, 1 and the thresholds beteween. eg c(0,0,3,0,7,1)). If missing, defaults to c(0, the incidence, 1) + + +**Input** + +The observed and predicted outcomes from two models to be compared. + +**Output** + +Lots'o'measures (TM) see pickering_paper_ for details + +**Attributions** + +pickering_paper_ is the paper the caclulations performed by this tool is based on + +pickering_code_ is the R function from John Pickering exposed by this Galaxy tool with minor modifications and hacks by Ross Lazarus. + +Galaxy_ (that's what you are using right now!) for gluing everything together + +Otherwise, all code and documentation comprising this tool was written by Ross Lazarus and is +licensed to you under the LGPL_ like other rgenetics artefacts + +.. _LGPL: http://www.gnu.org/copyleft/lesser.html +.. _pickering_code: http://www.researchgate.net/publication/264672640_R_function_for_Risk_Assessment_Plot__reclassification_metrics_NRI_IDI_cfNRI +.. _pickering_paper: http://cjasn.asnjournals.org/content/early/2012/05/24/CJN.09590911.full +.. _Galaxy: http://getgalaxy.org + + +</help> + +<configfiles> +<configfile name="runme"> + +<![CDATA[ + +### http://www.researchgate.net/publication/264672640_R_function_for_Risk_Assessment_Plot__reclassification_metrics_NRI_IDI_cfNRI code +### http://cjasn.asnjournals.org/content/early/2012/05/24/CJN.09590911.full is the reference +### lots of little tweaks and but fixes. Using t as a variable name seems fraught to me. +### Ross Lazarus october 2014 for a Galaxy tool wrapper using the toolfactory infrastucture + +############################################################################# +###Functions to create risk assessment plots and associated summary statistics +############################################################################# +### +### (c) 2012 Dr John W Pickering, john.pickering@otago.ac.nz, and Dr David Cairns +### Last modified August 2014 +### +### Redistribution and use in source and binary forms, with or without +### modification, are permitted provided that the following conditions are met: +### * Redistributions of source code must retain the above copyright +### notice, this list of conditions and the following disclaimer. +### * Redistributions in binary form must reproduce the above copyright +### notice, this list of conditions and the following disclaimer in +### the documentation and/or other materials provided with the distribution + +### FUNCTIONS +### raplot +### Produces a Risk Assessment Plot and outputs the coordinates of the four curves +### Based on: Pickering, J. W. and Endre, Z. H. (2012). New Metrics for Assessing Diagnostic Potential of +### Candidate Biomarkers. Clinical Journal of the American Society of Nephrology, 7, 1355–1364. doi:10.2215/CJN.09590911 +### +### statistics.raplot +### Produces the NRIs, IDIs, IS, IP, AUCs. +### Based on: Pencina, M. J., D'Agostino, R. B. and Steyerberg, E. W. (2011). +### Extensions of net reclassification improvement calculations to measure usefulness of new biomarkers. Statistics in Medicine, 30(1), 11–21. doi:10.1002/sim.4085 +### Pencina, M. J., D'Agostino, R. B. and Vasan, R. S. (2008). Evaluating the added predictive ability of a new marker: From area under the ROC curve to reclassification and beyond. +### Statistics in Medicine, 27(2), 157–172. doi:10.1002/sim.2929 +### DeLong, E., DeLong, D. and Clarke-Pearson, D. (1988). Comparing the areas under 2 or more correlated receiver operating characteristic curves - a nonparametric approach. +### Biometrics, 44(3), 837–845. +### +### summary.raplot +### Produces the NRIs, IDIs, IS, IP, AUCs with confidence intervals using a bootstrap or asymptotic procedure. (I prefer bootstrap which is chosed by cis=c("boot")) + + +### Required arguments for all functions: +### x1 is calculated risk (eg from a glm) for the null model, i.e. predict(,type="response") on a glm object +### x2 is calculated risk (eg from a glm) for the alternative model +### y is the case-control indicator (0 for controls, 1 for cases) +### Optional argument +### t are the boundaries of the risks for each group (ie 0, 1 and the thresholds beteween. eg c(0,0,3,0,7,1)). If missing, defaults to c(0, the incidence, 1) + + +### risk assessment plot + +library('e1071') +library('caret') +library('pROC') +library('Hmisc') +library('pracma') + +raplot = function(x1, x2, y, outplot, title) { + + roc.model1 = roc(y, x1) + roc.model2 = roc(y, x2) + sens.model1 = roc.model1\$sensitivities + spec.model1 = 1 - roc.model1\$specificities + n.model1 = length(sens.model1) + thresh.model1 = roc.model1\$thresholds + thresh.model1 = thresh.model1[c(-1,-n.model1)] + sens.model1 = sens.model1[c(-1,-n.model1)] + spec.model1 = spec.model1[c(-1,-n.model1)] + sens.model2 = roc.model2\$sensitivities + spec.model2 = 1 - roc.model2\$specificities + n.model2 = length(sens.model2) + thresh.model2 = roc.model2\$thresholds + thresh.model2[1]=0 + thresh.model2[length(thresh.model2)]=1 + thresh.model2 = thresh.model2[c(-1,-n.model2)] + sens.model2 = sens.model2[c(-1,-n.model2)] + spec.model2 = spec.model2[c(-1,-n.model2)] + + n.model1 = length(sens.model1) + n.model2 = length(sens.model2) + + ### actual plotting + pdf(outplot) + plot(thresh.model1, sens.model1, xlim = c(0, 1), ylim = c(0, 1), type = "n", + lty = 2, lwd = 2, xlab = "Risk of Event", ylab = "", col = "black", main=title) + grid() + + polygon(x = c(thresh.model1, thresh.model2[n.model2:1]), + y = c(sens.model1, sens.model2[n.model2:1]), border = NA, col = gray(0.8)) + polygon(x = c(thresh.model1, thresh.model2[n.model2:1]), + y = c(spec.model1, spec.model2[n.model2:1]), border = NA, col = gray(0.8)) + + lines(thresh.model1, sens.model1, type = "l", lty = 2, lwd = 2, col = "black") + lines(thresh.model2, sens.model2, type = "l", lty = 1, lwd = 2, col = "black") + + lines(thresh.model1, spec.model1, type = "l", lty = 2, lwd = 2, col = "red") + lines(thresh.model2, spec.model2, type = "l", lty = 1, lwd = 2, col = "red") + + text(x = -0.15, y = 0.4, labels = "Sensitivity, ", col = "black", xpd = TRUE, srt = 90) + text(x = -0.15, y = 0.4 + 0.175, labels = "1-Specificity", col = "red", xpd = TRUE, srt = 90) + legend("topleft", c("Event: New model", "Event: Baseline model", + "No Event: New model", "No Event: Baseline model"), + col = c("black", "black", "red", "red"), + lty = c(1,2, 1, 2), lwd = 2, bty = "n") + dev.off() + return(data.frame("Null.p.sens"=thresh.model1, + "Null.sens"=sens.model1, + "Null.p.1spec"=thresh.model1, + "Null.1spec"=sens.model1, + "Alt.p.sens"=thresh.model2, + "Alt.sens"=sens.model2, + "Alt.p.1spec"=thresh.model2, + "Alt.1spec"=sens.model2)) + +} + + + +### statistics from a raplot (is an adaptation of improveProb() from Hmisc) + +statistics.raplot = function(x1, x2, y, threshvec) +{ + + s = is.na(x1 + x2 + y) ###Remove rows with missing data + if (any(s)) { + smiss = sum(s) + s = !s + x1 = x1[s] + x2 = x2[s] + y = y[s] + print.noquote(paste('Warning: removed',smiss,'cases with missing values')) + } + n = length(y) + y = as.numeric(y) + u = sort(unique(y)) + if (length(u) != 2 || u[1] != 0 || u[2] != 1) { + print.noquote("INPUT ERROR: y must have only two values: 0 and 1") + sink() + quit(save="no",status=2) + } + r = range(x1, x2) + if (r[1] < 0 || r[2] > 1) { + print.noquote("INPUT ERROR: x1 and x2 must be in [0,1]") + sink() + quit(save="no",status=3) + } + incidence=sum(y)/n + if (missing(threshvec)) { + threshvec=c(0, incidence,1) + print(paste('threshvec missing. using',paste(threshvec,collapse=','))) + } + a = (y == 1) + b = (y == 0) + na = sum(a) + nb = sum(b) + d = x2 - x1 + ### NRI + n.thresh=length(threshvec)-1 + risk.class.x1.ev=cut2(x1[a],threshvec) + risk.class.x2.ev=cut2(x2[a],threshvec) + thresh=c() + lt = length(threshvec) + for (i in 1:(lt-1)) { + thresh[i] = paste("[",toString(threshvec[i]),",",toString(threshvec[i+1]),"]") + } + levels(risk.class.x1.ev)=thresh + levels(risk.class.x2.ev)=thresh + cM.ev=confusionMatrix(risk.class.x2.ev,risk.class.x1.ev) + pup.ev=0 + pdown.ev=0 + for (i in 1:(n.thresh-1)) { pup.ev = pup.ev + sum(cM.ev\$table[(i+1):n.thresh,i])} + for (i in 2:n.thresh) { pdown.ev = pdown.ev + sum(cM.ev\$table[1:(i-1),i])} + pup.ev=pup.ev/na + pdown.ev=pdown.ev/na + risk.class.x1.ne=cut2(x1[b],threshvec) + risk.class.x2.ne=cut2(x2[b],threshvec) + levels(risk.class.x1.ne)=thresh + levels(risk.class.x2.ne)=thresh + cM.ne=confusionMatrix(risk.class.x2.ne,risk.class.x1.ne) + pup.ne=0 + pdown.ne=0 + for (i in 1:(n.thresh-1)){pup.ne=pup.ev+sum(cM.ne\$table[(i+1):n.thresh,i])} + for (i in 2:n.thresh){pdown.ne=pdown.ne+sum(cM.ne\$table[1:(i-1),i])} + pdown.ne=pdown.ne/nb + pup.ne=pup.ne/nb + nri = pup.ev - pdown.ev - (pup.ne - pdown.ne) + se.nri = sqrt((pup.ev + pdown.ev)/na + (pup.ne + pdown.ne)/nb) + z.nri = nri/se.nri + nri.ev = pup.ev - pdown.ev + se.nri.ev = sqrt((pup.ev + pdown.ev)/na) + z.nri.ev = nri.ev/se.nri.ev + nri.ne = pdown.ne - pup.ne + se.nri.ne = sqrt((pdown.ne + pup.ne)/nb) + z.nri.ne = nri.ne/se.nri.ne + ### Category Free NRI calculations + cfpup.ev = mean(d[a] > 0) + cfpup.ne = mean(d[b] > 0) + cfpdown.ev = mean(d[a] < 0) + cfpdown.ne = mean(d[b] < 0) + cfnri = cfpup.ev - cfpdown.ev - (cfpup.ne - cfpdown.ne) + se.cfnri = sqrt((cfpup.ev + cfpdown.ev)/na + (cfpup.ne + cfpdown.ne)/nb) + z.cfnri = cfnri/se.cfnri + cfnri.ev = cfpup.ev - cfpdown.ev + se.cfnri.ev = sqrt((cfpup.ev + cfpdown.ev)/na) + z.cfnri.ev = cfnri.ev/se.cfnri.ev + cfnri.ne = cfpdown.ne - cfpup.ne + se.cfnri.ne = sqrt((cfpdown.ne + cfpup.ne)/nb) + z.cfnri.ne = cfnri.ne/se.cfnri.ne + ### IDI calculations + improveSens = sum(d[a])/na + improveSpec = -sum(d[b])/nb + idi.ev = mean(improveSens) + idi.ne = mean(improveSpec) + idi = idi.ev - idi.ne + var.ev = var(d[a])/na + se.idi.ev = sqrt(var.ev) + z.idi.ev = idi.ev/se.idi.ev + var.ne = var(d[b])/nb + se.idi.ne = sqrt(var.ne) + z.idi.ne = idi.ne/se.idi.ne + se.idi = sqrt(var.ev + var.ne) + z.idi = idi/se.idi + ### AUC calculations + roc.x1 = roc(y, x1) + auc.x1 = auc(roc.x1) + ci.auc.x1 = ci.auc(roc.x1) + se.auc.x1 = (ci.auc.x1[3] - auc.x1)/qnorm(0.975) + roc.x2 = roc(y, x2) + auc.x2 = auc(roc.x2) + ci.auc.x2 = ci.auc(roc.x2) + se.auc.x2 = (ci.auc.x2[3] - auc.x2)/qnorm(0.975) + roc.test.x1.x2 = roc.test(roc.x1, roc.x2) ###Uses the default Delong method + sens.x1 = roc.x1\$sensitivities + spec.x1 = 1 - roc.x1\$specificities + n.x1 = length(sens.x1) + x1 = roc.x1\$thresholds + x1 = x1[c(-1,-n.x1)] + sens.x1 = sens.x1[c(-1,-n.x1)] + spec.x1 = spec.x1[c(-1,-n.x1)] + sens.x2 = roc.x2\$sensitivities + spec.x2 = 1 - roc.x2\$specificities + n.x2 = length(sens.x2) + x2 = roc.x2\$thresholds + x2 = x2[c(-1,-n.x2)] + sens.x2 = sens.x2[c(-1,-n.x2)] + spec.x2 = spec.x2[c(-1,-n.x2)] + ### Integrated sensitivity and 1-specificity calculations + is.x1 = trapz(x = x1, y = sens.x1) ### area under curves (relates to integrated sens, 1-spec) + is.x2 = trapz(x = x2, y = sens.x2) + ip.x1 = trapz(x = x1, y = spec.x1) + ip.x2 = trapz(x = x2, y = spec.x2) + + ### Output + output = c(n, na, nb, pup.ev, pup.ne, pdown.ev, pdown.ne, nri, se.nri, z.nri, + nri.ev, se.nri.ev, z.nri.ev, nri.ne, se.nri.ne, z.nri.ne, + cfpup.ev, cfpup.ne, cfpdown.ev, cfpdown.ne, cfnri, se.cfnri, z.cfnri, + cfnri.ev, se.cfnri.ev, z.cfnri.ev, cfnri.ne, se.cfnri.ne, z.cfnri.ne, + improveSens, improveSpec, idi.ev, se.idi.ev, z.idi.ev, idi.ne, + se.idi.ne, z.idi.ne, idi, se.idi, z.idi, is.x1, NA, is.x2, NA, + ip.x1, NA, ip.x2, NA, auc.x1, se.auc.x1, auc.x2, se.auc.x2, + roc.test.x1.x2\$p.value,incidence) + names(output) = c("n", "na", "nb", "pup.ev", "pup.ne", "pdown.ev", "pdown.ne", + "nri", "se.nri", "z.nri", "nri.ev", "se.nri.ev", "z.nri.ev", + "nri.ne", "se.nri.ne", "z.nri.ne", + "cfpup.ev", "cfpup.ne", "cfpdown.ev", "cfpdown.ne", + "cfnri", "se.cfnri", "z.cfnri", "cfnri.ev", "se.cfnri.ev", "z.cfnri.ev", + "cfnri.ne", "se.cfnri.ne", "z.cfnri.ne", "improveSens", "improveSpec", + "idi.ev", "se.idi.ev", "z.idi.ev", "idi.ne", "se.idi.ne", + "z.idi.ne", "idi", "se.idi", "z.idi", "is.x1", "se.is.x1", + "is.x2", "se.is.x2", "ip.x1", "se.ip.x1", "ip.x2", "se.ip.x2", + "auc.x1", "se.auc.x1", "auc.x2", "se.auc.x2", + "roc.test.x1.x2.pvalue","incidence") + resdf = data.frame(N=n, Na=na, Nb=nb, pup.ev=pup.ev, pup.ne=pup.ne, pdown.ev=pdown.ev, pdown.ne=pdown.ne, NRI=nri, NRI.se=se.nri, NRI.z=z.nri, + NRI.ev=nri.ev, NRI.ev.se=se.nri.ev, NRI.ev.z=z.nri.ev, NRI.ne=nri.ne, NRI.ne.se=se.nri.ne, NRI.ne.z=z.nri.ne, + cfpup.ev=cfpup.ev, cfpup.ne=cfpup.ne, cfpdown.ev=cfpdown.ev, cfpdown.ne=cfpdown.ne, CFNRI=cfnri, CFNRI.se=se.cfnri, CFNRI.z=z.cfnri, + CFNRI.ev=cfnri.ev, CFNRI.ev.se=se.cfnri.ev, CFNRI.ev.z=z.cfnri.ev, CFNRI.ne=cfnri.ne, CFNRI.ne.se=se.cfnri.ne, CFNRI.ne.z=z.cfnri.ne, + improvSens=improveSens, improvSpec=improveSpec, IDI.ev=idi.ev, IDI.ev.se=se.idi.ev, IDI.ev.z=z.idi.ev, IDI.ne=idi.ne, + IDI.ne.se=se.idi.ne, IDI.ne.z=z.idi.ne, IDI=idi, IDI.se=se.idi, IDI.z=z.idi, isx1=is.x1, isx2=is.x2, + ipxi=ip.x1, ipx2=ip.x2, AUC.x1=auc.x1, AUC.x1.se=se.auc.x1, AUC.x2=auc.x2, AUC.x2.se=se.auc.x2, + roctestpval=roc.test.x1.x2\$p.value,incidence=incidence) + tr = t(resdf) + tresdf = data.frame(measure=colnames(resdf),value=tr[,1]) + return(list(resdf=tresdf,output=output)) +} + + +### More comprehensive summary statistics from a raplot +### Choice of confidence intervals determined through asymptotics or bootstrapping (n.boot = ### of bootstrap resamples) +### dp is number of decimal places for results table + +summary.raplot = function(x1, x2, y, threshvec, cis = "boot", conf.level = 0.95, n.boot = 2000, dp = 4, stat_ra=NA) +{ + results = stat_ra + if (cis == "boot") { + print.noquote("Bootstrap estimates for SE") + results.boot = matrix(NA, n.boot, length(names(results))) + + colnames(results.boot) = names(results) + + for (i in 1:n.boot) { + ###boot.index = sample(length(cc.status), replace = TRUE) + ###risk.model1.boot = risk.model1[boot.index] + ###risk.model2.boot = risk.model2[boot.index] + ###cc.status.boot = cc.status[boot.index] + boot.index = sample(length(y), replace = TRUE) + risk.model1.boot = x1[boot.index] + risk.model2.boot = x2[boot.index] + cc.status.boot = y[boot.index] + r = statistics.raplot(x1 = risk.model1.boot, x2 = risk.model2.boot, y = cc.status.boot) + results.boot[i, ] = r\$output + } + + results.se.boot = apply(results.boot, 2, sd) + print(paste(results.se.boot,collapse=',')) + + + results[grep("se", names(results))] = results.se.boot[grep("se", names(results)) - 1] + + } + + + + ### calculate cis and return + + z = abs(qnorm((1 - conf.level)/2)) + + results.matrix = matrix(NA, 24, 2) + + results.matrix[1, ] = c("Total (n)", results["n"]) + results.matrix[2, ] = c("Events (n)", results["na"]) + results.matrix[3, ] = c("Non-events (n)", results["nb"]) + results.matrix[4, ] = c("Category free NRI and summary statistics","-------------------------") + results.matrix[5, ] = c("cfNRI events (%)", + paste(round(100*results["cfnri.ev"], dp-2), " (", + round(100*results["cfnri.ev"] - z * 100*results["se.cfnri.ev"], dp-2), + " to ", round(100*results["cfnri.ev"] + + z * 100*results["se.cfnri.ev"], dp-2), ")", sep = "")) + results.matrix[6, ] = c("cfNRI non-events (%)", + paste(round(100*results["cfnri.ne"], dp-2), " (", + round(100*results["cfnri.ne"] - z * 100*results["se.cfnri.ne"], dp)-2, + " to ", round(100*results["cfnri.ne"] + z * 100*results["se.cfnri.ne"], + dp-2), ")", sep = "")) + results.matrix[7, ] = c("cfNRI (%)", + paste(round(100*results["cfnri"], dp-2), " (", + round(100*results["cfnri"] - z * 100*results["se.cfnri"], dp-2), + " to ", round(100*results["cfnri"] + z * 100*results["se.cfnri"], + dp-2), ")", sep = "")) + results.matrix[8, ] = c("NRI and summary statistics","-------------------------") + results.matrix[9, ] = c("NRI events (%)", + paste(round(100*results["nri.ev"], dp-2), " (", + round(100*results["nri.ev"] - z * 100*results["se.nri.ev"], dp-2), + " to ", round(100*results["nri.ev"] + + z * 100*results["se.nri.ev"], dp-2), ")", sep = "")) + results.matrix[10, ] = c("NRI non-events (%)", + paste(round(100*results["nri.ne"], dp-2), " (", + round(100*results["nri.ne"] - z * 100*results["se.nri.ne"], dp-2), + " to ", round(100*results["nri.ne"] + z * 100*results["se.nri.ne"], + dp-2), ")", sep = "")) + results.matrix[11, ] = c("NRI (%)", + paste(round(100*results["nri"], dp-2), " (", + round(100*results["nri"] - z * 100*results["se.nri"], dp-2), + " to ", round(100*results["nri"] + z * 100*results["se.nri"], + dp-2), ")", sep = "")) + results.matrix[12, ] = c("IDI and summary statistics","-------------------------") + results.matrix[13, ] = c("IDI events", + paste(round(results["idi.ev"], dp), " (", + round(results["idi.ev"] - z * results["se.idi.ev"], dp), + " to ", round(results["idi.ev"] + z * results["se.idi.ev"], + dp), ")", sep = "")) + results.matrix[14, ] = c("IDI non-events", + paste(round(results["idi.ne"], dp), " (", + round(results["idi.ne"] - z * results["se.idi.ne"], dp), + " to ", round(results["idi.ne"] + z * results["se.idi.ne"], + dp), ")", sep = "")) + results.matrix[15, ] = c("IDI", + paste(round(results["idi"], dp), " (", + round(results["idi"] - z * results["se.idi"], dp), + " to ", round(results["idi"] + z * results["se.idi"], + dp), ")", sep = "")) + results.matrix[16, ] = c("IS (null model)", + paste(round(results["is.x1"], dp), " (", + round(results["is.x1"] - z * results["se.is.x1"], dp), + " to ", round(results["is.x1"] + z * results["se.is.x1"], + dp), ")", sep = "")) + results.matrix[17, ] = c("IS (alt model)", + paste(round(results["is.x2"], dp), " (", + round(results["is.x2"] - z * results["se.is.x2"], dp), + " to ", round(results["is.x2"] + z * results["se.is.x2"], + dp), ")", sep = "")) + results.matrix[18, ] = c("IP (null model)", + paste(round(results["ip.x1"], dp), " (", + round(results["ip.x1"] - z * results["se.ip.x1"], dp), + " to ", round(results["ip.x1"] + z * results["se.ip.x1"], + dp), ")", sep = "")) + results.matrix[19, ] = c("IP (alt model)", + paste(round(results["ip.x2"], dp), " (", + round(results["ip.x2"] - z * results["se.ip.x2"], dp), + " to ", round(results["ip.x2"] + z * results["se.ip.x2"], + dp), ")", sep = "")) + results.matrix[20, ] = c("AUC","-------------------------") + results.matrix[21, ] = c("AUC (null model)", + paste(round(results["auc.x1"], dp), " (", + round(results["auc.x1"] - z * results["se.auc.x1"], dp), + " to ", round(results["auc.x1"] + z * results["se.auc.x1"], + dp), ")", sep = "")) + results.matrix[22, ] = c("AUC (alt model)", + paste(round(results["auc.x2"], dp), " (", + round(results["auc.x2"] - z * results["se.auc.x2"], dp), + " to ", round(results["auc.x2"] + z * results["se.auc.x2"], + dp), ")", sep = "")) + results.matrix[23, ] = c("difference (P)", round(results["roc.test.x1.x2.pvalue"], dp)) + results.matrix[24, ] = c("Incidence", round(results["incidence"], dp)) + + return(results.matrix) +} + + + +]]> + +options(width=120) +options(digits=5) +logf = file("rgNRI.log", open = "a") +sink(logf,type = c("output", "message")) +Out_Dir = "$html_file.files_path" +Input1 = "$input1" +Input2 = "$input2" +myTitle = "$title" +outtab = "$nri_file" +input1_obs = $input1_observed +input1_pred = $input1_predicted +input1_id = $input1_id +input2_obs = $input2_observed +input2_pred = $input2_predicted +input2_id = $input2_id +in1 = read.table(Input1,head=T,sep='\t') +in2 = read.table(Input2,head=T,sep='\t') +id1 = in1[,input1_id] +id2 = in2[,input2_id] +useme1 = in1[which(id1 %in% id2),] +useme2 = in2[which(id2 %in% id1),] +id1 = useme1[,input1_id] +id2 = useme2[,input2_id] +useme1 = useme1[order(id1),] +useme2 = useme2[order(id2),] +x1 = useme1[,input1_pred] +x2 = useme2[,input2_pred] +y1 = useme1[,input1_obs] +y2 = useme2[,input2_obs] +n.boot = $CImeth.nboot +conf.level = 0.95 +cis = "$CImeth.cis" +digits = 4 +nydiff = sum(y1 != y2) +if (nydiff > 0) { + print.noquote(paste('Input error: observed status column has',nydiff,'differences - cannot reliably proceed')) + quit(save="no",status=1) + } +y = y2 +outplot = 'rgNRI_EventRisk.pdf' +res = raplot(x1=x1, x2=x2, y=y, outplot=outplot,title=myTitle) + +stats = statistics.raplot(x1=x1, x2=x2, y=y) +res1 = stats\$resdf +out1 = stats\$output +print.noquote('Results:') +print.noquote(res1,digits=4) +res2 = summary.raplot(x1=x1, x2=x2, y=y, cis = cis, conf.level = conf.level, n.boot = n.boot, dp = digits, stat_ra=out1) +print.noquote('Summary:') +print.noquote(res2,digits=4) +write.table(format(res1,digits=4),outtab,quote=F, col.names=F,sep="\t",row.names=F) +print.noquote('SessionInfo for this R session:') +sessionInfo() +print.noquote('warnings for this R session:') +warnings() +sink() +</configfile> +</configfiles> +<citations> + <citation type="doi">doi: 10.2215/​CJN.09590911</citation> +</citations> +</tool> + + +
--- a/rglasso_cox.xml Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,547 +0,0 @@ -<tool id="rglasso_cox" name="Lasso" version="0.03"> - <description>and cox regression using elastic net</description> - <requirements> - <requirement type="package" version="3.1.1">R_3_1_1</requirement> - <requirement type="package" version="1.3.18">graphicsmagick</requirement> - <requirement type="package" version="9.10">ghostscript</requirement> - <requirement type="package" version="2.14">glmnet_lars_2_14</requirement> - </requirements> - <command interpreter="python"> - rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rglasso" - --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" - </command> - <inputs> - <param name="title" type="text" value="lasso test" size="80" label="Title for job outputs" help="Typing a short, meaningful text here will help remind you (and explain to others) what the outputs represent"> - <sanitizer invalid_char=""> - <valid initial="string.letters,string.digits"><add value="_" /> </valid> - </sanitizer> - </param> - <param name="input1" type="data" format="tabular" label="Select an input tabular text file from your history. Rows represent samples; Columns are measured phenotypes" - multiple='False' help="Tabular text data with samples as rows, phenotypes as columns with a header row of column identifiers" /> - <param name="xvar_cols" label="Select columns containing numeric variables to use as predictor (x) variables" type="data_column" data_ref="input1" numerical="True" - multiple="True" use_header_names="True" /> - <param name="force_xvar_cols" label="Select numeric columns containing variables ALWAYS included as predictors in cross validation" type="data_column" data_ref="input1" numerical="True" - multiple="True" use_header_names="True" optional="True"/> - <conditional name="model"> - <param name="fam" type="select" label="GLM Link function for models" - help="Binary dependant variables will automatically be set to Binomial no matter what this is set to"> - "gaussian","binomial","poisson","multinomial","cox","mgaussian" - <option value="gaussian" selected="true">Gaussian - continuous dependent (y)</option> - <option value="binomial">Binomial dependent variables</option> - <option value="poisson">Poisson (eg counts)</option> - <option value="cox">Cox models - require special setup for y variables - see below</option> - </param> - <when value="gaussian"> - <param name="yvar_cols" label="Select numeric columns containing variables to use as the dependent (y) in elasticnet" type="data_column" data_ref="input1" numerical="True" - multiple="True" use_header_names="True" help = "If multiple, each will be modelled against all the x variables and reported separately."/> - </when> - <when value="binomial"> - <param name="yvar_cols" label="Select numeric columns containing variables to use as the dependent (y) in elasticnet" type="data_column" data_ref="input1" numerical="True" - multiple="True" use_header_names="True" help = "If multiple, each will be modelled against all the x variables and reported separately."/> - </when> - <when value="poisson"> - <param name="yvar_cols" label="Select columns containing variables to use as the dependent (y) in elasticnet" type="data_column" data_ref="input1" numerical="True" - multiple="True" use_header_names="True" help = "If multiple, each will be modelled against all the x variables and reported separately."/> - </when> - <when value="cox"> - <param name="cox_time" label="Select column containing time under observation for Cox regression" - type="data_column" data_ref="input1" numerical="True" optional="False" - multiple="False" use_header_names="True" help = "This MUST contain a time period - eg continuous years or days to failure or right censoring"/> - <param name="cox_status" label="Select column containing status = 1 for outcome of interest at the end of the time under observation or 0 for right censoring" - type="data_column" data_ref="input1" numerical="True" optional="False" - multiple="False" use_header_names="True" help = "This MUST contain 1 for subjects who had an event at that time or 0 for a right censored observation"/> - <param name="output_full" type="select" label="Create a tabular output with coefficients for all predictors" > - <option value="F" selected="true">No full model output file</option> - <option value="T">Create a full model output file</option> - </param> - </when> - </conditional> - <param name="logtrans" type="select" label="Perform a log transformation on all predictors" help="A tiny number (1e-10) will be added to prevent taking logs of zero"> - <option value="False" selected="true">No log transformation of predictors</option> - <option value="True">Log 2 transform predictors before modelling - will fail if negative values!</option> - </param> - <param name="do_standard" type="select" label="Standardise x vars" - help="If all measurements on same scale, may not be needed. Coefficients are always returned on the original scale."> - <option value="False" selected="true">No standardisation of predictors</option>l - - <option value="True">Standardise predictors before model</option> - </param> - <param name="alpha" type="float" value="0.95" size="5" min="0.01" max="1.0" label="Alpha - see glmnet docs. 1 for pure lasso. 0.0 for pure ridge regression" - help="Default 0.95 allows lasso to cope better with expected predictor collinearity. Use (eg) 0.5 for hybrid regularised regression or (eg) 0.025 for ridge regression"/> - <param name="nfold" type="integer" value="10" size="5" label="Number of folds for internal cross validation" - help="Default of 10 is usually ok"/> - </inputs> - <outputs> - <data format="html" name="html_file" label="${title}.html"/> - <data format="tabular" name="model_file" label="${title}_modelres.xls"/> - <data format="tabular" name="full_cox" label="${title}_full_cox_model.xls"> - <filter>model['output_full'] == 'T'</filter> - </data> - </outputs> - <tests> - <test> - <param name='input1' value='cox_test.xls' ftype='tabular' /> - <param name='treatment_name' value='case' /> - <param name='title' value='Cox glmnet test' /> - <param name='nfold' value='10' /> - <param name='logtrans' value='False' /> - <param name='alpha' value='0.95' /> - <param name='do_standard' value="True" /> - <param name='cox_time' value='1' /> - <param name='cox_status' value='2' /> - <param name='fam' value='cox' /> - <param name='model.yvar_cols' value='' /> - <param name='xvar_cols' value='3,4,5' /> - <param name='force_xvar_cols' value='3' /> - <param name='output_full' value='0' /> - <output name='model_file'> - <assert_contents> - <has_text text="rhubarb" /> - <has_text text="TRUE" /> - <!-- 	 is XML escape code for tab --> - <has_line line="regulator	partial_likelihood	forced_in	glmnet_model	best_lambda" /> - <has_n_columns n="5" /> - </assert_contents> - </output> - <output name='html_file' file='coxlassotest.html' compare='diff' lines_diff='10' /> - </test> -</tests> -<help> - -**Before you start** - -Please read the glmnet documentation @ glmnet_ - -This Galaxy wrapper merely exposes that code and the glmnet_ documentation is essential reading -before getting useful results here. - -**What it does** - -From documentation at glmnet_ :: - - Glmnet is a package that fits a generalized linear model via penalized maximum likelihood. - The regularization path is computed for the lasso or elasticnet penalty at a grid of values for the regularization parameter lambda. - The algorithm is extremely fast, and can exploit sparsity in the input matrix x. - It fits linear, logistic and multinomial, poisson, and Cox regression models. - A variety of predictions can be made from the fitted models. - -All solutions starting with all zero regression parameters, allowing one to become non-zero at a time -in a forward stepwise manner, are estimated by coordinate descent, each point corresponding to a value for lambda. - -Internal cross validation is used to optimise the choice of lambda based on CV AUC for logistic (binomial outcome) models, or CV mse for gaussian. - -Sufficiently big jobs will take a while (eg each lasso regression with 20k features on 1k samples takes about 2-3 minutes on our aged cluster) - -**Input** - -Assuming you have more measurements than samples, you supply data as a tabular text file where each row is a sample and columns -are variables. You specify which columns are dependent (predictors) and which are observations for each sample. Each of multiple -dependent variable columns will be run and reported independently. Prectors can be forced in to the model. - -**Output** - -For each selected dependent regression variable, a brief report of the model coefficients predicted at the -'optimal' nfold CV value of lambda. - -**Attributions** - -glmnet_ is the R package exposed by this Galaxy tool. - -Galaxy_ (that's what you are using right now!) for gluing everything together - -Otherwise, all code and documentation comprising this tool was written by Ross Lazarus and is -licensed to you under the LGPL_ like other rgenetics artefacts - -.. _LGPL: http://www.gnu.org/copyleft/lesser.html -.. _glmnet: http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html -.. _Galaxy: http://getgalaxy.org - - -</help> - -<configfiles> -<configfile name="runme"> -<![CDATA[ -library('glmnet') -library('lars') - - -mdsPlot = function(dm,myTitle,groups,outpdfname,transpose=T) -{ - - samples = colnames(dm) - mt = paste(unlist(strsplit(myTitle,' ')),collapse=" ") - pcols=c('maroon') - if (! is.na(groups)) - { - gu = unique(groups) - colours = rainbow(length(gu),start=0.1,end=0.9) - pcols = colours[match(groups,gu)] - } - mydata = dm - if (transpose==T) - { - mydata = t(dm) - } - d = dist(mydata) - fit = cmdscale(d,eig=TRUE, k=10) - - x = fit\$points[,1] - y = fit\$points[,2] - pdf(outpdfname) - plot(x, y, xlab="Dimension 1", ylab="Dimension 2", - main=paste(mt,"MDS Plot"),type="n", col=pcols, cex=0.3) - text(x, y, labels = row.names(mydata), cex=0.3, col=pcols) - grid(col="lightgray",lty="dotted") - dev.off() -} - - -dolasso_cox = function(x,y,debugOn=F,maxsteps=10000,nfold=10,xcolnames,ycolnames,out_full=F,out_full_file=NA, - descr='Cox test',logtransform=F,do_standard=F,alpha=0.9,penalty=c()) -{ - logf = file("cox_rglasso.log", open = "wt") - sink(logf,type = c("output", "message")) - res = NULL - - ok = complete.cases(x,y) - if (sum(! ok) > 0) { - y = y[(ok),] - x = x[(ok),] - print(paste("Removed",sum(! ok),'cases with any missing variables')) - } - if (sum(ok) == 0 ) { - print(paste("No complete cases found for cox models in input x dim =",paste(dim(xm),collapse=','),"length y=",length(y))) - sink - return(NA) - } - do_standard = do_standard - standardize = do_standard - normalize = do_standard - larsres = try(glmnet(x,y,family='cox',standardize=standardize,alpha=alpha ),T) - if (class(larsres) == "try-error") - { - print.noquote('Unable to run cox glmnet on your data') - print(larsres) - sink() - return(NA) - } - if (out_full == T) - { - b = as.matrix(larsres\$beta) - nb = length(colnames(b)) - bcoef = b[,nb] - lastl = larsres\$lambda[length(larsres\$lambda)] - allres = data.frame(x=rownames(b),beta=bcoef,lambda=lastl) - write.table(allres,out_full_file,quote=FALSE, sep="\t",row.names=F) - } - - outpdf = paste('cox',descr,'glmnetdev.pdf',sep='_') - try( - { - pdf(outpdf) - plot(larsres,main='cox glmnet',label=T) - grid() - dev.off() - },T) - - larscv = NA - - larscv = try(cv.glmnet(x,y,family=fam,type.measure='deviance',penalty=penalty),T) - if (class(larscv) == "try-error") { return(NA) } - outpdf = paste('cox',descr,'glmnet_cvdeviance.pdf',sep='_') - - try( - { - pdf(outpdf) - plot(larscv,main='Deviance',label=T) - grid() - dev.off() - },T) - best_lambda = larscv\$lambda.min - bestcoef = coef(larscv, s = "lambda.min") - inmodel = which(bestcoef != 0) - coefs = bestcoef[inmodel,1] - preds = rownames(bestcoef)[inmodel] - names(coefs) = preds - pen = as.logical( ! penalty[inmodel]) - if (debugOn) { - print.noquote(paste('best_lambda=',best_lambda,'saving cox respreds=',paste(names(coefs),collapse=','),'as predictors of survival. Coefs=',paste(coefs,collapse=','))) - } - res = try(data.frame(regulator=names(coefs),partial_likelihood=coefs,forced_in=pen,glmnet_model='cox',best_lambda=best_lambda),T) - if (class(res) == "try-error") { return(NA) } - if (debugOn) { - print.noquote('CV results:') - print.noquote(res) - } - sink() - return(res) - -} - -do_lasso = function(x=NA,y=NA,do_standard=T,debugOn=T,defaultFam="gaussian",logtransform=T,descr='description', - indx=1,target='target',sane=F,alpha=1.0,nfold=10,penalty=c()) -{ - logf = file(paste(target,"rglasso.log",sep='_'), open = "wt") - sink(logf,type = c("output", "message")) - res = NA - phe_is_bin = (length(unique(y)) == 2) - forcedin = paste(colnames(x)[which(penalty == 0)],collapse=',') - print.noquote(paste('i=',indx,'target=',target,'logtrans=',logtransform,'is binary=',phe_is_bin,'dim(x)=',paste(dim(x),collapse=','),'length(y)=',length(y),'force=',forcedin)) - fam = "gaussian" - if (defaultFam %in% c("poisson","binomial","gaussian","multinomial")) fam=defaultFam - if (phe_is_bin == T) { - fam = "binomial" - } - standardize = do_standard - if (fam == "binomial") - { - larsres = try(glmnet(x,y,family=fam,standardize=standardize,maxit=10000,alpha=alpha,type.logistic = "modified.Newton" ),T) - } else { - larsres = try(glmnet(x,y,family=fam,standardize=standardize,maxit=10000,alpha=alpha,type.gaussian="covariance"),T) - } - if (class(larsres) == "try-error") - { - print.noquote(paste('Unable to run glmnet on your data for',target)) - sink() - return(NA) - } - mt = paste(descr,'glmnet on',target) - outpdf = paste(target,descr,'glmnetPath.pdf',sep='_') - pdf(outpdf) - try( - { - plot(larsres,main=mt,label=T) - grid() - dev.off() - },T) - outpdf = paste(target,descr,'glmnetDeviance.pdf',sep='_') - pdf(outpdf) - mt2 = paste(descr,'Deviance for',target) - try( { - plot(larsres,xvar="dev",main=mt2,label=T) - grid() - dev.off() - },T) - larscv = NA - if (fam=="binomial") { - tmain = paste(target,'AUC') - outpdf = paste(target,descr,'glmnetCV_AUC.pdf',sep='_') - larscv = try(cv.glmnet(x=x,y=y,family=fam,type.measure='auc',penalty=penalty),T) - } else { - tmain = paste(target,'CV MSE') - outpdf = paste(target,descr,'glmnetCV_MSE.pdf',sep='_') - larscv = try(cv.glmnet(x,y,family=fam,type.measure='mse',penalty=penalty),T) - } - if (class(larscv) == "try-error") - { - print(paste('ERROR: unable to run cross validation for target',target,'error=',larscv)) - sink() - return(NA) - } - pdf(outpdf) - try( - { - plot(larscv,main=tmain) - grid() - dev.off() - },T) - ipenalty = c(0,penalty) ### must include intercept which is always forcedq - best_lambda = larscv\$lambda.min - bestpred = as.matrix(coef(larscv, s = "lambda.min")) - inmodel = which(bestpred != 0) - coefs = bestpred[inmodel,1] - preds = rownames(bestpred)[inmodel] - iforced = ipenalty[inmodel] - forced = ! as.logical(iforced) - names(coefs) = preds - if (debugOn) {cat(indx,'best_lambda=',best_lambda,'saving',fam,'respreds=',names(coefs),'as predictors of',target,'coefs=',coefs,'\n')} - res = try(data.frame(i=indx,pred=target,regulator=names(coefs),coef=coefs,forced_in=forced,glmnet_model=fam),T) - if (class(res) == "try-error") { - sink() - return(NA) } - print.noquote(res) - sink() - return(res) -} - -dolasso_generic = function(predvars=NA,depvars=NA,debugOn=T,maxsteps=100, alpha=alpha,nfold=10,xcolnames=c(),ycolnames=c(), - descr="describe me",logtransform=F,do_standard=F,defaultFam="gaussian",penalty=penalty) -{ - logf = file("rglasso.log", open = "a") - sink(logf,type = c("output", "message")) - xdat = predvars - if (logtransform) { - small = 1e-10 - try( { lpred = log(predvars) },T) - if (class(lpred) == "try-error") - try( { lpred = log(predvars+small) }, T) - if (class(lpred) == "try-error") - { - print.noquote('Unable to logtransform your data') - return(NA) - } - descr = paste('logx',descr,sep='_') - } - - xm = as.matrix(xdat) - res = NA - depnames = ycolnames - ndep = length(depnames) - dm = xm - outpdfname = 'rglasso_x_in_sample_space_MDS.pdf' - mdsPlot(dm,'measurements in sample space',groups=NA,outpdfname=outpdfname,transpose=T) - outpdfname = 'rglasso_samples_in_x_space_MDS.pdf' - mdsPlot(dm,'samples in measurement space',groups=NA,outpdfname=outpdfname,transpose=F) - for (i in c(1:max(1,ndep))) { - target = depnames[i] - if (length(target) < 1) { target='y' } - if (i %% 100 == 0) { cat(i,target,'\n') } - if (ndep <= 1) { - y=depvars - } else { - y = depvars[,i] - } - x = xm - ok = complete.cases(x,y) - if (sum(! ok) > 0) { - y = y[(ok)] - x = xm[(ok),] - } - if (sum(ok) == 0 ) { - print(paste("No complete cases found for",target,"in input x dim =",paste(dim(xm),collapse=','),"length y=",length(y))) - } else { - regres = do_lasso(x=x,y=y,do_standard=do_standard,debugOn=debugOn,defaultFam=defaultFam, - logtransform=logtransform,descr=descr,indx=i,target=target,alpha=alpha,nfold=nfold,penalty=penalty) - if (! is.na(regres)) { res = rbind(res,regres) } - } - - } - if (! is.na(res)) { - rownames(res) = c(1:length(res\$i)) - } - sink() - return(res) -} - - - - -corPlot=function(xdat=c(),main='main title',is_raw=T) -{ - library(pheatmap) - library(gplots) - if (is_raw) { - cxdat = cor(xdat,method="spearman",use="pairwise.complete.obs") - } else { - cxdat=xdat - } - xro = nrow(cxdat) - if (xro > 1000) stop("Too many rows for heatmap, who can read?!") - fontsize_col = 5.0 - pheatmap(cxdat, main=main, show_colnames = F, width=30, height=30, - fontsize_row=fontsize_col, border_color=NA) -} - - - - -runTest = function(n=10) -{ - set.seed (NULL) - Y = data.frame(y1=runif (n),y2=runif(n)) - Xv <- runif(n*n) - X <- matrix(Xv, nrow = n, ncol = n) - - mydf <- data.frame(Y, X) - - regres_out = dolasso_generic(predvars=X,depvars=Y,debugOn=T,p.cutoff = 0.05,maxsteps=10000,nfold=10, - descr='randomdata',logtransform=F,do_standard=do_standard,defaultFam="gaussian",alpha=0.05) - return(regres_out) -} -]]> -options(width=120) -alpha = $alpha -nfold = $nfold -Out_Dir = "$html_file.files_path" -Input = "$input1" -myTitle = "$title" -outtab = "$model_file" -logtrans = as.logical("$logtrans") -do_standard = as.logical("$do_standard") -fam = "$model.fam" -xvar_cols_in = "$xvar_cols" -force_xvar_cols_in = "$force_xvar_cols" -xvar_cols = as.numeric(strsplit(xvar_cols_in,",")[[1]]) -force_xvar_cols = c() -if (force_xvar_cols_in > "") -{ - force_xvar_cols = as.numeric(strsplit(force_xvar_cols_in,",")[[1]]) - allx = c(xvar_cols,force_xvar_cols) - xvar_cols = unique(allx) - xvar_cols = xvar_cols[order(xvar_cols)] -} - -indat = read.table(Input,head=T,sep='\t') -datcols = colnames(indat) - -x = indat[,xvar_cols] -xcolnames = datcols[xvar_cols] -penalties = rep(1,length(datcols)) -penalty = penalties[xvar_cols] -if (force_xvar_cols_in > "") -{ - penalties[force_xvar_cols] = 0 - penalty = penalties[xvar_cols] - forcedin = paste(datcols[which(penalties == 0)],collapse=',') -} - -if (file.exists(Out_Dir) == F) dir.create(Out_Dir) -out_full = F -out_full_file = NA -#if $model.fam == "cox": - cox_time = $model.cox_time - cox_status = $model.cox_status - out_full = as.logical("$model.output_full") - if (out_full == T) { out_full_file="$full_cox" } - yvar_cols = c(cox_time,cox_status) - ycolnames = c('time','status') - istat = as.double(indat[,cox_status]) - itime = as.double(indat[,cox_time]) - y = data.frame(time = itime, status = istat) - ustat = unique(istat) - if ((length(ustat) != 2) | (! 1 %in% ustat ) | (! 0 %in% ustat)) - { - print.noquote(paste('INPUT ERROR: status must have 0 (censored) or 1 (event) but found',paste(ustat,collapse=',') )) - quit(save='no',status=1) - } - y = as.matrix(y) - x = as.matrix(x) - forced = paste(force_xvar_cols,collapse=',') - print.noquote(paste('Cox model: Got yvar=',paste(ycolnames,collapse=','),'cols=',paste(yvar_cols,collapse=','),'n preds=',length(xcolnames), - 'forced in=',forced)) - regres_out = dolasso_cox(x=x,y=y,debugOn=F,maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames,out_full=out_full,out_full_file=out_full_file, - descr=myTitle,logtransform=logtrans,do_standard=do_standard,alpha=alpha,penalty=penalty) -#else: - yvar_cols = "$model.yvar_cols" - yvar_cols = as.numeric(strsplit(yvar_cols,",")[[1]]) - y = indat[,yvar_cols] - ycolnames = colnames(indat)[yvar_cols] - print.noquote(paste('Got yvar=',paste(ycolnames,collapse=','),'cols',paste(yvar_cols,collapse=','),'n preds=',length(xcolnames),'forced in=',paste(force_xvar_cols,collapse=','))) - regres_out = dolasso_generic(predvars=x,depvars=y,debugOn=F, maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames, - descr=myTitle,logtransform=logtrans,do_standard=do_standard,defaultFam=fam,alpha=alpha,penalty=penalty) -#end if -print.noquote('Results preview:') -print.noquote(regres_out) -write.table(regres_out,outtab,quote=FALSE, sep="\t",row.names=F) -print.noquote('SessionInfo for this R session:') -sessionInfo() -warnings() - -</configfile> -</configfiles> - -</tool> - - -
--- a/test-data/cox_test.xls Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,502 +0,0 @@ -time status rhubarb vegemite apple -575.966708 0 452.405468 30.339584 32286.089057 -1539.245319 1 329.689929 30.603839 15735.863202 -2072.798422 0 534.379263 32.474983 22639.136685 -1638.450154 0 362.522161 30.925996 31108.364370 -1771.625630 0 417.167751 31.444652 31317.491931 -4413.484706 1 459.947526 30.069033 11232.037882 -1519.376431 1 576.868414 30.102709 16539.948172 -4805.620290 0 554.756164 30.685550 9979.125726 -657.497124 0 439.481951 32.305768 22301.057122 -181.498205 1 581.498091 32.549519 10499.737992 -1630.552909 0 580.670671 30.825057 22321.821037 -2285.043061 1 588.034529 32.799742 8067.710675 -3697.405336 0 400.195630 31.826099 20975.730587 -3821.226480 0 561.466202 32.333230 17887.547850 -4810.837196 0 447.004938 32.547717 25383.925946 -3195.877681 1 482.169788 31.420740 27659.409774 -4556.750833 1 415.144072 31.960725 18401.138992 -3797.192061 1 484.492359 30.789882 7431.214452 -3734.484446 1 587.339066 32.052766 29443.923798 -567.665985 1 377.599777 31.203617 30148.076752 -2551.367122 0 447.623689 31.795708 15922.972458 -3.724234 1 534.080187 31.175131 8259.786167 -3528.906130 0 471.293653 31.225006 18307.543198 -3649.459162 1 443.804109 32.390564 9720.974163 -1490.640766 1 456.837379 32.099378 15292.782705 -4023.629633 1 443.869307 30.022985 23145.513801 -1389.288557 0 497.850711 30.361120 28757.227200 -535.427950 1 437.634304 32.307966 18289.733136 -3115.819953 0 587.593182 31.317353 21533.733397 -1561.895566 1 472.250365 31.550546 8724.067134 -2579.785003 0 331.955830 32.992176 12105.988070 -1592.819387 0 361.242026 31.011367 12823.339723 -4473.715814 0 334.539193 30.448049 9621.600487 -4328.198133 0 537.782739 32.683122 19090.292579 -3761.244687 1 392.916387 32.665889 22119.375302 -4638.730675 0 514.168146 30.548414 20651.266975 -1531.289211 0 451.404672 32.437627 21748.461852 -3844.140593 0 427.404840 32.894070 25824.198908 -4542.766339 0 366.992470 30.814372 6090.765585 -3296.236419 0 575.910564 32.509193 4050.864668 -710.915095 0 395.062092 32.940952 9086.464980 -3211.769396 0 428.608771 30.918481 29996.819084 -4448.442217 0 300.153203 31.160941 8735.458385 -4215.329607 0 504.064526 32.040819 29036.000397 -3249.725224 0 489.538009 32.863547 27520.919952 -3388.033079 0 485.938597 30.920035 20395.348731 -1465.239279 0 435.641977 30.193093 20469.748911 -2022.519649 0 362.011214 32.756787 18367.094394 -1474.662644 0 448.727985 31.194883 18456.157990 -2680.178585 0 553.521680 32.339108 17386.094437 -254.953131 1 467.930054 32.705913 28823.817908 -2080.749522 0 387.588403 30.520279 3389.418182 -4609.398123 1 471.568096 31.101077 27574.777543 -2506.254687 1 466.233702 32.928613 24693.847389 -1418.686322 0 519.402535 32.685931 9055.300824 -1716.501599 0 305.671024 31.453871 25382.505028 -2080.668376 0 589.514396 31.143248 24546.670311 -4051.441999 1 431.028785 30.854820 32291.736106 -221.739141 1 435.123612 30.465214 19077.175842 -3945.564918 0 559.566964 32.871624 18320.616591 -3543.190633 1 548.960006 30.269711 6259.113819 -4141.263169 1 418.013926 30.159280 12806.181141 -4613.643735 1 467.593153 30.584036 28998.377681 -2526.696757 1 442.602416 31.075723 24392.589804 -661.056610 0 306.186743 30.071900 4931.492293 -420.059232 0 557.043309 30.889234 30410.871000 -3777.621646 1 564.397202 30.810439 25359.485010 -4210.010862 0 541.943504 31.109403 18240.759045 -2251.752128 1 442.108369 31.308836 14136.141189 -4037.391365 1 467.892619 32.089289 14540.468274 -2478.012769 0 305.600832 30.653583 7117.580065 -2435.693015 0 551.918277 30.355254 29969.743529 -3816.690227 0 314.516061 31.511033 9880.067119 -531.400047 0 565.141705 30.722177 11252.478408 -4367.861436 0 421.500739 31.696812 16879.486148 -634.857316 1 476.587731 30.741441 31351.695315 -2757.953233 1 484.171370 31.510266 7767.961089 -676.813083 0 412.948756 30.174301 11123.425270 -3672.369695 0 374.771725 30.531392 32349.940422 -2572.842389 0 554.137199 32.777709 12727.962812 -4960.699040 0 560.164589 30.206776 3453.492367 -4316.475328 1 484.168551 32.069106 31984.273406 -4015.881070 1 400.274241 30.658886 11783.352436 -3052.451190 1 398.216589 30.358316 27928.803920 -1398.823206 0 446.521178 30.185715 31193.500034 -3922.909840 1 575.609185 32.724842 10233.363838 -182.803149 0 370.875903 30.856525 23954.116536 -2632.176915 1 426.948756 30.641145 30860.455084 -1060.390433 0 300.347626 31.143117 23615.898639 -3178.728283 1 514.261332 31.356224 26944.147061 -763.218442 1 463.367562 31.274306 7509.665226 -953.708474 0 384.301310 31.065789 15102.331493 -4434.527995 1 431.967331 32.369140 7177.785265 -623.914857 0 495.381873 32.915097 22739.985454 -3157.386533 1 515.293513 32.210421 31073.191185 -521.293142 0 593.449301 30.668076 29606.157493 -2845.587324 1 367.180408 32.424788 7052.350189 -4832.119270 1 465.679571 30.019029 5249.534096 -2856.445302 1 591.090965 31.133381 32498.377584 -2572.527257 1 409.021421 31.885768 13743.016360 -1413.194580 1 330.844512 32.740409 5441.606816 -2918.271731 0 304.169439 31.256918 22384.990257 -1674.024061 1 589.100300 32.763039 13305.351548 -2808.138084 0 491.357574 31.150814 17202.472869 -1424.663181 0 359.990029 31.343898 28564.024266 -2718.446177 1 580.778432 31.842978 29231.311831 -3735.730841 0 573.549637 30.079704 31952.172593 -4641.585855 0 365.172295 32.894223 17289.384320 -4071.637795 0 505.251583 32.472570 9633.518042 -471.788564 1 322.257982 32.546866 30076.426907 -1266.346159 1 335.749836 31.646670 32477.736291 -2695.370600 1 509.475837 31.947731 22351.940960 -3631.475811 0 584.575140 30.366829 13700.704821 -488.174184 0 522.383785 30.637536 29035.941989 -2645.805895 1 390.848623 31.639290 3004.701971 -4037.703358 1 425.444485 32.775358 28720.752294 -2196.610999 1 547.619331 32.640598 29108.134173 -1081.264261 1 350.044271 32.021783 19844.288986 -1074.305164 1 497.881588 32.048148 16325.089356 -1643.851420 0 372.156064 30.815571 12924.845649 -4919.287039 0 470.345027 31.043515 6953.198088 -3961.891671 1 353.750322 32.386195 10074.185995 -1643.427847 0 521.758887 31.784944 10303.308946 -2714.177818 0 363.519932 31.954710 4089.970927 -1958.737870 0 372.752456 32.778447 25735.088968 -72.430537 1 374.015311 30.406160 3519.910494 -666.824775 0 443.814245 31.382885 32222.429084 -179.203304 0 339.335011 30.911475 29294.393947 -3193.179459 1 329.239965 32.270432 7539.157039 -1405.955780 1 575.062041 31.769778 26056.208652 -4378.622335 1 583.779711 30.280724 12958.062814 -4933.985803 0 474.998284 30.835755 8298.116085 -479.970673 0 311.469610 31.110600 16498.771798 -194.917302 1 533.735116 32.682195 17363.293404 -3126.230367 1 568.043474 32.209623 22525.398240 -3605.775657 0 519.068951 30.750083 21588.128190 -68.236883 1 316.956971 32.735304 5303.699685 -2956.632516 1 359.835070 31.735527 26701.676218 -1273.924283 1 346.804469 31.185254 29977.854210 -1530.707021 1 568.442488 31.510433 20716.251168 -3847.696838 1 487.304670 32.242322 5276.023316 -3121.510941 0 518.612631 30.448810 10702.449428 -2693.401144 1 533.008711 32.964097 31390.744850 -1774.742268 1 588.145493 30.487710 21800.531169 -2216.272824 0 467.093275 31.533770 7991.454914 -1526.314833 0 433.703579 32.809884 11209.353005 -2089.304275 0 321.084472 31.021678 27856.425967 -3046.990363 1 460.307711 32.631130 7182.074497 -1317.038946 0 593.108863 31.275244 26560.416891 -1352.101840 0 306.529955 30.979817 21571.506540 -3163.866962 0 304.262526 30.303302 24717.381911 -3417.487105 1 367.150160 31.345893 13987.483107 -4582.950578 0 478.091198 32.280854 23336.546777 -1183.214096 1 435.761798 31.182880 32631.545653 -549.532217 0 436.107384 32.534391 23409.562965 -4018.079504 1 514.150908 32.582211 11007.862101 -4537.862661 0 474.720396 31.380999 6160.620917 -2023.882772 1 576.659734 31.873967 24564.400593 -301.339657 0 525.192309 30.612934 19844.142010 -989.521636 0 323.610331 30.471812 6456.387005 -3647.414661 1 353.160773 30.687803 20893.860619 -184.390450 0 367.556798 32.894773 7004.911101 -2851.709060 1 420.820133 32.238140 4715.183972 -3988.318894 0 352.860472 31.705620 18426.123706 -2244.796525 0 332.927674 31.617697 20768.384056 -726.779082 0 370.641547 32.668547 5629.086601 -2379.683675 1 566.970531 30.287987 26943.017249 -2039.968848 0 550.042019 31.203815 6686.226792 -949.545757 1 302.786734 30.845201 7393.079915 -2509.519553 1 378.264518 31.980310 16762.542739 -2949.998835 0 472.701800 32.896473 12708.457939 -3958.575946 1 402.750276 30.318068 10667.350272 -550.167142 0 394.589505 30.431805 32262.087079 -1462.278890 1 534.727448 31.514861 27759.667272 -2696.454522 1 535.656182 32.096399 32431.581354 -3756.029506 1 332.505291 32.486460 26887.712404 -1948.023625 1 470.099353 31.100519 31437.137850 -2008.758594 1 489.846790 31.656842 18513.611970 -889.052994 0 378.045215 30.360086 30999.351259 -956.860624 0 498.874188 31.238091 32986.760129 -1596.180150 1 468.477647 32.876509 10708.862352 -2044.623310 0 341.364451 30.831209 4377.873160 -1670.382802 0 344.655922 31.528413 25861.700033 -2751.419144 0 525.374048 31.427506 23893.998314 -1139.794772 0 510.822662 32.192052 27853.949410 -1951.275598 1 355.760660 30.362940 26312.111688 -3652.432308 0 330.985338 32.585428 25167.952966 -500.427055 1 453.806271 30.774837 16661.345691 -2884.409660 0 529.815414 31.683479 3317.970560 -2365.322179 0 540.382433 31.999929 29498.364669 -3853.672492 1 367.171163 32.671321 7512.370365 -938.719173 0 569.038125 31.387703 23444.121607 -3058.327259 1 411.675057 31.700811 28738.835254 -4215.442937 0 425.942827 31.617841 19505.555179 -1639.143130 1 310.802062 31.996376 24378.789768 -3015.816015 1 397.865213 30.774914 32208.055822 -1361.137094 1 314.093792 30.698064 21667.130306 -1900.651819 1 582.449408 30.643572 26755.743826 -2038.834435 1 359.535511 31.758170 26389.582514 -4217.395801 1 532.389233 32.367523 28128.901119 -4048.824843 1 415.408543 30.661216 6695.688224 -3278.940683 0 590.282168 31.148497 17706.150830 -2590.370816 1 467.475790 32.894899 31984.209283 -2961.953200 1 409.350826 30.412768 26800.587858 -757.305077 0 531.036597 32.692111 23056.009036 -4823.679979 0 547.031467 31.700254 4679.396986 -2159.887871 1 461.722894 31.554845 23118.215328 -1990.158592 0 436.583574 32.276433 22516.752446 -4010.528172 0 315.169873 31.546055 4521.472077 -527.527936 0 429.812828 32.161079 26933.467576 -2517.395024 1 442.643281 32.903345 11566.976782 -3973.530406 0 444.298671 30.860264 30738.467879 -909.115735 1 548.278268 31.969666 3061.067051 -1003.138129 0 599.492409 32.731826 5268.191640 -496.214281 1 525.539783 32.288968 18314.877509 -3321.884314 1 518.880235 31.593901 27585.403619 -531.060305 0 311.564795 30.660864 26630.390327 -1985.865186 0 468.954849 30.022861 7126.223629 -1075.080978 0 566.983564 32.548614 4189.576013 -2899.914447 0 501.981481 31.522821 3699.352956 -1209.686792 1 397.096875 30.589148 26885.709293 -1061.484419 0 399.802412 31.376845 21631.836258 -2341.148192 0 599.041273 32.705976 7273.381892 -2350.754927 1 578.118121 32.505017 16776.353449 -3715.290307 1 374.328239 31.642496 24018.734911 -353.266473 1 306.976330 30.935877 17228.695180 -739.123546 1 529.022431 30.659845 18290.776722 -465.537407 0 308.998940 31.824504 3898.463033 -3524.840397 1 385.275512 30.348328 6065.321632 -751.616003 0 322.281940 30.576017 20828.524508 -4989.405682 1 381.074571 30.437314 16844.697669 -488.563077 0 318.125000 31.138994 28854.306278 -4558.130890 0 479.058472 31.076337 3808.209751 -3644.754920 0 437.984557 30.038643 24302.584238 -4595.447449 0 323.287824 31.904699 24911.538572 -1064.168831 1 541.625023 30.168168 12380.733046 -4197.068049 0 418.764938 31.133507 14481.902057 -3593.303853 0 545.685419 30.904012 27860.611155 -3842.381434 0 596.665096 31.650039 32071.065698 -781.763222 0 544.786508 32.895055 4905.803177 -3941.386557 1 362.501808 32.489031 8386.976863 -4285.917626 1 373.105515 31.909823 18733.487831 -4152.435642 1 440.036083 30.439913 4669.690889 -2326.656610 0 409.061810 31.219523 17268.378337 -887.243489 1 458.991888 32.904742 32560.558550 -2383.154926 0 402.990252 31.501777 32165.195285 -1575.020117 1 567.576879 31.223480 11851.574021 -2681.793474 0 382.275765 30.672297 31505.220294 -3745.755933 1 524.831384 31.755858 20277.860327 -2310.372497 0 381.887619 32.549019 23884.208743 -1934.629896 0 430.992059 30.090286 16064.349955 -228.815046 1 510.344577 30.272399 13314.791074 -106.494754 0 421.895617 30.249703 6022.257615 -107.758807 1 352.048161 31.015999 27010.353099 -2366.497161 1 589.257754 32.599191 29060.870687 -3652.314844 1 325.392214 30.378902 9562.454775 -4628.382437 0 389.680904 30.488690 21333.408550 -3992.522719 1 460.945057 30.259246 21474.935828 -3116.879118 0 470.551054 30.384387 23511.934422 -3669.891947 1 468.454570 30.792238 6696.682439 -2455.615625 1 536.995869 30.130133 27385.537567 -334.254171 0 510.509005 31.195189 15166.910603 -2392.466005 1 503.464732 31.609503 28757.028622 -563.741633 0 535.281334 31.242527 11455.193241 -2683.204132 0 530.400221 30.122120 21995.611523 -116.891969 0 458.723324 30.381523 5852.247596 -2817.127744 1 547.927894 32.877328 22278.789600 -323.704077 0 432.630710 31.643437 8793.259513 -2223.194666 1 584.593612 31.798636 4232.455807 -4836.806252 0 328.278741 31.551587 28407.563886 -2302.391657 1 304.186017 30.683099 21760.214902 -1597.778913 1 593.734794 32.275770 26311.401980 -3231.228313 0 469.955959 31.155902 21789.637085 -4329.078171 1 326.299894 30.629070 28598.928501 -3004.531335 0 345.517768 31.874809 20075.187439 -3426.658667 0 407.558711 32.882949 32961.004750 -2942.293151 1 474.436479 30.681717 10297.918976 -2790.908728 1 372.634696 31.883438 18617.426989 -3645.736003 0 338.768900 30.271779 13340.959790 -4954.230936 0 475.204218 32.305344 31368.726270 -1374.169666 0 322.794825 32.054615 29507.251607 -952.311143 0 466.066676 31.283949 31527.207980 -1582.935333 1 430.497522 32.636945 6365.677867 -1574.488078 1 512.547475 32.000036 3214.800133 -3675.314941 0 420.357568 30.365188 24034.010725 -3737.900636 0 310.463467 32.944674 19212.098849 -453.453905 1 307.998144 31.734006 6696.312537 -180.430188 0 500.121289 31.647589 25439.918801 -2223.386978 0 342.168201 31.413275 10321.858985 -3347.152707 0 450.255234 32.419409 21979.139997 -2032.988552 1 359.268392 32.444386 4151.429949 -1855.874006 1 486.902129 32.194183 28734.066535 -576.472090 0 491.131111 30.250436 21183.624322 -2629.739140 1 490.880496 30.706030 20741.966223 -3234.449151 0 591.702238 31.565719 9467.501732 -4568.779861 0 391.821980 31.778473 12171.711242 -1369.135728 1 352.424832 32.501343 24914.054783 -1203.818730 1 453.408627 30.021814 9344.375802 -1328.193901 1 326.714975 31.957412 22506.010391 -3143.392191 0 488.722719 30.906428 20404.329854 -1516.300551 1 524.551137 30.783381 26784.582992 -4401.773609 1 490.129039 31.474209 27673.989380 -471.342147 0 345.653763 30.714091 27614.068357 -2846.556659 0 387.773914 32.101213 25758.358083 -3113.332960 0 340.017571 31.250862 8517.797718 -569.863929 1 507.285752 32.723247 7350.729977 -1351.433366 1 355.255156 31.629009 10035.728199 -2057.479962 0 504.044203 31.725231 14699.764489 -2029.207221 1 326.501058 31.883646 25956.067323 -4987.734533 0 391.140321 30.918262 15723.785267 -2728.957211 1 498.936107 32.017269 21531.521853 -2885.926679 0 458.513681 31.221876 23111.352960 -1214.796671 0 531.593343 32.409879 26260.051485 -4073.228494 1 563.045005 31.629048 5742.795393 -632.393708 1 428.087852 32.360714 20236.813182 -4548.671698 0 387.645357 30.048297 28095.140888 -735.589017 0 510.876342 30.202972 28887.002260 -2683.928468 1 352.751328 31.385745 3968.105767 -3761.692766 1 470.650685 31.713693 21645.943688 -3949.678238 0 572.913454 32.230970 18919.434641 -1695.346973 0 387.076598 32.045632 8265.498458 -4805.534106 0 322.040787 30.256378 20922.671974 -3834.130476 0 548.536590 30.581470 12530.815196 -3700.095266 0 434.154537 30.131024 28914.087338 -2083.132166 1 435.468758 31.961937 24155.150686 -2566.436697 1 455.368285 30.113525 23903.526223 -3372.416477 1 439.390127 32.380970 8312.312507 -3341.243804 1 496.573292 32.055384 9484.908468 -1289.638171 1 346.709851 31.116968 13853.905474 -441.798476 1 395.115019 32.820316 27178.509065 -3876.357092 1 501.568592 30.100469 3872.161975 -3082.204963 0 445.140351 31.403349 21141.060076 -1222.633709 1 510.019589 32.832816 12522.872061 -3257.092925 1 567.760441 30.935411 17780.698060 -4978.318869 0 488.434181 31.349966 20399.317688 -268.617967 0 563.410236 30.677835 13353.114481 -3668.370880 1 432.604476 31.336602 18298.862821 -1616.920156 1 528.362539 30.664065 18290.365790 -1582.470140 1 566.928554 31.829775 31526.762426 -1894.952274 0 576.611659 32.930078 5451.111476 -130.599064 1 357.672426 30.192175 5127.046275 -1658.869611 0 518.418202 30.222626 21787.150970 -1930.327624 0 539.744396 31.409771 11998.451350 -1360.478477 1 357.741711 32.358872 4894.722814 -2762.405397 1 498.202228 32.110196 20926.903902 -66.089100 0 525.721147 30.742496 28773.910408 -4617.260796 0 388.807122 31.052491 19010.672351 -4046.674710 1 363.050835 32.600458 8208.205191 -103.101768 0 401.874765 30.885852 29276.307064 -2870.170146 1 479.961929 31.997551 19841.555525 -1069.648823 1 598.655109 31.936390 10655.006702 -1988.650561 0 483.709224 31.937744 11014.379454 -2983.772162 0 328.380916 32.023817 8974.076899 -4534.002381 0 311.260565 30.136494 25947.551507 -3064.855784 1 369.080218 31.164070 29129.542680 -1746.878916 0 313.512424 30.992758 20731.078669 -4018.774695 1 497.866125 31.834090 3971.821700 -3742.871097 1 353.541863 30.115952 20815.885581 -2828.315773 1 454.655490 32.424599 26891.927439 -637.033003 1 424.553960 31.428681 20585.183759 -3087.856465 1 536.195443 31.597533 12661.823700 -604.873967 1 378.515493 30.925004 22439.100233 -3748.531033 1 515.115463 32.625710 6408.313067 -1558.535013 1 433.937532 32.380082 8383.657584 -2365.992186 0 409.216710 31.115409 32511.614098 -1183.022513 0 362.544664 30.959448 5065.906749 -3789.375438 0 595.724663 30.914010 20244.524890 -850.273991 1 387.136972 31.618840 5094.567581 -1231.544018 1 322.683917 31.457545 20767.599384 -1581.181089 1 305.296620 31.856051 17539.248499 -4880.147078 1 316.598007 31.317619 11086.097072 -4013.027211 0 485.648876 31.953258 30138.577997 -31.538821 0 486.966298 30.424936 3940.149283 -1735.735564 1 453.551097 30.972695 29273.785203 -4682.372724 1 426.064369 30.915967 29848.714453 -1419.673303 1 410.849400 31.957769 20423.937410 -4735.693777 0 331.994083 30.314568 21510.813276 -796.562868 0 502.136561 31.584336 14765.778268 -2637.389588 1 550.253910 31.759164 29206.835372 -1596.336472 0 367.679269 32.352327 5586.166324 -4673.788569 0 380.575280 32.897426 12723.911545 -2339.522676 0 405.444884 30.501805 20468.692308 -2935.072522 0 458.797029 30.585909 9955.177397 -2865.796979 0 497.765242 31.556959 27770.325805 -1032.502528 1 390.204545 31.050812 26330.543896 -1571.491134 1 454.897576 32.799730 4600.369582 -453.244362 0 513.834562 31.693168 18448.030313 -2566.490430 1 370.231825 31.371622 3530.831525 -2470.800433 1 506.755104 32.707049 15365.573438 -1674.578908 1 484.850818 31.657126 31034.118967 -1569.707807 0 380.339060 31.625216 26765.694375 -2048.318765 1 502.814703 32.061583 26794.094870 -4878.355942 1 308.718223 31.785469 9019.516491 -3394.948104 0 354.869800 31.553766 18461.635011 -3784.598905 0 506.145921 32.828899 22372.102070 -3648.947834 0 503.931222 32.468882 27102.757907 -1722.690965 1 389.086536 31.876640 5309.945334 -1266.660624 0 452.533007 30.468222 28572.903114 -794.701939 1 459.408945 32.075307 21311.040821 -1461.841041 0 370.527899 32.775555 10211.935095 -1976.852418 0 492.898069 30.817372 17854.725954 -859.962977 1 304.267257 30.348074 3439.856318 -2900.907710 1 396.919281 31.272958 4858.314485 -3693.615130 1 597.052502 30.903533 25716.531228 -1281.197465 1 470.854828 30.164368 9315.105349 -126.181622 1 381.237663 30.660067 28411.425041 -2517.244075 0 559.007071 30.373940 11135.538898 -2069.720704 1 544.708646 31.610539 16441.322232 -578.677623 1 591.867328 32.146388 21419.935158 -1918.809273 1 308.991680 31.457030 31410.129144 -1074.094412 0 336.860228 30.308287 17215.627286 -3034.994121 0 560.991229 31.862778 17867.418001 -4794.119298 1 480.292221 30.685158 11676.367480 -3307.911419 0 389.291016 30.558138 31108.361766 -457.891802 1 588.926068 32.330020 21119.054157 -1258.403059 0 332.016717 31.309245 16741.712001 -3312.429446 0 568.695418 31.010073 14155.107714 -3839.897879 0 370.810835 31.722406 16960.037012 -2068.214455 0 392.107828 30.058746 32853.457910 -2526.811696 0 510.287902 30.913056 22333.505644 -914.090787 1 323.087744 32.697177 14094.308393 -1890.199376 0 597.011219 31.259415 19326.502842 -1234.225123 0 536.105198 30.519191 9869.430841 -350.596239 1 320.162343 30.373749 5065.307636 -2846.841176 0 403.168597 31.624641 25354.766351 -1853.530932 1 564.965357 30.095729 11551.597590 -1183.876207 1 572.411649 32.013383 17167.971369 -3953.600477 0 404.253211 30.364392 17515.400202 -3289.696058 1 576.771217 31.015263 14654.467190 -2975.113909 1 587.011738 31.223106 7437.225919 -1404.141852 0 372.793625 30.325109 13500.525950 -1634.182114 0 569.705673 32.338410 27233.776345 -462.646125 1 467.892571 30.013450 27556.195432 -4138.443866 0 318.632678 30.348871 21795.067964 -4139.367443 1 304.173328 32.473154 29126.297266 -969.741373 1 599.886523 32.409973 16299.878199 -1494.463287 1 429.995117 30.484126 17461.996880 -4336.770459 1 416.241963 30.152207 32032.419174 -3064.782012 1 428.313584 31.688383 31731.608033 -2581.013362 0 490.480997 32.990983 17803.475008 -413.230864 1 430.932119 30.358681 18338.089371 -4469.602591 0 544.949864 30.599950 17579.676286 -2299.774314 1 532.263729 32.087371 25491.390739 -3570.957635 0 429.164688 30.798996 21783.084247 -2730.725224 0 434.512962 32.353024 8584.173078 -1859.613968 0 465.494240 30.853468 24199.875583 -3239.212454 1 585.752947 31.392592 17429.754381 -1626.243010 0 548.642457 32.586064 28320.002213 -4552.584886 0 411.908494 32.541147 18556.711124 -719.687465 1 437.155107 32.439905 18592.734838 -4528.629292 1 598.029535 30.261328 28073.717443 -1445.996134 1 407.243137 30.545202 28549.910302 -2143.373187 1 581.134422 30.728743 26291.822745 -4543.013656 0 406.670256 30.716994 24689.785142 -3479.794580 0 585.305633 31.528087 18146.824396 -2206.562888 0 526.292818 32.250928 6586.281804 -1796.678722 0 441.391846 32.716749 16294.189483 -2705.347154 0 551.894682 31.578422 11390.827042 -172.103511 1 320.919515 31.863888 29876.119334 -315.943273 1 512.197902 32.302445 17856.311999 -2750.695088 1 565.301177 31.514433 19939.221077 -4398.626158 0 418.977196 32.337330 4024.240756 -1522.050792 0 554.491924 32.249187 30588.435646 -3385.739992 1 496.030916 31.480018 17815.393181 -955.825566 1 393.353810 31.434487 26026.721071 -3971.619403 1 552.518065 30.319670 5068.002389 -2010.848206 1 569.817785 31.925723 26117.233780 -2181.752401 1 535.925058 30.952605 3509.328188 -905.194864 1 588.512406 30.229715 15643.423750 -3274.240917 0 590.044688 32.642379 13352.913630 -2061.196669 0 457.883033 30.745459 19655.487967 -4056.062516 0 354.337157 30.418176 24440.653038 -611.999600 0 432.598681 30.201091 19387.208118 -45.061679 1 437.997979 31.939322 21882.265038 -2091.842258 1 445.122275 32.963802 4496.459612 -2785.821696 1 535.953261 32.295237 5166.117813 -1272.363761 0 498.289820 31.101707 13536.602174 -2161.870929 0 487.379529 30.927239 22791.304718 -1230.121506 0 326.397921 31.953971 21880.979367 -2065.316190 0 405.291555 30.660933 28025.816633 -3230.173916 1 365.768097 30.055539 23327.051800 -4572.451029 1 323.399890 31.386897 21097.082477 -4916.172714 0 384.428209 32.754127 10981.900793 -4882.524497 0 399.434058 30.795923 18117.477126 -1457.322605 1 412.012101 30.834341 20230.063544 -1988.228627 0 566.906278 31.371606 15413.227181 -4337.262787 1 316.056609 32.556132 11299.235579 -2762.610907 1 454.177332 32.779747 14436.150973 -4756.268480 1 562.493863 32.727774 26209.569931 -3031.941096 0 313.505196 31.386484 9197.983181 -4586.200757 0 329.983728 30.115395 27039.250998 -300.892434 0 520.668472 32.713615 14579.471742 -392.962829 0 513.293270 30.568538 10311.182381 -2116.838847 1 355.318371 30.883049 18887.246656 -779.846218 1 407.952315 30.145534 22063.619725 -4236.655978 0 326.344600 32.797447 16004.751246 -4008.051049 1 319.382840 30.262423 20092.936132 -1550.891721 1 570.189663 32.238898 11686.528747 -2616.214276 1 448.726264 31.633295 26368.520969 -2775.847491 1 532.436984 30.938190 21751.684582 -
--- a/test-data/coxlassotest.html Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,140 +0,0 @@ -<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> - <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"> - <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> - <meta name="generator" content="Galaxy rgToolFactory.py tool output - see http://g2.trac.bx.psu.edu/" /> - <title></title> - <link rel="stylesheet" href="/static/style/base.css" type="text/css" /> - </head> - <body> - <div class="toolFormBody"> - -<div class="infomessage">Galaxy Tool "rglasso" run at 19/10/2014 13:11:14</div><br/> -<div class="toolFormTitle">cox images and outputs</div> -(Click on a thumbnail image to download the corresponding original PDF image)<br/> -<div><table class="simple" cellpadding="2" cellspacing="2"> -<tr> -<td><a href="cox_Coxglmnettest_glmnet_cvdeviance.pdf"><img src="cox_Coxglmnettest_glmnet_cvdeviance.png" title="Click to download a PDF of cox_Coxglmnettest_glmnet_cvdeviance.pdf" hspace="5" width="400" - alt="Image called cox_Coxglmnettest_glmnet_cvdeviance.pdf"/></a></td> - -<td><a href="cox_Coxglmnettest_glmnetdev.pdf"><img src="cox_Coxglmnettest_glmnetdev.png" title="Click to download a PDF of cox_Coxglmnettest_glmnetdev.pdf" hspace="5" width="400" - alt="Image called cox_Coxglmnettest_glmnetdev.pdf"/></a></td> -</tr> - -</table></div> - -<div class="toolFormTitle">cox log output</div> - -<pre> - -[1] "All partial likelihood coefficients from Cox model at minimum lambda of 0.00138091491730048:" - -3 x 1 sparse Matrix of class "dgCMatrix" - - 1 - -rhubarb 1.124517e-03 - -vegemite 1.444990e-01 - -apple 3.344565e-06 - - -</pre> - -<div class="toolFormTitle">rglasso log output</div> - -<pre> - -Loading required package: Matrix - -Loading required package: methods - -Loaded glmnet 1.9-8 - -Loaded lars 1.2 - -Warning messages: - -1: In if (class(larsres) == "try-error") { : - - the condition has length > 1 and only the first element will be used - -2: In plot.window(...) : "label" is not a graphical parameter - -3: In plot.xy(xy, type, ...) : "label" is not a graphical parameter - -4: In axis(side = side, at = at, labels = labels, ...) : - - "label" is not a graphical parameter - -5: In axis(side = side, at = at, labels = labels, ...) : - - "label" is not a graphical parameter - -6: In box(...) : "label" is not a graphical parameter - -7: In title(...) : "label" is not a graphical parameter - - -</pre> - -<div class="toolFormTitle">Other log output</div> - -<pre> - -## Toolfactory generated command line = Rscript - None None - -[1] Cox model: Got yvar= time,status cols 1 n preds= 3 forced in= 3 - -[2] Cox model: Got yvar= time,status cols 2 n preds= 3 forced in= 3 - -[1] Results preview: - - regulator partial_likelihood forced_in glmnet_model best_lambda - -rhubarb rhubarb 0.001180324 TRUE cox 0.01508333 - -vegemite vegemite 0.093827440 FALSE cox 0.01508333 - -[1] SessionInfo for this R session: - -R version 3.1.0 (2014-04-10) - -Platform: x86_64-unknown-linux-gnu (64-bit) - -locale: - - [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8 - - [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8 LC_PAPER=en_AU.UTF-8 LC_NAME=C - - [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C - -attached base packages: - -[1] methods stats graphics grDevices utils datasets base - -other attached packages: - -[1] lars_1.2 glmnet_1.9-8 Matrix_1.1-4 - -loaded via a namespace (and not attached): - -[1] grid_3.1.0 lattice_0.20-29 - - -</pre> - -<div class="toolFormTitle">All output files available for downloading</div> - -<div><table class="colored" cellpadding="3" cellspacing="3"><tr><th>Output File Name (click to view)</th><th>Size</th></tr> - -<tr><td><a href="cox_Coxglmnettest_glmnet_cvdeviance.pdf">cox_Coxglmnettest_glmnet_cvdeviance.pdf</a></td><td>6.4 KB</td></tr> -<tr class="odd_row"><td><a href="cox_Coxglmnettest_glmnetdev.pdf">cox_Coxglmnettest_glmnetdev.pdf</a></td><td>5.2 KB</td></tr> -<tr><td><a href="cox_rglasso.log">cox_rglasso.log</a></td><td>228 B</td></tr> -<tr class="odd_row"><td><a href="rglasso.Rscript">rglasso.Rscript</a></td><td>9.8 KB</td></tr> -<tr><td><a href="rglasso_error.log">rglasso_error.log</a></td><td>671 B</td></tr> -<tr class="odd_row"><td><a href="rglasso_runner.log">rglasso_runner.log</a></td><td>1.1 KB</td></tr> -</table></div><br/> -</div></body></html> -
--- a/test-data/coxlassotest_modelres.xls Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -regulator partial_likelihood forced_in glmnet_model best_lambda -rhubarb 0.0011803238918088 TRUE cox 0.0150833254995026 -vegemite 0.0938274399971897 FALSE cox 0.0150833254995026
--- a/tool_dependencies.xml Tue Oct 21 22:53:57 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -<?xml version="1.0"?> -<tool_dependency> - <package name="R_3_1_1" version="3.1.1"> - <repository changeset_revision="6ca21e466ef6" name="package_r_3_1_1" owner="fubar" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> - <package name="graphicsmagick" version="1.3.18"> - <repository changeset_revision="2fd4eb971ba5" name="package_graphicsmagick_1_3" owner="iuc" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> - <package name="ghostscript" version="9.10"> - <repository changeset_revision="9345d2740f0c" name="package_ghostscript_9_10" owner="devteam" prior_installation_required="True" toolshed="https://testtoolshed.g2.bx.psu.edu" /> - </package> - <package name="glmnet_lars_2_14" version="2.14"> - <install version="1.0"> - <actions> - <action type="setup_r_environment"> - <repository changeset_revision="6ca21e466ef6" name="package_r_3_1_1" owner="fubar" toolshed="https://testtoolshed.g2.bx.psu.edu"> - <package name="R_3_1_1" version="3.1.1" /> - </repository> - <package>https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/glmnet_1.9-8.tar.gz?raw=true</package> - <package>https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/lars_1.2.tar.gz?raw=true</package> - <package>https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/lattice_0.20-29.tar.gz?raw=true</package> - </action> - </actions> - </install> - <readme> - Yeee Haaa! - Lasso for Galaxy - </readme> - </package> -</tool_dependency>