# HG changeset patch # User fubar # Date 1446269667 14400 # Node ID 21b12c7c52e4a30ba905b684b64adc79b34246c1 # Parent 7d66bfa4fd56d8a6a748ca9d973e5ff913dced6e Fixes to paths in git for deps diff -r 7d66bfa4fd56 -r 21b12c7c52e4 .shed.yml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.shed.yml Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,9 @@ +categories: +- Statistics +description: glmnet and lars for galaxy +long_description: | + Allows hybrid lasso-ridge regression using glmnet +name: rglasso_1_9_8 +owner: fubar +remote_repository_url: https://github.com/galaxyproject/tools-iuc/tree/master/tools/rglasso +type: unrestricted diff -r 7d66bfa4fd56 -r 21b12c7c52e4 readme.rst --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/readme.rst Sat Oct 31 01:34:27 2015 -0400 @@ -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 + diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rgToolFactory.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rgToolFactory.py Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,649 @@ +# 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): + """ + 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 + 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.html = [] + self.cl = (opts.interpreter,self.sfile) + 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 + + + a tabular file + + reverse.py --script_path "$runMe" --interpreter "python" + --tool_name "reverse" --input_tab "$input1" --output_tab "$tab_file" + + + + + + + + + + + +**What it Does** + +Reverse the columns in a tabular file + + + + + +# 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() + + + + + + + """ + newXML=""" + %(tooldesc)s + %(command)s + + %(inputs)s + + + %(outputs)s + + + + %(script)s + + + %(tooltests)s + + %(help)s + + """ # needs a dict with toolname, toolid, interpreter, scriptname, command, inputs as a multi line string ready to write, outputs ditto, help ditto + + newCommand=""" + %(toolname)s.py --script_path "$runMe" --interpreter "%(interpreter)s" + --tool_name "%(toolname)s" %(command_inputs)s %(command_outputs)s + """ # may NOT be an input or htmlout + tooltestsTabOnly = """ + + + + + """ + tooltestsHTMLOnly = """ + + + + + """ + tooltestsBoth = """ + + + + + + """ + 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'] = '%s' % 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'] = ' \n' % self.inputFormats + else: + xdict['command_inputs'] = '' # assume no input - eg a random data generator + xdict['inputs'] = '' + xdict['inputs'] += ' \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'] += ' \n' + if self.opts.output_tab <> 'None': + xdict['command_outputs'] += ' --output_tab "$tab_file"' + xdict['outputs'] += ' \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 = """ + + + + + + + +
+ """ + galhtmlattr = """
This tool (%s) was generated by the Galaxy Tool Factory

""" + galhtmlpostfix = """
\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('
Galaxy Tool "%s" run at %s

' % (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('%s%s' % (fname,fname,sfsize)) + else: + fhtml.append('%s%s' % (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('
%s images and outputs
' % sectionname) + html.append('(Click on a thumbnail image to download the corresponding original PDF image)
') + ntogo = nacross # counter for table row padding with empty cells + html.append('
\n') + for i,paths in enumerate(ourpdfs): + fname,thumb = paths + s= """\n""" % (fname,thumb,fname,width,fname) + if ((i+1) % nacross == 0): + s += '\n' + ntogo = 0 + if i < (npdf - 1): # more to come + s += '' + ntogo = nacross + else: + ntogo -= 1 + html.append(s) + if html[-1].strip().endswith(''): + html.append('
Image called %s
\n') + else: + if ntogo > 0: # pad + html.append(' '*ntogo) + html.append('\n') + logt = open(logfname,'r').readlines() + logtext = [x for x in logt if x.strip() > ''] + html.append('
%s log output
' % sectionname) + if len(logtext) > 1: + html.append('\n
\n')
+                    html += logtext
+                    html.append('\n
\n') + else: + html.append('%s is empty
' % logfname) + if len(fhtml) > 0: + fhtml.insert(0,'
\n') + fhtml.append('
Output File Name (click to view)Size

') + html.append('
All output files available for downloading
\n') + html += fhtml # add all non-pdf files to the end of the display + else: + html.append('
### Error - %s returned no files - please confirm that parameters are sane
' % 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.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,cwd=self.opts.output_dir,env=my_env) + retval = p.wait() + 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() + else: + p = subprocess.Popen(self.cl,shell=False,env=my_env) + retval = p.wait() + if self.opts.make_HTML: + self.makeHtml() + return retval + + + def remove_me_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() + return retval + + +def main(): + u = """ + This is a Galaxy wrapper. It expects to be called by a special purpose tool.xml as (eg): + rgToolFactory.py --script_path "$scriptPath" --tool_name "foo" --interpreter "Rscript" + + The tool writes a script to a scriptPath using a configfile. + Some things in the script are templates. + The idea here is that depending on how this code is called, it uses the specified interpreter + to run a (hopefully correct) script/template. Optionally generates a clone of itself + which will run that same script/template as a toolshed repository tarball for uploading to a toolshed. + There's now an updated version which allows arbitrary parameters. + And so it goes. + """ + 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() + + diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rg_nri.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rg_nri.xml Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,633 @@ + + and other model improvement measures + + R_3_2_2 + graphicsmagick + ghostscript + glmnet_lars_3_2 + + + rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rg_NRI" + --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" + + + + + 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() + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**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 + + + + + + doi: 10.2215/​CJN.09590911 + + + + + diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rglasso_cox.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rglasso_cox.xml Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,910 @@ + + and cox regression using elastic net + + R_3_2_2 + graphicsmagick + ghostscript + glmnet_lars_3_2 + + + rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rglasso" + --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" + + + + 0 ) { + predcols = inmodel + xmat = as.matrix(x[,predcols]) + colnames(xmat) = preds + bestpred = getpredp_cox(x=xmat,time=y[,'time'],status=y[,'status'],id=cox_id, predict_at=predict_at) + pred = data.frame(responsep=bestpred, best_lambda=best_lambda,lamchoice=optLambda,alpha=alpha) + write.table(pred,out_pred_file,quote=FALSE, sep="\t",row.names=F) + } else { print.noquote('WARNING: No coefficients in selected model to predict with - no predictions made') } + } + if (debugOn) { + print.noquote(paste('best_lambda=',best_lambda,'saving cox respreds=',paste(names(coefs),collapse=','),'as predictors of survival. Coefs=',paste(coefs,collapse=','))) + } + p = try({res = data.frame(regulator=names(coefs),partial_likelihood=coefs,forced_in=pen,glmnet_model='cox',best_lambda=best_lambda, + lambdaChoice=optLambda,alpha=alpha)},T) + if (class(p) == "try-error") { + message(paste('@@@ unable to return a dataframe',p)) + sink() + return(NA) + } + print.noquote('@@@ Results preview:') + print.noquote(res,digits=5) + sink() + return(res) + +} + + +do_lasso = function(x=NA,y=NA,do_standard=T,debugOn=T,defaultFam="gaussian",optLambda='minLambda',descr='description', indx=1,target='target',sane=F, + alpha=0.9,nfold=10,penalty=c(),out_pred=F,out_pred_file='outpred',yvarname='yvar',id=c(),mdsplots=F) +{ + logf = file(paste(target,"rglasso.log",sep='_'), open = "a") + sink(logf,type = c("output", "message")) + res = NA + phe_is_bin = (length(unique(y)) == 2) + forcedin = paste(colnames(x)[which(penalty == 0)],collapse=',') + fam = "gaussian" + if (defaultFam %in% c("poisson","binomial","gaussian","multinomial")) fam=defaultFam + if (phe_is_bin == T) { + fam = "binomial" + } + print.noquote(paste('target=',target,'is binary=',phe_is_bin,'dim(x)=',paste(dim(x),collapse=','),'length(y)=',length(y),'force=',forcedin,'fam=',fam)) + standardize = do_standard + p = try({larsres = glmnet(x,y,family=fam,standardize=standardize,maxit=10000,alpha=alpha,penalty.factor=penalty) },T) + if (class(p) == "try-error") + { + print(paste('ERROR: unable to run glmnet for target',target,'error=',p)) + sink() + return(NA) + } + + mt = paste('Glmnet fraction deviance for',target) + outpdf = paste(target,'glmnetPath.pdf',sep='_') + pdf(outpdf) + plot(larsres,main=mt,label=T) + grid() + dev.off() + + outpdf = paste(target,'glmnetDeviance.pdf',sep='_') + + mt2 = paste('Glmnet lambda for',target) + + pdf(outpdf) + plot(larsres,xvar="lambda",main=mt2,label=T) + grid() + dev.off() + + larscv = NA + if (fam=="binomial") { + tmain = paste(target,'AUC') + outpdf = paste(target,'glmnetCV_AUC.pdf',sep='_') + p = try({larscv = ross.cv.glmnet(x=x,y=y,family=fam,type.measure='auc')},T) + } else { + tmain = paste(target,'CV MSE') + outpdf = paste(target,'glmnetCV_MSE.pdf',sep='_') + p = try({larscv = ross.cv.glmnet(x=x,y=y,family=fam,type.measure='mse')},T) + } + if (class(p) == "try-error") + { + print(paste('ERROR: unable to run cross validation for target',target,'error=',p)) + sink() + return(NA) + } + + pdf(outpdf) + plot(larscv,main=tmain) + grid() + dev.off() + + lse = larscv\$cvhits.1se + lmin = larscv\$cvhits.min + tot = lse + lmin + allhits = data.frame(pred=colnames(x),hits_lambda_1se = lse,hits_lambda_min = lmin) + nzhits = allhits[which(tot != 0),] + message('Total hit count for each predictor over all CV models (excluding zero count predictors)') + print.noquote(nzhits) + out_nz_file = paste(target,'cross_validation_model_counts.xls',sep='_') + write.table(nzhits,out_nz_file,quote=FALSE, sep="\t",row.names=F) + + ipenalty = c(0,penalty) + if (optLambda == 'lambda.min') { + best_lambda = larscv\$lambda.min + bestpred = as.matrix(coef(larscv, s = "lambda.min")) + } else { + best_lambda = larscv\$lambda.1se + bestpred = as.matrix(coef(larscv, s = "lambda.1se")) + } + inmodel = which(bestpred != 0) + coefs = bestpred[inmodel,1] + preds = rownames(bestpred)[inmodel] + iforced = ipenalty[inmodel] + forced = ! as.logical(iforced) + names(coefs) = preds + ncoef = length(coefs) - 1 + if (out_pred==T && fam=="binomial") + { + print.noquote(paste('Predicting',target,'probabilities from binomial glmnet at alpha',alpha,'and lambda',best_lambda)) + bestpred = predict.glmnet(larsres,s=best_lambda,newx=x,type="response") + bestpred = exp(bestpred)/(1+exp(bestpred)) + pred = data.frame(id=id,y=y,predp=as.vector(bestpred), best_lambda=best_lambda) + write.table(pred,out_pred_file,quote=FALSE, sep="\t",row.names=F) + } + 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,ncoef=ncoef, + best_lambda=best_lambda,lambdaChoice=optLambda,alpha=alpha),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=0.9,nfold=10,xcolnames=c(),ycolnames=c(),optLambda='minLambda', out_pred_file=NA, + descr="describe me",do_standard=F,defaultFam="gaussian",penalty=c(),out_pred=F,cox_id=c(),mdsplots=F,xfilt=0.95) +{ + logf = file("rglasso.log", open = "a") + sink(logf,type = c("output", "message")) + xdat = predvars + xm = data.matrix(xdat) + res = NULL + id = cox_id + depnames = ycolnames + ndep = length(depnames) + if (mdsplots==T) { + outpdfname = 'rglasso_x_in_sample_space_MDS.pdf' + p = try({ mdsPlot(xm,'measurements in sample space',groups=NA,outpdfname=outpdfname,transpose=T) },T) + if (class(p) == "try-error") + { + print.noquote(paste('Unable to produce predictors in sample space mds plot. Error:',p)) + } + outpdfname = 'rglasso_samples_in_x_space_MDS.pdf' + p = try({mdsPlot(xm,'samples in measurement space',groups=NA,outpdfname=outpdfname,transpose=F) },T) + if (class(p) == "try-error") + { + print.noquote(paste('Unable to produce samples in measurement space mds plot. Error:',p)) + } + } + ndat = nrow(xm) + cfracs = colSums(! is.na(xm))/ndat + keepme = (cfracs >= xfilt) + print.noquote(paste('Removing', sum(! keepme), 'xvars with more than',xfilt,'fraction missing')) + vars = apply(xm,2,var,na.rm=T) + xm = xm[,keepme] + 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] + } + if (fam == "binomial") {y = as.factor(y)} + x = xm + id = cox_id + if (fam != "cox") { + ok = complete.cases(x,y) + if (sum(! ok) > 0) { + message(paste('@@@ Removing',sum(! ok),'cases with missing y of',length(y),'@@@')) + y = y[(ok)] + x = x[(ok),] + id = id[(ok)] + } + } + ok = complete.cases(y) + 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 { + if (i == 1) { outpred = out_pred_file + } else { + outpred = paste(target,'predicted_output.xls') + } + regres = do_lasso(x=x,y=y,do_standard=do_standard,debugOn=debugOn,defaultFam=defaultFam,optLambda=optLambda,out_pred_file=outpred, + descr=descr,indx=i,target=target,alpha=alpha,nfold=nfold,penalty=penalty,out_pred=out_pred,yvarname=target,id=id,mdsplots=mdsplots) + if (! is.na(regres)) { res = rbind(res,regres) } + } + } + print.noquote('@@@ Results preview:') + print.noquote(res,digits=5) + 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',do_standard=do_standard,defaultFam="gaussian",alpha=0.05) + return(regres_out) +} +]]> +options(width=512) +options(digits=5) +alpha = $alpha +nfold = $nfold +optLambda = "$optLambda" +Out_Dir = "$html_file.files_path" +Input = "$input1" +indat = read.table(Input,head=T,sep='\t') +datcols = colnames(indat) +myTitle = "$title" +outtab = "$model_file" +do_standard = as.logical("$do_standard") +mdsplots = as.logical("$mdsplots") +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() +penalties = rep(1,length(datcols)) +forced_in = NA + +logxform = "$logxform_cols" +if (logxform != "None") { + logxform_cols = as.numeric(strsplit(logxform,",")[[1]]) + if (length(logxform_cols) > 0) { + small = 1e-10 + sset = indat[,logxform_cols] + zeros = which(sset==0,arr.ind=T) + nz = nrow(zeros) + if (nz > 0) { + message(paste('Log transforming encountered',nz,'zeros - added 1e-10')) + sset[zeros] = sset[zeros] + small + lset = log(sset) + indat[,logxform_cols] = lset + } + } +} +if (force_xvar_cols_in != "None") +{ + 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)] + penalties[force_xvar_cols] = 0 +} +penalty = penalties[xvar_cols] +forcedin = paste(datcols[which(penalties == 0)],collapse=',') +cox_id_col = NA +cox_id = NA + +message(paste('@@@ Using alpha =',alpha,'for all models')) +x = indat[,xvar_cols] +nx = nrow(x) +cx = ncol(x) +message(paste('@@@@ Input has',nx,'samples and',cx,'predictors')) +if (cx > nx) { +message('@@@ WARNING: Models will have more variables than cases so glmnet will likely return one of many possible solutions! Please DO NOT expect reliable results - glmnet is clever but not magical @@@') +} + +xcolnames = datcols[xvar_cols] + +if (file.exists(Out_Dir) == F) dir.create(Out_Dir) +out_full = F +out_full_file = NA +out_pred_file = "" +out_pred = as.logical("$model.output_pred") + +#if $model.fam == "binomial" or $model.fam == "cox": + cox_id_col = $model.cox_id + cox_id = indat[,cox_id_col] + if (out_pred == T) { + out_pred_file="$output_pred_file" + rownames(x) = cox_id + } +#end if +#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="$output_full_file" } + yvar_cols = c(cox_time,cox_status) + ycolnames = c('time','status') + istat = as.double(indat[,cox_status]) + itime = as.double(indat[,cox_time]) + predict_at = quantile(itime) + if ("$model.predict_at" > "") + { + pa = "$model.predict_at" + predict_at = as.numeric(strsplit(pa,",")[[1]]) + } + 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) + print.noquote(paste('@@@ Cox model will predict yvar=',datcols[cox_status],'using cols=',paste(xcolnames,collapse=','),'n preds=',length(xcolnames), + 'forced in=',forcedin)) + regres_out = dolasso_cox(x=x,y=y,debugOn=F,maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames,optLambda=optLambda,out_full=out_full,out_full_file=out_full_file, + out_pred=out_pred,out_pred_file=out_pred_file,cox_id=cox_id,descr=myTitle,do_standard=do_standard,alpha=alpha,penalty=penalty,predict_at=predict_at,mdsplots=mdsplots) +#else: + yvar_cols = "$model.yvar_cols" + yvar_cols = as.numeric(strsplit(yvar_cols,",")[[1]]) + ycolnames = datcols[yvar_cols] + print.noquote(paste('@@@',fam,'model will predict yvar=',paste(ycolnames,collapse=','),'using cols=',paste(xcolnames,collapse=','),'n preds=',length(xcolnames), + 'forced in=',forcedin)) + y = data.matrix(indat[,yvar_cols]) + print.noquote(paste('Model will use',fam,'link function to predict yvar=',paste(ycolnames,collapse=','),'n preds=',length(xcolnames),'forced in=',forcedin)) + regres_out = dolasso_generic(predvars=x,depvars=y,debugOn=F, maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames,optLambda=optLambda,out_pred_file=out_pred_file, + descr=myTitle,do_standard=do_standard,defaultFam=fam,alpha=alpha,penalty=penalty,out_pred=out_pred,cox_id=cox_id,mdsplots=mdsplots) +#end if + +write.table(format(regres_out,digits=5),outtab,quote=FALSE, sep="\t",row.names=F) +print.noquote('@@@ SessionInfo for this R session:') +sessionInfo() +warnings() + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + l + + + + l + + + + + + + + + + model['output_full'] == 'T' + + + model['output_pred'] == 'T' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**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. + +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. + +**Warning about the tyrany of dimensionality** + +Yes, this package will select 'optimal' models even when you (optimistically) supply more predictors than you have cases. +The model returned is unlikely to represent the only informative regularisation path through your data - if you run repeatedly with +exactly the same settings, you will probably see many different models being selected. +This is not a software bug - the real problem is that you just don't have enough information in your data. + +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. Predictors 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. + +**Predicted event probabilities for Cox and Logistic models** + +If you want to compare (eg) two competing clinical predictions, there's a companion generic NRI tool +for predicted event probabilities. Estimates dozens of measures of improvement in prediction. Currently only works for identical id subjects +but can probably be extended to independent sample predictions. + +Given a model, we can generate a predicted p (for status 1) in binomial or cox frameworks so models can be evaluated in terms of NRI. +Of course, estimates are likely substantially inflated over 'real world' performance by being estimated from the same sample - but you probably +already knew that since you were smart enough to reach this far down into the on screen help. The author salutes you, intrepid reader! + +It may seem an odd thing to do, but we can predict p for an event for each subject from our original data, given a parsimonious model. Doing +this for two separate models (eg, forcing in an additional known explanatory measurement to the new model) allows comparison of the two models +predicted status for each subject, or the same model in independent populations to see how badly it does + +**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 + + + + +@Article{Friedman2010, title = {Regularization Paths for Generalized Linear Models via Coordinate Descent}, + author = {Jerome Friedman and Trevor Hastie and Robert Tibshirani}, + journal = {Journal of Statistical Software}, + year = {2010}, + volume = {33}, + number = {1}, + pages = {1--22}, + url = {http://www.jstatsoft.org/v33/i01/} + } + + +10.1093/bioinformatics/bts573 + + + diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_coxlassotest_glmnet_cvdeviance.pdf Binary file test-data/cox_coxlassotest_glmnet_cvdeviance.pdf has changed diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_coxlassotest_glmnetdev.pdf Binary file test-data/cox_coxlassotest_glmnetdev.pdf has changed diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_test.xls --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/cox_test.xls Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,502 @@ +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 + diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/coxlassotest.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/coxlassotest.html Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,167 @@ + + + + + + + + +
+ +
Galaxy Tool "rglasso" run at 18/02/2015 22:06:08

+
cox images and outputs
+(Click on a thumbnail image to download the corresponding original PDF image)
+
+ + + + + + +
Image called cox_Coxglmnettest_glmnet_cvdeviance.pdfImage called cox_Coxglmnettest_glmnetdev.pdf
+ +
cox log output
+ +
+
+[1] @@@ Cox model will be predicted at times = 3.724234,1325.40516225,2373.090418,3650.1730825,4989.405682
+
+[1] Times each predictor was selected in CV models (excluding zero count predictors)
+
+         hits_lambda_1se hits_lambda_min
+
+rhubarb               10              10
+
+vegemite               6               6
+
+[1] @@@ Results preview:
+
+        regulator partial_likelihood forced_in glmnet_model best_lambda lambdaChoice alpha
+
+rhubarb   rhubarb          0.0012215      TRUE          cox     0.04197   lambda.1se  0.95
+
+
+
+ +
rglasso log output
+ +
+
+Loading required package: Matrix
+
+Loading required package: methods
+
+Loaded glmnet 1.9-8
+
+Loaded lars 1.2
+
+Loading required package: splines
+
+Warning messages:
+
+1: In if (is.na(predict_at)) { :
+
+  the condition has length > 1 and only the first element will be used
+
+2: In if (class(p) == "try-error") { :
+
+  the condition has length > 1 and only the first element will be used
+
+3: In plot.window(...) : "label" is not a graphical parameter
+
+4: In plot.xy(xy, type, ...) : "label" is not a graphical parameter
+
+5: In axis(side = side, at = at, labels = labels, ...) :
+
+  "label" is not a graphical parameter
+
+6: In axis(side = side, at = at, labels = labels, ...) :
+
+  "label" is not a graphical parameter
+
+7: In box(...) : "label" is not a graphical parameter
+
+8: In title(...) : "label" is not a graphical parameter
+
+
+
+ +
Other log output
+ +
+
+## Toolfactory running rglasso as Rscript script
+
+[1] @@@ Using alpha = 0.95 for all models
+
+[1] @@@@ Input has 500 samples and 3 predictors
+
+[1] @@@ Cox model will predict yvar= status using cols= rhubarb,vegemite,apple n preds= 3 forced in= rhubarb
+
+[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     LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8    LC_PAPER=en_AU.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       
+
+attached base packages:
+
+[1] splines   methods   stats     graphics  grDevices utils     datasets  base     
+
+other attached packages:
+
+[1] pec_2.4.4       survival_2.37-7 lars_1.2        glmnet_1.9-8    Matrix_1.1-5   
+
+loaded via a namespace (and not attached):
+
+[1] codetools_0.2-10 foreach_1.4.2    grid_3.1.0       iterators_1.0.7  lattice_0.20-29  lava_1.3         prodlim_1.5.1   
+
+Warning messages:
+
+1: In if (is.na(predict_at)) { ... :
+
+  the condition has length > 1 and only the first element will be used
+
+2: In if (class(p) == "try-error") { ... :
+
+  the condition has length > 1 and only the first element will be used
+
+3: In plot.window(...) : "label" is not a graphical parameter
+
+4: In plot.xy(xy, type, ...) : "label" is not a graphical parameter
+
+5: In axis(side = side, at = at, labels = labels, ...) :
+
+  "label" is not a graphical parameter
+
+6: In axis(side = side, at = at, labels = labels, ...) :
+
+  "label" is not a graphical parameter
+
+7: In box(...) : "label" is not a graphical parameter
+
+8: In title(...) : "label" is not a graphical parameter
+
+
+
+ +
All output files available for downloading
+ +
+ + + + + + + + +
Output File Name (click to view)Size
cox_Coxglmnettest_glmnet_cvdeviance.pdf6.4 KB
cox_Coxglmnettest_glmnetdev.pdf5.3 KB
cox_cross_validation_model_counts.xls42 B
cox_rglasso.log522 B
rglasso.Rscript21.5 KB
rglasso_error.log803 B
rglasso_runner.log1.7 KB

+
+ diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/coxlassotest_modelres.xls --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/coxlassotest_modelres.xls Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,2 @@ +regulator partial_likelihood forced_in glmnet_model best_lambda lambdaChoice alpha +rhubarb 0.0012215 TRUE cox 0.04197 lambda.1se 0.95 diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/genTest.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/genTest.R Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,10 @@ +ids=c(1:50) +io1 = rep(c(0,0,0,0,1),10) +ip2 = runif(50)+0.1 +ip2[which(ip2>1.0)] = 1.0 +ip1 = runif(50)+0.05 +ip1[which(ip1>1.0)] = 1.0 +df=data.frame(id=ids,input1_observed=io1,input1_predicted=ip1,input2_predicted=ip2) +fout='test-data/nri_test1.xls' +write.table(df,file=fout, quote=FALSE, sep="\t",row.names=F) +# planemo test --job_output_files /home/rlazarus/tmp --test_output /home/rlazarus/tmp/startest/foo.html --update_test_data --galaxy_root /home/rlazarus/galaxy rg_nri.xml diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/nri_test1.xls --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/nri_test1.xls Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,51 @@ +id input1_observed input1_predicted input2_predicted +1 0 0.431491391919553 0.786344331596047 +2 0 0.235500678373501 0.37771512279287 +3 0 0.458084875205532 0.343597211781889 +4 0 0.400509147858247 0.50323077281937 +5 1 0.463376078708097 0.28904068809934 +6 0 0.638350193109363 0.745497886231169 +7 0 0.909975615050644 1 +8 0 0.0712537909392267 0.613464191928506 +9 0 0.475019342219457 0.757945845043287 +10 1 0.146289554610848 0.948540200153366 +11 0 0.947521302569658 0.452743433788419 +12 0 0.512624199455604 0.107415618747473 +13 0 0.260238706320524 0.165854234388098 +14 0 0.252348658815026 0.866530107706785 +15 1 0.552841167151928 0.781692814920098 +16 0 0.203356911847368 0.397802447108552 +17 0 0.138653858751059 0.966900746012107 +18 0 0.237485485337675 0.126699626818299 +19 0 0.317045996850356 0.646952681383118 +20 1 0.271802965737879 0.828731852956116 +21 0 0.522744563594461 0.421856845868751 +22 0 0.819749742234126 0.628745510149747 +23 0 0.27876661689952 0.826974717108533 +24 0 0.077877387823537 0.12339704008773 +25 1 0.911161796143279 1 +26 0 0.567042255960405 1 +27 0 0.223712538788095 0.738010874623433 +28 0 0.306199585506693 0.388159935269505 +29 0 0.305726604955271 0.618934138352051 +30 1 0.919456032756716 1 +31 0 0.493467168603092 0.470883537689224 +32 0 0.626042458415031 0.753118774993345 +33 0 0.739080456737429 0.812446624180302 +34 0 0.786149174161255 0.768807894224301 +35 1 0.65491428244859 0.855509085487574 +36 0 0.62289561489597 0.585996081028134 +37 0 0.297479379642755 0.357051227381453 +38 0 0.796507286094129 0.775113077275455 +39 0 0.581927774194628 0.569900369411334 +40 1 0.76229884494096 0.877794308261946 +41 0 0.845583727071062 0.804134764615446 +42 0 0.79128703516908 1 +43 0 0.0507813768461347 0.97111975508742 +44 0 0.49284805515781 0.101141295302659 +45 1 0.665627157501876 0.121852198988199 +46 0 0.059045681450516 0.996109608234838 +47 0 0.0814651515800506 0.852028644317761 +48 0 0.162336849560961 1 +49 0 0.577456943178549 0.343083152920008 +50 1 0.905623372131959 0.651731353392825 diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/nri_test1_out.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/nri_test1_out.html Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,34 @@ + + + + + + + + +
+ +
Galaxy Tool "rg_NRI" run at 08/01/2015 16:12:57

+
rg log output
+ +
+
+Error in library("e1071") : there is no package called ‘e1071’
+
+Execution halted
+
+
+
+ +
Other log output
+/tmp/tmpq72Dni/job_working_directory/000/2/dataset_2_files/rg_NRI_runner.log is empty
+
All output files available for downloading
+ +
+ + + + +
Output File Name (click to view)Size
rg_NRI.Rscript18.1 KB
rg_NRI_error.log84 B
rg_NRI_runner.log48 B

+
+ diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/nri_test1_out.xls diff -r 7d66bfa4fd56 -r 21b12c7c52e4 tool_dependencies.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Sat Oct 31 01:34:27 2015 -0400 @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/RColorBrewer_1.1-2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/e1071_1.6-7.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/Rcpp_0.12.1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/plyr_1.8.3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/digest_0.6.8.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/gtable_0.1.2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/stringi_1.0-1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/magrittr_1.5.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/stringr_1.0.0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/reshape2_1.4.1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/dichromat_2.0-0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/colorspace_1.2-6.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/munsell_0.4.2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/labeling_0.3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/scales_0.3.0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/proto_0.3-10.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/ggplot2_1.0.1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/minqa_1.2.4.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/nloptr_1.0.4.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/lme4_1.1-10.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/pbkrtest_0.4-2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/SparseM_1.7.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/MatrixModels_0.4-1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/quantreg_5.19.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/car_2.1-0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/caret_6.0-58.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/iterators_1.0.8.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/foreach_1.4.3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/pROC_1.8.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/Formula_1.2-1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/latticeExtra_0.6-26.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/acepack_1.3-3.3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/gridExtra_2.0.0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/Hmisc_3.17-0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/pracma_1.8.6.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/survival_2.38-3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/lars_1.2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/glmnet_2.0-2.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/numDeriv_2014.2-1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/lava_1.4.1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/prodlim_1.5.5.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/polspline_1.1.12.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/mvtnorm_1.0-3.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/TH.data_1.0-6.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/zoo_1.7-12.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/sandwich_2.3-4.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/multcomp_1.4-1.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/rms_4.4-0.tar.gz?raw=true + https://github.com/fubar2/galaxy_tool_source/blob/master/RELEASE_2_14/pec_2.4.7.tar.gz?raw=true + + + + + Yee-Haw! Lasso for Galaxy! + + +