# HG changeset patch
# User richjp
# Date 1413201436 14400
# Node ID 810bc97ddeb2cc94eb8c3fc07a99282209b2f411
# Parent c3956b7cb71f394b2a355c58ea34b89c7c9d261d
Uploaded
diff -r c3956b7cb71f -r 810bc97ddeb2 Mydynamicscript/Mydynamicscript.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Mydynamicscript/Mydynamicscript.py Mon Oct 13 07:57:16 2014 -0400
@@ -0,0 +1,742 @@
+# Mydynamicscript/Mydynamicscript.py - a self annotated version of rgToolFactory.py generated by running rgToolFactory.py
+# to make a new Galaxy tool called Mydynamicscript
+# User r.poole@ucl.ac.uk at 13/10/2014 12:55:18
+# 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
+#
+# August 2014
+# merged John Chilton's citation addition and ideas from Marius van den Beek to enable arbitrary
+# data types for input and output - thanks!
+#
+# march 2014
+# had to remove dependencies because cross toolshed dependencies are not possible - can't pre-specify a toolshed url for graphicsmagick and ghostscript
+# grrrrr - night before a demo
+# added dependencies to a tool_dependencies.xml if html page generated so generated tool is properly portable
+#
+# added ghostscript and graphicsmagick as dependencies
+# fixed a wierd problem where gs was trying to use the new_files_path from universe (database/tmp) as ./database/tmp
+# errors ensued
+#
+# 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 = 'V001.1 March 2014'
+verbose = False
+debug = False
+toolFactoryURL = 'https://bitbucket.org/fubar/galaxytoolfactory'
+
+# if we do html we need these dependencies specified in a tool_dependencies.xml file and referred to in the generated
+# tool xml
+toolhtmldepskel = """
+
+
+
+
+
+
+
+
+ %s
+
+
+"""
+
+protorequirements = """
+ ghostscript
+ graphicsmagick
+ """
+
+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
+
+
+def parse_citations(citations_text):
+ """
+ """
+ citations = [c for c in citations_text.split("**ENTRY**") if c.strip()]
+ citation_tuples = []
+ for citation in citations:
+ if citation.startswith("doi"):
+ citation_tuples.append( ("doi", citation[len("doi"):].strip() ) )
+ else:
+ citation_tuples.append( ("bibtex", citation[len("bibtex"):].strip() ) )
+ return citation_tuples
+
+
+class ScriptRunner:
+ """class is a wrapper for an arbitrary script
+ """
+
+ def __init__(self,opts=None,treatbashSpecial=True):
+ """
+ cleanup inputs, setup some outputs
+
+ """
+ self.useGM = cmd_exists('gm')
+ self.useIM = cmd_exists('convert')
+ self.useGS = cmd_exists('gs')
+ self.temp_warned = False # we want only one warning if $TMP not set
+ self.treatbashSpecial = treatbashSpecial
+ if opts.output_dir: # simplify for the tool tarball
+ os.chdir(opts.output_dir)
+ self.thumbformat = 'png'
+ self.opts = opts
+ self.toolname = re.sub('[^a-zA-Z0-9_]+', '', opts.tool_name) # a sanitizer now does this but..
+ self.toolid = self.toolname
+ self.myname = sys.argv[0] # get our name because we write ourselves out as a tool later
+ self.pyfile = self.myname # crude but efficient - the cruft won't hurt much
+ self.xmlfile = '%s.xml' % self.toolname
+ s = open(self.opts.script_path,'r').readlines()
+ s = [x.rstrip() for x in s] # remove pesky dos line endings if needed
+ self.script = '\n'.join(s)
+ fhandle,self.sfile = tempfile.mkstemp(prefix=self.toolname,suffix=".%s" % (opts.interpreter))
+ tscript = open(self.sfile,'w') # use self.sfile as script source for Popen
+ tscript.write(self.script)
+ tscript.close()
+ self.indentedScript = '\n'.join([' %s' % html_escape(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)
+ if opts.output_dir: # may not want these complexities
+ self.tlog = os.path.join(self.opts.output_dir,"%s_runner.log" % self.toolname)
+ art = '%s.%s' % (self.toolname,opts.interpreter)
+ artpath = os.path.join(self.opts.output_dir,art) # need full path
+ artifact = open(artpath,'w') # use self.sfile as script source for Popen
+ artifact.write(self.script)
+ artifact.close()
+ self.cl = []
+ self.html = []
+ a = self.cl.append
+ a(opts.interpreter)
+ if self.treatbashSpecial and opts.interpreter in ['bash','sh']:
+ a(self.sfile)
+ else:
+ a('-') # stdin
+ a(opts.input_tab)
+ a(opts.output_tab)
+ self.outputFormat = self.opts.output_format
+ self.inputFormats = self.opts.input_formats
+ 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
+%(requirements)s
+
+%(command)s
+
+
+%(inputs)s
+
+
+%(outputs)s
+
+
+
+%(script)s
+
+
+
+%(tooltests)s
+
+
+
+%(help)s
+
+
+
+ %(citations)s
+ 10.1093/bioinformatics/bts573
+
+ """ # 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 - appended later
+ tooltestsTabOnly = """
+
+
+
+
+
+
+
+
+ """
+ tooltestsHTMLOnly = """
+
+
+
+
+
+
+
+
+ """
+ tooltestsBoth = """
+
+
+
+
+
+
+
+
+ """
+ xdict = {}
+ xdict['outputFormat'] = self.outputFormat
+ xdict['inputFormats'] = self.inputFormats
+ xdict['requirements'] = ''
+ if self.opts.make_HTML:
+ if self.opts.include_dependencies == "yes":
+ xdict['requirements'] = protorequirements
+ 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:
+ helptext = open(self.opts.help_text,'r').readlines()
+ helptext = [html_escape(x) for x in helptext] # must html escape here too - thanks to Marius van den Beek
+ xdict['help'] = ''.join([x for x in helptext])
+ else:
+ xdict['help'] = 'Please ask the tool author (%s) for help as none was supplied at tool generation\n' % (self.opts.user_email)
+ if self.opts.citations:
+ citationstext = open(self.opts.citations,'r').read()
+ citation_tuples = parse_citations(citationstext)
+ citations_xml = ""
+ for citation_type, citation_content in citation_tuples:
+ citation_xml = """%s """ % (citation_type, html_escape(citation_content))
+ citations_xml += citation_xml
+ xdict['citations'] = citations_xml
+ else:
+ xdict['citations'] = ""
+ coda = ['**Script**','Pressing execute will run the following code over your input file and generate some outputs in your history::']
+ coda.append('\n')
+ coda.append(self.indentedScript)
+ coda.append('\n**Attribution**\nThis Galaxy tool was created by %s at %s\nusing the Galaxy Tool Factory.\n' % (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\n')
+ 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'
+ else:
+ xdict['command_outputs'] += ' --output_dir "./"'
+ if self.opts.output_tab <> 'None':
+ xdict['command_outputs'] += ' --output_tab "$tab_file"'
+ xdict['outputs'] += ' \n' % self.outputFormat
+ 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)
+ tdir = self.toolname
+ os.mkdir(tdir)
+ self.makeXML()
+ if self.opts.make_HTML:
+ if self.opts.help_text:
+ hlp = open(self.opts.help_text,'r').read()
+ else:
+ hlp = 'Please ask the tool author for help as none was supplied at tool generation\n'
+ if self.opts.include_dependencies:
+ tooldepcontent = toolhtmldepskel % hlp
+ depf = open(os.path.join(tdir,'tool_dependencies.xml'),'w')
+ depf.write(tooldepcontent)
+ depf.write('\n')
+ depf.close()
+ 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'))
+ outpif = '%s.py' % self.toolname # new name
+ outpiname = os.path.join(tdir,outpif) # 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' % (outpiname,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)
+ hlog = os.path.join(self.opts.output_dir,"compress_%s.txt" % os.path.basename(inpdf))
+ sto = open(hlog,'a')
+ our_env = os.environ.copy()
+ our_tmp = our_env.get('TMP',None)
+ if not our_tmp:
+ our_tmp = our_env.get('TEMP',None)
+ if not (our_tmp and os.path.exists(our_tmp)):
+ newtmp = os.path.join(self.opts.output_dir,'tmp')
+ try:
+ os.mkdir(newtmp)
+ except:
+ sto.write('## WARNING - cannot make %s - it may exist or permissions need fixing\n' % newtmp)
+ our_env['TEMP'] = newtmp
+ if not self.temp_warned:
+ sto.write('## WARNING - no $TMP or $TEMP!!! Please fix - using %s temporarily\n' % newtmp)
+ self.temp_warned = True
+ 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)
+ 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:
+ 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 = """
"""
+ 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('
\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,'Output File Name (click to view) Size \n')
+ fhtml.append('
')
+ 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!
+ """
+ if self.treatbashSpecial and self.opts.interpreter in ['bash','sh']:
+ retval = self.runBash()
+ else:
+ if self.opts.output_dir:
+ ste = open(self.elog,'w')
+ sto = open(self.tlog,'w')
+ sto.write('## Toolfactory generated command line = %s\n' % ' '.join(self.cl))
+ sto.flush()
+ p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=ste,stdin=subprocess.PIPE,cwd=self.opts.output_dir)
+ else:
+ p = subprocess.Popen(self.cl,shell=False,stdin=subprocess.PIPE)
+ p.stdin.write(self.script)
+ p.stdin.close()
+ retval = p.wait()
+ if self.opts.output_dir:
+ sto.close()
+ ste.close()
+ err = open(self.elog,'r').readlines()
+ if retval <> 0 and err: # problem
+ print >> sys.stderr,err
+ if self.opts.make_HTML:
+ self.makeHtml()
+ return retval
+
+ def runBash(self):
+ """
+ cannot use - for bash so use self.sfile
+ """
+ if self.opts.output_dir:
+ s = '## Toolfactory generated command line = %s\n' % ' '.join(self.cl)
+ sto = open(self.tlog,'w')
+ sto.write(s)
+ sto.flush()
+ p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=sto,cwd=self.opts.output_dir)
+ else:
+ p = subprocess.Popen(self.cl,shell=False)
+ retval = p.wait()
+ if self.opts.output_dir:
+ sto.close()
+ if self.opts.make_HTML:
+ self.makeHtml()
+ return retval
+
+
+def main():
+ u = """
+ This is a Galaxy wrapper. It expects to be called by a special purpose tool.xml as:
+ rgBaseScriptWrapper.py --script_path "$scriptPath" --tool_name "foo" --interpreter "Rscript"
+
+ """
+ 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='./')
+ a('--output_html',default=None)
+ a('--input_tab',default="None")
+ a('--input_formats',default="tabular,text")
+ a('--output_tab',default="None")
+ a('--output_format',default="tabular")
+ 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('--citations',default=None)
+ a('--tool_desc',default=None)
+ a('--new_tool',default=None)
+ a('--tool_version',default=None)
+ a('--include_dependencies',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 c3956b7cb71f -r 810bc97ddeb2 Mydynamicscript/Mydynamicscript.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Mydynamicscript/Mydynamicscript.xml Mon Oct 13 07:57:16 2014 -0400
@@ -0,0 +1,1297 @@
+
+
+
+ ghostscript
+ graphicsmagick
+
+
+
+ Mydynamicscript.py --script_path "$runMe" --interpreter "Rscript"
+ --tool_name "Mydynamicscript" --input_tab "$input1" --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes"
+
+
+
+
+
+
+
+
+
+
+
+
+#VDM-mapping script
+#input file: Q200 variants no wildtype
+ourargs = commandArgs(TRUE)
+inf = ourargs[1]
+outf = ourargs[2]
+Q200var = read.table(inf)
+
+
+het <- Q200var[grep("AF=0.5",Q200var\$V8), ]
+hom <- Q200var[grep("AF=1.0",Q200var\$V8), ]
+
+get_hom1 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr1")
+ het <- subset(x2,x2\$V1=="chr1")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom2 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr2")
+ het <- subset(x2,x2\$V1=="chr2")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom3 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr3")
+ het <- subset(x2,x2\$V1=="chr3")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom4 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr4")
+ het <- subset(x2,x2\$V1=="chr4")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom5 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr5")
+ het <- subset(x2,x2\$V1=="chr5")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom6 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr6")
+ het <- subset(x2,x2\$V1=="chr6")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom7 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr7")
+ het <- subset(x2,x2\$V1=="chr7")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom8 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr8")
+ het <- subset(x2,x2\$V1=="chr8")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom9 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr9")
+ het <- subset(x2,x2\$V1=="chr9")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom10 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr10")
+ het <- subset(x2,x2\$V1=="chr10")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom11 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr11")
+ het <- subset(x2,x2\$V1=="chr11")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom12 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr12")
+ het <- subset(x2,x2\$V1=="chr12")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom13 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr13")
+ het <- subset(x2,x2\$V1=="chr13")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom14 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr14")
+ het <- subset(x2,x2\$V1=="chr14")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom15 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr15")
+ het <- subset(x2,x2\$V1=="chr15")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom16 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr16")
+ het <- subset(x2,x2\$V1=="chr16")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom17 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr17")
+ het <- subset(x2,x2\$V1=="chr17")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom18 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr18")
+ het <- subset(x2,x2\$V1=="chr18")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom19 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr19")
+ het <- subset(x2,x2\$V1=="chr19")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom20 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr20")
+ het <- subset(x2,x2\$V1=="chr20")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom21 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr21")
+ het <- subset(x2,x2\$V1=="chr21")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom22 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr22")
+ het <- subset(x2,x2\$V1=="chr22")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom23 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr23")
+ het <- subset(x2,x2\$V1=="chr23")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom24 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr24")
+ het <- subset(x2,x2\$V1=="chr24")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom25 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr25")
+ het <- subset(x2,x2\$V1=="chr25")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+h1 <- get_hom1(hom,het)
+h2 <- get_hom2(hom,het)
+h3 <- get_hom3(hom,het)
+h4 <- get_hom4(hom,het)
+h5 <- get_hom5(hom,het)
+h6 <- get_hom6(hom,het)
+h7 <- get_hom7(hom,het)
+h8 <- get_hom8(hom,het)
+h9 <- get_hom9(hom,het)
+h10 <- get_hom10(hom,het)
+h11 <- get_hom11(hom,het)
+h12 <- get_hom12(hom,het)
+h13 <- get_hom13(hom,het)
+h14 <- get_hom14(hom,het)
+h15 <- get_hom15(hom,het)
+h16 <- get_hom16(hom,het)
+h17 <- get_hom17(hom,het)
+h18 <- get_hom18(hom,het)
+h19 <- get_hom19(hom,het)
+h20 <- get_hom20(hom,het)
+h21 <- get_hom21(hom,het)
+h22 <- get_hom22(hom,het)
+h23 <- get_hom23(hom,het)
+h24 <- get_hom24(hom,het)
+h25 <- get_hom25(hom,het)
+
+
+
+#Extracting maxima (VDM-mapping)
+
+zero_crossings_maxima <- function(x,y){
+
+ dy <- diff(y)
+ xn <- x[2:length(x)]
+
+ n <- length(dy)
+ t1 <- dy[1:n-1];
+ t2 <- dy[2:n];
+ tt <- t1*t2
+ ind <- which(tt < 0)
+
+ # print out the maxima
+ d2y <- diff(dy)
+ vind <- c()
+ vy <- c()
+ vxn <- c()
+
+ xnn <- x[3:length(x)]
+ for(i in 1:length(ind)){
+ ### cat( ind[i], d2y[ind[i]], '\n')
+ if( d2y[ind[i]] < 0 ){
+ ### cat(ind[i], y[ind[i]], xn[ind[i]], "\n" )
+ vind <- c(vind, ind[i])
+ vy <- c(vy, y[ind[i]])
+ vxn <- c(vxn, xn[ind[i]])
+ }
+ }
+
+ maxima <- data.frame(ind=vind, y=vy, x=vxn)
+ maxima <- maxima[order(maxima\$y),]
+
+ #print( maxima )
+
+ results <- list()
+ results\$dy <- dy
+ results\$xn <- xn
+ results\$maxima <- maxima[ order(maxima\$y, decreasing=T), ]
+
+ return(results)
+}
+
+res1 <- zero_crossings_maxima(h1\$x,h1\$ratio)
+res2 <- zero_crossings_maxima(h2\$x,h2\$ratio)
+res3 <- zero_crossings_maxima(h3\$x,h3\$ratio)
+res4 <- zero_crossings_maxima(h4\$x,h4\$ratio)
+res5 <- zero_crossings_maxima(h5\$x,h5\$ratio)
+res6 <- zero_crossings_maxima(h6\$x,h6\$ratio)
+res7 <- zero_crossings_maxima(h7\$x,h7\$ratio)
+res8 <- zero_crossings_maxima(h8\$x,h8\$ratio)
+res9 <- zero_crossings_maxima(h9\$x,h9\$ratio)
+res10 <- zero_crossings_maxima(h10\$x,h10\$ratio)
+res11 <- zero_crossings_maxima(h11\$x,h11\$ratio)
+res12 <- zero_crossings_maxima(h12\$x,h12\$ratio)
+res13 <- zero_crossings_maxima(h13\$x,h13\$ratio)
+res14 <- zero_crossings_maxima(h14\$x,h14\$ratio)
+res15 <- zero_crossings_maxima(h15\$x,h15\$ratio)
+res16 <- zero_crossings_maxima(h16\$x,h16\$ratio)
+res17 <- zero_crossings_maxima(h17\$x,h17\$ratio)
+res18 <- zero_crossings_maxima(h18\$x,h18\$ratio)
+res19 <- zero_crossings_maxima(h19\$x,h19\$ratio)
+res20 <- zero_crossings_maxima(h20\$x,h20\$ratio)
+res21 <- zero_crossings_maxima(h21\$x,h21\$ratio)
+res22 <- zero_crossings_maxima(h22\$x,h22\$ratio)
+res23 <- zero_crossings_maxima(h23\$x,h23\$ratio)
+res24 <- zero_crossings_maxima(h24\$x,h24\$ratio)
+res25 <- zero_crossings_maxima(h25\$x,h25\$ratio)
+
+
+#plotting (including maxima and adjusting to max ratio)
+# all on one plot
+pdf("allchrs.pdf",height=10, width=20)
+library(Hmisc)
+ymax <- 1.05*max( c(h1\$ratio,h2\$ratio,h3\$ratio,h4\$ratio,h5\$ratio,h6\$ratio,h7\$ratio,h8\$ratio,h9\$ratio,h10\$ratio,h11\$ratio,h12\$ratio,h13\$ratio,h14\$ratio,h15\$ratio,h16\$ratio,h17\$ratio,h18\$ratio,h19\$ratio,h20\$ratio,h21\$ratio,h22\$ratio,h23\$ratio,h24\$ratio,h25\$ratio))
+par(mfrow=c(5,5))
+plot(h1\$x, h1\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 1")
+points( res1\$maxima\$x[1], res1\$maxima\$y[1], col="red", pch=19 )
+abline(v=res1\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res1\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h2\$x, h2\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 2")
+points( res2\$maxima\$x[1], res2\$maxima\$y[1], col="red", pch=19 )
+abline(v=res2\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res2\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h3\$x, h3\$ratio,type="l", lwd=3, col="red", xlim=c(0,64000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 3")
+points( res3\$maxima\$x[1], res3\$maxima\$y[1], col="red", pch=19 )
+abline(v=res3\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res3\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h4\$x, h4\$ratio,type="l", lwd=3, col="red",, xlim=c(0,64000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 4")
+points( res4\$maxima\$x[1], res4\$maxima\$y[1], col="red", pch=19 )
+abline(v=res4\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res4\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h5\$x, h5\$ratio,type="l", lwd=3, col="red",, xlim=c(0,76000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 5")
+points( res5\$maxima\$x[1], res5\$maxima\$y[1], col="red", pch=19 )
+abline(v=res5\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res5\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h6\$x, h6\$ratio,type="l", lwd=3, col="red",, xlim=c(0,60000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 6")
+points( res6\$maxima\$x[1], res6\$maxima\$y[1], col="red", pch=19 )
+abline(v=res6\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res6\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h7\$x, h7\$ratio,type="l", lwd=3, col="red",, xlim=c(0,78000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 7")
+points( res7\$maxima\$x[1], res7\$maxima\$y[1], col="red", pch=19 )
+abline(v=res7\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res7\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h8\$x, h8\$ratio,type="l", lwd=3, col="red",, xlim=c(0,57000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 8")
+points( res8\$maxima\$x[1], res8\$maxima\$y[1], col="red", pch=19 )
+abline(v=res8\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res8\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h9\$x, h9\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 9")
+points( res9\$maxima\$x[1], res9\$maxima\$y[1], col="red", pch=19 )
+abline(v=res9\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res9\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h10\$x, h10\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 10")
+points( res10\$maxima\$x[1], res10\$maxima\$y[1], col="red", pch=19 )
+abline(v=res10\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res10\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h11\$x, h11\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 11")
+points( res11\$maxima\$x[1], res11\$maxima\$y[1], col="red", pch=19 )
+abline(v=res11\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res11\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h12\$x, h12\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 12")
+points( res12\$maxima\$x[1], res12\$maxima\$y[1], col="red", pch=19 )
+abline(v=res12\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res12\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h13\$x, h13\$ratio,type="l", lwd=3, col="red",, xlim=c(0,55000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 13")
+points( res13\$maxima\$x[1], res13\$maxima\$y[1], col="red", pch=19 )
+abline(v=res13\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res13\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h14\$x, h14\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 14")
+points( res14\$maxima\$x[1], res14\$maxima\$y[1], col="red", pch=19 )
+abline(v=res14\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res14\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h15\$x, h15\$ratio,type="l", lwd=3, col="red",, xlim=c(0,48000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 15")
+points( res15\$maxima\$x[1], res15\$maxima\$y[1], col="red", pch=19 )
+abline(v=res15\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res15\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h16\$x, h16\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 16")
+points( res16\$maxima\$x[1], res16\$maxima\$y[1], col="red", pch=19 )
+abline(v=res16\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res16\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h17\$x, h17\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000),, ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 17")
+points( res17\$maxima\$x[1], res17\$maxima\$y[1], col="red", pch=19 )
+abline(v=res17\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res17\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h18\$x, h18\$ratio,type="l", lwd=3, col="red", xlim=c(0,50000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 18")
+points( res18\$maxima\$x[1], res18\$maxima\$y[1], col="red", pch=19 )
+abline(v=res18\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res18\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h19\$x, h19\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", , main="Chromosome 19")
+points( res19\$maxima\$x[1], res19\$maxima\$y[1], col="red", pch=19 )
+abline(v=res19\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res19\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h20\$x, h20\$ratio,type="l", lwd=3, col="red",, xlim=c(0,56000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 20")
+points( res20\$maxima\$x[1], res20\$maxima\$y[1], col="red", pch=19 )
+abline(v=res20\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res20\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h21\$x, h21\$ratio,type="l", lwd=3, col="red",, xlim=c(0,45000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 21")
+points( res21\$maxima\$x[1], res21\$maxima\$y[1], col="red", pch=19 )
+abline(v=res21\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res21\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h22\$x, h22\$ratio,type="l", lwd=3, col="red",, xlim=c(0,43000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 22")
+points( res22\$maxima\$x[1], res22\$maxima\$y[1], col="red", pch=19 )
+abline(v=res22\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res22\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h23\$x, h23\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 23")
+points( res23\$maxima\$x[1], res23\$maxima\$y[1], col="red", pch=19 )
+abline(v=res23\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res23\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h24\$x, h24\$ratio,type="l", lwd=3, col="red",, xlim=c(0,44000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 24")
+points( res24\$maxima\$x[1], res24\$maxima\$y[1], col="red", pch=19 )
+abline(v=res24\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res24\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+plot(h25\$x, h25\$ratio,type="l", lwd=3, col="red",, xlim=c(0,39000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 25")
+points( res25\$maxima\$x[1], res25\$maxima\$y[1], col="red", pch=19 )
+abline(v=res25\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", paste(format(round(res25\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(10)
+dev.off()
+
+# seperate plots
+pdf("chr1.pdf",height=5, width=10)
+plot(h1\$x, h1\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 1")
+points( res1\$maxima\$x[1], res1\$maxima\$y[1], col="red", pch=19 )
+abline(v=res1\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res1\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr2.pdf",height=5, width=10)
+plot(h2\$x, h2\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 2")
+points( res2\$maxima\$x[1], res2\$maxima\$y[1], col="red", pch=19 )
+abline(v=res2\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res2\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr3.pdf",height=5, width=10)
+plot(h3\$x, h3\$ratio,type="l", lwd=3, col="red", xlim=c(0,64000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 3")
+points( res3\$maxima\$x[1], res3\$maxima\$y[1], col="red", pch=19 )
+abline(v=res3\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res3\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr4.pdf",height=5, width=10)
+plot(h4\$x, h4\$ratio,type="l", lwd=3, col="red",, xlim=c(0,64000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 4")
+points( res4\$maxima\$x[1], res4\$maxima\$y[1], col="red", pch=19 )
+abline(v=res4\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res4\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr5.pdf",height=5, width=10)
+plot(h5\$x, h5\$ratio,type="l", lwd=3, col="red",, xlim=c(0,76000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 5")
+points( res5\$maxima\$x[1], res5\$maxima\$y[1], col="red", pch=19 )
+abline(v=res5\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res5\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr6.pdf",height=5, width=10)
+plot(h6\$x, h6\$ratio,type="l", lwd=3, col="red",, xlim=c(0,60000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 6")
+points( res6\$maxima\$x[1], res6\$maxima\$y[1], col="red", pch=19 )
+abline(v=res6\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res6\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr7.pdf",height=5, width=10)
+plot(h7\$x, h7\$ratio,type="l", lwd=3, col="red",, xlim=c(0,78000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 7")
+points( res7\$maxima\$x[1], res7\$maxima\$y[1], col="red", pch=19 )
+abline(v=res7\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res7\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr8.pdf",height=5, width=10)
+plot(h8\$x, h8\$ratio,type="l", lwd=3, col="red",, xlim=c(0,57000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 8")
+points( res8\$maxima\$x[1], res8\$maxima\$y[1], col="red", pch=19 )
+abline(v=res8\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res8\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr9.pdf",height=5, width=10)
+plot(h9\$x, h9\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 9")
+points( res9\$maxima\$x[1], res9\$maxima\$y[1], col="red", pch=19 )
+abline(v=res9\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res9\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr10.pdf",height=5, width=10)
+plot(h10\$x, h10\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 10")
+points( res10\$maxima\$x[1], res10\$maxima\$y[1], col="red", pch=19 )
+abline(v=res10\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res10\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr11.pdf",height=5, width=10)
+plot(h11\$x, h11\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 11")
+points( res11\$maxima\$x[1], res11\$maxima\$y[1], col="red", pch=19 )
+abline(v=res11\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res11\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr12.pdf",height=5, width=10)
+plot(h12\$x, h12\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 12")
+points( res12\$maxima\$x[1], res12\$maxima\$y[1], col="red", pch=19 )
+abline(v=res12\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res12\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr13.pdf",height=5, width=10)
+plot(h13\$x, h13\$ratio,type="l", lwd=3, col="red",, xlim=c(0,55000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 13")
+points( res13\$maxima\$x[1], res13\$maxima\$y[1], col="red", pch=19 )
+abline(v=res13\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res13\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr14.pdf",height=5, width=10)
+plot(h14\$x, h14\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 14")
+points( res14\$maxima\$x[1], res14\$maxima\$y[1], col="red", pch=19 )
+abline(v=res14\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res14\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr15.pdf",height=5, width=10)
+plot(h15\$x, h15\$ratio,type="l", lwd=3, col="red",, xlim=c(0,48000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 15")
+points( res15\$maxima\$x[1], res15\$maxima\$y[1], col="red", pch=19 )
+abline(v=res15\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res15\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr16.pdf",height=5, width=10)
+plot(h16\$x, h16\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 16")
+points( res16\$maxima\$x[1], res16\$maxima\$y[1], col="red", pch=19 )
+abline(v=res16\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res16\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr17.pdf",height=5, width=10)
+plot(h17\$x, h17\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000),, ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 17")
+points( res17\$maxima\$x[1], res17\$maxima\$y[1], col="red", pch=19 )
+abline(v=res17\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res17\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr18.pdf",height=5, width=10)
+plot(h18\$x, h18\$ratio,type="l", lwd=3, col="red", xlim=c(0,50000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 18")
+points( res18\$maxima\$x[1], res18\$maxima\$y[1], col="red", pch=19 )
+abline(v=res18\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res18\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr19.pdf",height=5, width=10)
+plot(h19\$x, h19\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", , main="Chromosome 19")
+points( res19\$maxima\$x[1], res19\$maxima\$y[1], col="red", pch=19 )
+abline(v=res19\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res19\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr20.pdf",height=5, width=10)
+plot(h20\$x, h20\$ratio,type="l", lwd=3, col="red",, xlim=c(0,56000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 20")
+points( res20\$maxima\$x[1], res20\$maxima\$y[1], col="red", pch=19 )
+abline(v=res20\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res20\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr21.pdf",height=5, width=10)
+plot(h21\$x, h21\$ratio,type="l", lwd=3, col="red",, xlim=c(0,45000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 21")
+points( res21\$maxima\$x[1], res21\$maxima\$y[1], col="red", pch=19 )
+abline(v=res21\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res21\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr22.pdf",height=5, width=10)
+plot(h22\$x, h22\$ratio,type="l", lwd=3, col="red",, xlim=c(0,43000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 22")
+points( res22\$maxima\$x[1], res22\$maxima\$y[1], col="red", pch=19 )
+abline(v=res22\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res22\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr23.pdf",height=5, width=10)
+plot(h23\$x, h23\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 23")
+points( res23\$maxima\$x[1], res23\$maxima\$y[1], col="red", pch=19 )
+abline(v=res23\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res23\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr24.pdf",height=5, width=10)
+plot(h24\$x, h24\$ratio,type="l", lwd=3, col="red",, xlim=c(0,44000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 24")
+points( res24\$maxima\$x[1], res24\$maxima\$y[1], col="red", pch=19 )
+abline(v=res24\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", legend=paste(format(round(res24\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(nx=10)
+dev.off()
+pdf("chr25.pdf",height=5, width=10)
+plot(h25\$x, h25\$ratio,type="l", lwd=3, col="red",, xlim=c(0,39000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 25")
+points( res25\$maxima\$x[1], res25\$maxima\$y[1], col="red", pch=19 )
+abline(v=res25\$maxima\$x[1], lwd=1)
+legend("topleft", bty="n", paste(format(round(res25\$maxima\$x[1] ,0),nsmall=0) ) )
+minor.tick(10)
+dev.off()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**What it Does**
+
+**Script**
+Pressing execute will run the following code over your input file and generate some outputs in your history::
+
+
+ #VDM-mapping script
+ #input file: Q200 variants no wildtype
+ ourargs = commandArgs(TRUE)
+ inf = ourargs[1]
+ outf = ourargs[2]
+ Q200var = read.table(inf)
+
+
+ het <- Q200var[grep("AF=0.5",Q200var\$V8), ]
+ hom <- Q200var[grep("AF=1.0",Q200var\$V8), ]
+
+ get_hom1 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr1")
+ het <- subset(x2,x2\$V1=="chr1")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom2 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr2")
+ het <- subset(x2,x2\$V1=="chr2")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom3 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr3")
+ het <- subset(x2,x2\$V1=="chr3")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom4 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr4")
+ het <- subset(x2,x2\$V1=="chr4")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom5 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr5")
+ het <- subset(x2,x2\$V1=="chr5")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom6 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr6")
+ het <- subset(x2,x2\$V1=="chr6")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom7 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr7")
+ het <- subset(x2,x2\$V1=="chr7")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom8 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr8")
+ het <- subset(x2,x2\$V1=="chr8")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom9 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr9")
+ het <- subset(x2,x2\$V1=="chr9")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom10 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr10")
+ het <- subset(x2,x2\$V1=="chr10")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom11 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr11")
+ het <- subset(x2,x2\$V1=="chr11")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom12 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr12")
+ het <- subset(x2,x2\$V1=="chr12")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom13 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr13")
+ het <- subset(x2,x2\$V1=="chr13")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom14 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr14")
+ het <- subset(x2,x2\$V1=="chr14")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom15 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr15")
+ het <- subset(x2,x2\$V1=="chr15")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom16 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr16")
+ het <- subset(x2,x2\$V1=="chr16")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom17 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr17")
+ het <- subset(x2,x2\$V1=="chr17")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom18 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr18")
+ het <- subset(x2,x2\$V1=="chr18")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom19 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr19")
+ het <- subset(x2,x2\$V1=="chr19")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom20 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr20")
+ het <- subset(x2,x2\$V1=="chr20")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom21 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr21")
+ het <- subset(x2,x2\$V1=="chr21")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom22 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr22")
+ het <- subset(x2,x2\$V1=="chr22")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom23 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr23")
+ het <- subset(x2,x2\$V1=="chr23")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom24 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr24")
+ het <- subset(x2,x2\$V1=="chr24")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ get_hom25 <- function(x1,x2){
+ hom <- subset(x1,x1\$V1=="chr25")
+ het <- subset(x2,x2\$V1=="chr25")
+ d1 <- density(hom\$V2)
+ d2 <- density(het\$V2)
+ return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+
+ h1 <- get_hom1(hom,het)
+ h2 <- get_hom2(hom,het)
+ h3 <- get_hom3(hom,het)
+ h4 <- get_hom4(hom,het)
+ h5 <- get_hom5(hom,het)
+ h6 <- get_hom6(hom,het)
+ h7 <- get_hom7(hom,het)
+ h8 <- get_hom8(hom,het)
+ h9 <- get_hom9(hom,het)
+ h10 <- get_hom10(hom,het)
+ h11 <- get_hom11(hom,het)
+ h12 <- get_hom12(hom,het)
+ h13 <- get_hom13(hom,het)
+ h14 <- get_hom14(hom,het)
+ h15 <- get_hom15(hom,het)
+ h16 <- get_hom16(hom,het)
+ h17 <- get_hom17(hom,het)
+ h18 <- get_hom18(hom,het)
+ h19 <- get_hom19(hom,het)
+ h20 <- get_hom20(hom,het)
+ h21 <- get_hom21(hom,het)
+ h22 <- get_hom22(hom,het)
+ h23 <- get_hom23(hom,het)
+ h24 <- get_hom24(hom,het)
+ h25 <- get_hom25(hom,het)
+
+
+
+ #Extracting maxima (VDM-mapping)
+
+ zero_crossings_maxima <- function(x,y){
+
+ dy <- diff(y)
+ xn <- x[2:length(x)]
+
+ n <- length(dy)
+ t1 <- dy[1:n-1];
+ t2 <- dy[2:n];
+ tt <- t1*t2
+ ind <- which(tt < 0)
+
+ # print out the maxima
+ d2y <- diff(dy)
+ vind <- c()
+ vy <- c()
+ vxn <- c()
+
+ xnn <- x[3:length(x)]
+ for(i in 1:length(ind)){
+ ### cat( ind[i], d2y[ind[i]], '\n')
+ if( d2y[ind[i]] < 0 ){
+ ### cat(ind[i], y[ind[i]], xn[ind[i]], "\n" )
+ vind <- c(vind, ind[i])
+ vy <- c(vy, y[ind[i]])
+ vxn <- c(vxn, xn[ind[i]])
+ }
+ }
+
+ maxima <- data.frame(ind=vind, y=vy, x=vxn)
+ maxima <- maxima[order(maxima\$y),]
+
+ #print( maxima )
+
+ results <- list()
+ results\$dy <- dy
+ results\$xn <- xn
+ results\$maxima <- maxima[ order(maxima\$y, decreasing=T), ]
+
+ return(results)
+ }
+
+ res1 <- zero_crossings_maxima(h1\$x,h1\$ratio)
+ res2 <- zero_crossings_maxima(h2\$x,h2\$ratio)
+ res3 <- zero_crossings_maxima(h3\$x,h3\$ratio)
+ res4 <- zero_crossings_maxima(h4\$x,h4\$ratio)
+ res5 <- zero_crossings_maxima(h5\$x,h5\$ratio)
+ res6 <- zero_crossings_maxima(h6\$x,h6\$ratio)
+ res7 <- zero_crossings_maxima(h7\$x,h7\$ratio)
+ res8 <- zero_crossings_maxima(h8\$x,h8\$ratio)
+ res9 <- zero_crossings_maxima(h9\$x,h9\$ratio)
+ res10 <- zero_crossings_maxima(h10\$x,h10\$ratio)
+ res11 <- zero_crossings_maxima(h11\$x,h11\$ratio)
+ res12 <- zero_crossings_maxima(h12\$x,h12\$ratio)
+ res13 <- zero_crossings_maxima(h13\$x,h13\$ratio)
+ res14 <- zero_crossings_maxima(h14\$x,h14\$ratio)
+ res15 <- zero_crossings_maxima(h15\$x,h15\$ratio)
+ res16 <- zero_crossings_maxima(h16\$x,h16\$ratio)
+ res17 <- zero_crossings_maxima(h17\$x,h17\$ratio)
+ res18 <- zero_crossings_maxima(h18\$x,h18\$ratio)
+ res19 <- zero_crossings_maxima(h19\$x,h19\$ratio)
+ res20 <- zero_crossings_maxima(h20\$x,h20\$ratio)
+ res21 <- zero_crossings_maxima(h21\$x,h21\$ratio)
+ res22 <- zero_crossings_maxima(h22\$x,h22\$ratio)
+ res23 <- zero_crossings_maxima(h23\$x,h23\$ratio)
+ res24 <- zero_crossings_maxima(h24\$x,h24\$ratio)
+ res25 <- zero_crossings_maxima(h25\$x,h25\$ratio)
+
+
+ #plotting (including maxima and adjusting to max ratio)
+ # all on one plot
+ pdf("allchrs.pdf",height=10, width=20)
+ library(Hmisc)
+ ymax <- 1.05*max( c(h1\$ratio,h2\$ratio,h3\$ratio,h4\$ratio,h5\$ratio,h6\$ratio,h7\$ratio,h8\$ratio,h9\$ratio,h10\$ratio,h11\$ratio,h12\$ratio,h13\$ratio,h14\$ratio,h15\$ratio,h16\$ratio,h17\$ratio,h18\$ratio,h19\$ratio,h20\$ratio,h21\$ratio,h22\$ratio,h23\$ratio,h24\$ratio,h25\$ratio))
+ par(mfrow=c(5,5))
+ plot(h1\$x, h1\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 1")
+ points( res1\$maxima\$x[1], res1\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res1\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res1\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h2\$x, h2\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 2")
+ points( res2\$maxima\$x[1], res2\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res2\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res2\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h3\$x, h3\$ratio,type="l", lwd=3, col="red", xlim=c(0,64000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 3")
+ points( res3\$maxima\$x[1], res3\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res3\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res3\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h4\$x, h4\$ratio,type="l", lwd=3, col="red",, xlim=c(0,64000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 4")
+ points( res4\$maxima\$x[1], res4\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res4\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res4\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h5\$x, h5\$ratio,type="l", lwd=3, col="red",, xlim=c(0,76000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 5")
+ points( res5\$maxima\$x[1], res5\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res5\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res5\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h6\$x, h6\$ratio,type="l", lwd=3, col="red",, xlim=c(0,60000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 6")
+ points( res6\$maxima\$x[1], res6\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res6\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res6\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h7\$x, h7\$ratio,type="l", lwd=3, col="red",, xlim=c(0,78000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 7")
+ points( res7\$maxima\$x[1], res7\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res7\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res7\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h8\$x, h8\$ratio,type="l", lwd=3, col="red",, xlim=c(0,57000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 8")
+ points( res8\$maxima\$x[1], res8\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res8\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res8\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h9\$x, h9\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 9")
+ points( res9\$maxima\$x[1], res9\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res9\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res9\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h10\$x, h10\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 10")
+ points( res10\$maxima\$x[1], res10\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res10\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res10\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h11\$x, h11\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 11")
+ points( res11\$maxima\$x[1], res11\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res11\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res11\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h12\$x, h12\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 12")
+ points( res12\$maxima\$x[1], res12\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res12\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res12\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h13\$x, h13\$ratio,type="l", lwd=3, col="red",, xlim=c(0,55000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 13")
+ points( res13\$maxima\$x[1], res13\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res13\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res13\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h14\$x, h14\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 14")
+ points( res14\$maxima\$x[1], res14\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res14\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res14\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h15\$x, h15\$ratio,type="l", lwd=3, col="red",, xlim=c(0,48000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 15")
+ points( res15\$maxima\$x[1], res15\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res15\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res15\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h16\$x, h16\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 16")
+ points( res16\$maxima\$x[1], res16\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res16\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res16\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h17\$x, h17\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000),, ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 17")
+ points( res17\$maxima\$x[1], res17\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res17\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res17\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h18\$x, h18\$ratio,type="l", lwd=3, col="red", xlim=c(0,50000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 18")
+ points( res18\$maxima\$x[1], res18\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res18\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res18\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h19\$x, h19\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", , main="Chromosome 19")
+ points( res19\$maxima\$x[1], res19\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res19\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res19\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h20\$x, h20\$ratio,type="l", lwd=3, col="red",, xlim=c(0,56000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 20")
+ points( res20\$maxima\$x[1], res20\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res20\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res20\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h21\$x, h21\$ratio,type="l", lwd=3, col="red",, xlim=c(0,45000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 21")
+ points( res21\$maxima\$x[1], res21\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res21\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res21\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h22\$x, h22\$ratio,type="l", lwd=3, col="red",, xlim=c(0,43000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 22")
+ points( res22\$maxima\$x[1], res22\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res22\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res22\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h23\$x, h23\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 23")
+ points( res23\$maxima\$x[1], res23\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res23\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res23\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h24\$x, h24\$ratio,type="l", lwd=3, col="red",, xlim=c(0,44000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 24")
+ points( res24\$maxima\$x[1], res24\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res24\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res24\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ plot(h25\$x, h25\$ratio,type="l", lwd=3, col="red",, xlim=c(0,39000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 25")
+ points( res25\$maxima\$x[1], res25\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res25\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", paste(format(round(res25\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(10)
+ dev.off()
+
+ # seperate plots
+ pdf("chr1.pdf",height=5, width=10)
+ plot(h1\$x, h1\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 1")
+ points( res1\$maxima\$x[1], res1\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res1\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res1\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr2.pdf",height=5, width=10)
+ plot(h2\$x, h2\$ratio,type="l", lwd=3, col="red", xlim=c(0,61000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 2")
+ points( res2\$maxima\$x[1], res2\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res2\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res2\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr3.pdf",height=5, width=10)
+ plot(h3\$x, h3\$ratio,type="l", lwd=3, col="red", xlim=c(0,64000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 3")
+ points( res3\$maxima\$x[1], res3\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res3\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res3\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr4.pdf",height=5, width=10)
+ plot(h4\$x, h4\$ratio,type="l", lwd=3, col="red",, xlim=c(0,64000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 4")
+ points( res4\$maxima\$x[1], res4\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res4\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res4\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr5.pdf",height=5, width=10)
+ plot(h5\$x, h5\$ratio,type="l", lwd=3, col="red",, xlim=c(0,76000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 5")
+ points( res5\$maxima\$x[1], res5\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res5\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res5\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr6.pdf",height=5, width=10)
+ plot(h6\$x, h6\$ratio,type="l", lwd=3, col="red",, xlim=c(0,60000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 6")
+ points( res6\$maxima\$x[1], res6\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res6\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res6\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr7.pdf",height=5, width=10)
+ plot(h7\$x, h7\$ratio,type="l", lwd=3, col="red",, xlim=c(0,78000000),ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 7")
+ points( res7\$maxima\$x[1], res7\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res7\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res7\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr8.pdf",height=5, width=10)
+ plot(h8\$x, h8\$ratio,type="l", lwd=3, col="red",, xlim=c(0,57000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 8")
+ points( res8\$maxima\$x[1], res8\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res8\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res8\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr9.pdf",height=5, width=10)
+ plot(h9\$x, h9\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio" , main="Chromosome 9")
+ points( res9\$maxima\$x[1], res9\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res9\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res9\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr10.pdf",height=5, width=10)
+ plot(h10\$x, h10\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 10")
+ points( res10\$maxima\$x[1], res10\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res10\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res10\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr11.pdf",height=5, width=10)
+ plot(h11\$x, h11\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 11")
+ points( res11\$maxima\$x[1], res11\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res11\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res11\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr12.pdf",height=5, width=10)
+ plot(h12\$x, h12\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 12")
+ points( res12\$maxima\$x[1], res12\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res12\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res12\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr13.pdf",height=5, width=10)
+ plot(h13\$x, h13\$ratio,type="l", lwd=3, col="red",, xlim=c(0,55000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 13")
+ points( res13\$maxima\$x[1], res13\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res13\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res13\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr14.pdf",height=5, width=10)
+ plot(h14\$x, h14\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 14")
+ points( res14\$maxima\$x[1], res14\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res14\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res14\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr15.pdf",height=5, width=10)
+ plot(h15\$x, h15\$ratio,type="l", lwd=3, col="red",, xlim=c(0,48000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 15")
+ points( res15\$maxima\$x[1], res15\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res15\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res15\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr16.pdf",height=5, width=10)
+ plot(h16\$x, h16\$ratio,type="l", lwd=3, col="red",, xlim=c(0,59000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 16")
+ points( res16\$maxima\$x[1], res16\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res16\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res16\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr17.pdf",height=5, width=10)
+ plot(h17\$x, h17\$ratio,type="l", lwd=3, col="red",, xlim=c(0,54000000),, ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 17")
+ points( res17\$maxima\$x[1], res17\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res17\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res17\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr18.pdf",height=5, width=10)
+ plot(h18\$x, h18\$ratio,type="l", lwd=3, col="red", xlim=c(0,50000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 18")
+ points( res18\$maxima\$x[1], res18\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res18\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res18\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr19.pdf",height=5, width=10)
+ plot(h19\$x, h19\$ratio,type="l", lwd=3, col="red",, xlim=c(0,51000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", , main="Chromosome 19")
+ points( res19\$maxima\$x[1], res19\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res19\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res19\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr20.pdf",height=5, width=10)
+ plot(h20\$x, h20\$ratio,type="l", lwd=3, col="red",, xlim=c(0,56000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 20")
+ points( res20\$maxima\$x[1], res20\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res20\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res20\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr21.pdf",height=5, width=10)
+ plot(h21\$x, h21\$ratio,type="l", lwd=3, col="red",, xlim=c(0,45000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 21")
+ points( res21\$maxima\$x[1], res21\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res21\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res21\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr22.pdf",height=5, width=10)
+ plot(h22\$x, h22\$ratio,type="l", lwd=3, col="red",, xlim=c(0,43000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 22")
+ points( res22\$maxima\$x[1], res22\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res22\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res22\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr23.pdf",height=5, width=10)
+ plot(h23\$x, h23\$ratio,type="l", lwd=3, col="red",, xlim=c(0,47000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 23")
+ points( res23\$maxima\$x[1], res23\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res23\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res23\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr24.pdf",height=5, width=10)
+ plot(h24\$x, h24\$ratio,type="l", lwd=3, col="red",, xlim=c(0,44000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 24")
+ points( res24\$maxima\$x[1], res24\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res24\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", legend=paste(format(round(res24\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(nx=10)
+ dev.off()
+ pdf("chr25.pdf",height=5, width=10)
+ plot(h25\$x, h25\$ratio,type="l", lwd=3, col="red",, xlim=c(0,39000000), ylim=c(0,ymax), xlab="Location (Mb)", ylab="Adjusted homozygosity ratio", main="Chromosome 25")
+ points( res25\$maxima\$x[1], res25\$maxima\$y[1], col="red", pch=19 )
+ abline(v=res25\$maxima\$x[1], lwd=1)
+ legend("topleft", bty="n", paste(format(round(res25\$maxima\$x[1] ,0),nsmall=0) ) )
+ minor.tick(10)
+ dev.off()
+
+**Attribution**
+This Galaxy tool was created by r.poole@ucl.ac.uk at 13/10/2014 12:55:18
+using the Galaxy Tool Factory.
+
+See https://bitbucket.org/fubar/galaxytoolfactory for details of that project
+Please cite: Creating re-usable tools from scripts: The Galaxy Tool Factory. Ross Lazarus; Antony Kaspi; Mark Ziemann; The Galaxy Team.
+Bioinformatics 2012; doi: 10.1093/bioinformatics/bts573
+
+
+
+
+
+ 10.1093/bioinformatics/bts573
+
+
diff -r c3956b7cb71f -r 810bc97ddeb2 Mydynamicscript/test-data/Mydynamicscript_test1_input.xls
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Mydynamicscript/test-data/Mydynamicscript_test1_input.xls Mon Oct 13 07:57:16 2014 -0400
@@ -0,0 +1,35467 @@
+##fileformat=VCFv4.1
+##FILTER=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##GATKCommandLine=
+##GATKCommandLine=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=
+##contig=