# HG changeset patch
# User richjp
# Date 1413201221 14400
# Node ID c3956b7cb71f394b2a355c58ea34b89c7c9d261d
# Parent f6ea3e1cdfb96f0efe37ae99d176660a5484dd6c
Deleted selected files
diff -r f6ea3e1cdfb9 -r c3956b7cb71f cloudmap_zf_hom_mapping/cloudmap_zf_hom_mapping.py
--- a/cloudmap_zf_hom_mapping/cloudmap_zf_hom_mapping.py Fri Oct 10 12:36:10 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,742 +0,0 @@
-# cloudmap_zf_hom_mapping/cloudmap_zf_hom_mapping.py - a self annotated version of rgToolFactory.py generated by running rgToolFactory.py
-# to make a new Galaxy tool called cloudmap_zf_hom_mapping
-# User r.poole@ucl.ac.uk at 10/10/2014 17:34:44
-# 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 f6ea3e1cdfb9 -r c3956b7cb71f cloudmap_zf_hom_mapping/cloudmap_zf_hom_mapping.xml
--- a/cloudmap_zf_hom_mapping/cloudmap_zf_hom_mapping.xml Fri Oct 10 12:36:10 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,1298 +0,0 @@
-
-This tool plots adjusted homozygosity ratios for zebrafish bulk sergeant analysis
-
-
-
- cloudmap_zf_hom_mapping.py --script_path "$runMe" --interpreter "Rscript"
- --tool_name "cloudmap_zf_hom_mapping" --input_tab "$input1" --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes"
-
-
-
-
-
-
-
-
-
-
-
-
-ourargs = commandArgs(TRUE)
-#Homozygosity-mapping script
-#input file: Q200 variants no wildtype
-
-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**
-
-This tool takes as input a VCF list of SNPs from pooled bulk-segregant zebrafish WGS or RNASeq data. This VCF should be quality filtered (usually Q200) and have had wild type SNPs/indels subtracted. This file is generated when running the zebrafish WGS or RNA-seq CloudMap workflows. It plots adjusted homozygosity ratios as in ores to determine the likely position of a causative mutation. This tool was created by Richard J. Poole (UCL) using TooolFactory and is described in:
-
-IN PRESS
-
-**Script**
-Pressing execute will run the following code over your input file and generate some outputs in your history::
-
-
- ourargs = commandArgs(TRUE)
- #Homozygosity-mapping script
- #input file: Q200 variants no wildtype
-
- 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 10/10/2014 17:34:44
-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 f6ea3e1cdfb9 -r c3956b7cb71f cloudmap_zf_hom_mapping/test-data/cloudmap_zf_hom_mapping_test1_input.xls
--- a/cloudmap_zf_hom_mapping/test-data/cloudmap_zf_hom_mapping_test1_input.xls Fri Oct 10 12:36:10 2014 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,31522 +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=