# HG changeset patch
# User richjp
# Date 1413204139 14400
# Node ID a08d6f405830bb7d03c0ebe4c3c2e7049d586240
# Parent 89411f98d1b781be98dac535eb15ce813f77b0a1
Uploaded
diff -r 89411f98d1b7 -r a08d6f405830 CloudMapzebrafishhomozygositymapping/CloudMapzebrafishhomozygositymapping.py
--- a/CloudMapzebrafishhomozygositymapping/CloudMapzebrafishhomozygositymapping.py Mon Oct 13 08:26:32 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,742 +0,0 @@
-# CloudMapzebrafishhomozygositymapping/CloudMapzebrafishhomozygositymapping.py - a self annotated version of rgToolFactory.py generated by running rgToolFactory.py
-# to make a new Galaxy tool called CloudMapzebrafishhomozygositymapping
-# User r.poole@ucl.ac.uk at 13/10/2014 13:24:10
-# 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 89411f98d1b7 -r a08d6f405830 CloudMapzebrafishhomozygositymapping/CloudMapzebrafishhomozygositymapping.xml
--- a/CloudMapzebrafishhomozygositymapping/CloudMapzebrafishhomozygositymapping.xml Mon Oct 13 08:26:32 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,1301 +0,0 @@
-
-Map a zebrafish mutant using bulk sergeant homozygosity mapping
-
- ghostscript
- graphicsmagick
-
-
-
- CloudMapzebrafishhomozygositymapping.py --script_path "$runMe" --interpreter "Rscript"
- --tool_name "CloudMapzebrafishhomozygositymapping" --input_tab "$input1" --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes"
-
-
-
-
-
-
-
-
-
-
-
-
-#homozygosity-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 (homozygosity-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**
-
-This tool is part of an update to the CloudMap pipeline (Minevich et al. 2012) to incorporate homozygosity mapping in zebrafish. This tool facilitates the mapping of zebrafish mutants from WGS or RNA-Seq data using bulk segregant homozygosity mapping. As input it takes a VCF list of variants (usually quality filtered for Q200 with WT SNPs/indels removed). This file is produced when running the CloudMap zebrafish pipeline. It plots an adjusted homozygosity ratio such that a high ratio indicates mapping region. It outputs an HTML file with combined and individual chromosome plots as PDFs. The adjusted homozygosity ratio is calculated as described in:
-
-IN SUBMISSION
-
-**Script**
-Pressing execute will run the following code over your input file and generate some outputs in your history::
-
-
- #homozygosity-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 (homozygosity-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 13:24:10
-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.1534/genetics.112.144204
- 10.1093/bioinformatics/bts573
-
-
diff -r 89411f98d1b7 -r a08d6f405830 CloudMapzebrafishhomozygositymapping/test-data/CloudMapzebrafishhomozygositymapping_test1_input.xls
--- a/CloudMapzebrafishhomozygositymapping/test-data/CloudMapzebrafishhomozygositymapping_test1_input.xls Mon Oct 13 08:26:32 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,35467 +0,0 @@
-##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=