Mercurial > repos > richjp > cloudmap_zebrafish_hom_mapping

changeset 8:0ede96d29bd5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author richjp
date Mon, 13 Oct 2014 08:45:03 -0400
parents a08d6f405830
children 0fee75e8a18e
files cloudmap_zf_hom_mapping.xml
diffstat 1 files changed, 1301 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cloudmap_zf_hom_mapping.xml	Mon Oct 13 08:45:03 2014 -0400
@@ -0,0 +1,1301 @@
+<tool id="CloudMap: zebrafish homozygosity mapping" name="CloudMap: zebrafish homozygosity mapping" version="0.01">
+<description>Map a zebrafish mutant using bulk sergeant homozygosity mapping</description>
+<requirements>
+      <requirement type="package" version="9.10">ghostscript</requirement>
+      <requirement type="package" version="1.3.18">graphicsmagick</requirement>
+  </requirements>
+<command interpreter="python">
+
+        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" 
+</command>
+<inputs>
+<param name="input1"  type="data" format="vcf" label="Select a suitable input file from your history"/> 
+<param name="job_name" type="text" label="Supply a name for the outputs to remind you what they contain" value="cloudmap_zf_hom_mapping"/> 
+
+</inputs>
+<outputs>
+ <data format="html" name="html_file" label="${job_name}.html"/>
+
+</outputs>
+<configfiles>
+<configfile name="runMe">
+#homozygosity-mapping script
+#input file: Q200 variants no wildtype
+ourargs = commandArgs(TRUE)
+inf = ourargs[1]
+outf = ourargs[2]
+Q200var = read.table(inf)
+
+
+het &lt;- Q200var[grep("AF=0.5",Q200var\$V8), ]
+hom &lt;- Q200var[grep("AF=1.0",Q200var\$V8), ]
+
+get_hom1 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr1")
+  het &lt;- subset(x2,x2\$V1=="chr1")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom2 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr2")
+  het &lt;- subset(x2,x2\$V1=="chr2")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom3 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr3")
+  het &lt;- subset(x2,x2\$V1=="chr3")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom4 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr4")
+  het &lt;- subset(x2,x2\$V1=="chr4")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom5 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr5")
+  het &lt;- subset(x2,x2\$V1=="chr5")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom6 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr6")
+  het &lt;- subset(x2,x2\$V1=="chr6")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom7 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr7")
+  het &lt;- subset(x2,x2\$V1=="chr7")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom8 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr8")
+  het &lt;- subset(x2,x2\$V1=="chr8")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom9 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr9")
+  het &lt;- subset(x2,x2\$V1=="chr9")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom10 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr10")
+  het &lt;- subset(x2,x2\$V1=="chr10")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom11 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr11")
+  het &lt;- subset(x2,x2\$V1=="chr11")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom12 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr12")
+  het &lt;- subset(x2,x2\$V1=="chr12")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom13 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr13")
+  het &lt;- subset(x2,x2\$V1=="chr13")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom14 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr14")
+  het &lt;- subset(x2,x2\$V1=="chr14")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom15 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr15")
+  het &lt;- subset(x2,x2\$V1=="chr15")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom16 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr16")
+  het &lt;- subset(x2,x2\$V1=="chr16")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom17 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr17")
+  het &lt;- subset(x2,x2\$V1=="chr17")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom18 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr18")
+  het &lt;- subset(x2,x2\$V1=="chr18")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom19 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr19")
+  het &lt;- subset(x2,x2\$V1=="chr19")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom20 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr20")
+  het &lt;- subset(x2,x2\$V1=="chr20")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom21 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr21")
+  het &lt;- subset(x2,x2\$V1=="chr21")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom22 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr22")
+  het &lt;- subset(x2,x2\$V1=="chr22")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom23 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr23")
+  het &lt;- subset(x2,x2\$V1=="chr23")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom24 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr24")
+  het &lt;- subset(x2,x2\$V1=="chr24")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+get_hom25 &lt;- function(x1,x2){
+  hom &lt;- subset(x1,x1\$V1=="chr25")
+  het &lt;- subset(x2,x2\$V1=="chr25")
+  d1 &lt;- density(hom\$V2)
+  d2 &lt;- density(het\$V2)
+  return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+}
+
+h1 &lt;- get_hom1(hom,het)
+h2 &lt;- get_hom2(hom,het)
+h3 &lt;- get_hom3(hom,het)
+h4 &lt;- get_hom4(hom,het)
+h5 &lt;- get_hom5(hom,het)
+h6 &lt;- get_hom6(hom,het)
+h7 &lt;- get_hom7(hom,het)
+h8 &lt;- get_hom8(hom,het)
+h9 &lt;- get_hom9(hom,het)
+h10 &lt;- get_hom10(hom,het)
+h11 &lt;- get_hom11(hom,het)
+h12 &lt;- get_hom12(hom,het)
+h13 &lt;- get_hom13(hom,het)
+h14 &lt;- get_hom14(hom,het)
+h15 &lt;- get_hom15(hom,het)
+h16 &lt;- get_hom16(hom,het)
+h17 &lt;- get_hom17(hom,het)
+h18 &lt;- get_hom18(hom,het)
+h19 &lt;- get_hom19(hom,het)
+h20 &lt;- get_hom20(hom,het)
+h21 &lt;- get_hom21(hom,het)
+h22 &lt;- get_hom22(hom,het)
+h23 &lt;- get_hom23(hom,het)
+h24 &lt;- get_hom24(hom,het)
+h25 &lt;- get_hom25(hom,het)
+
+
+
+#Extracting maxima (homozygosity-mapping)
+
+zero_crossings_maxima &lt;- function(x,y){
+
+  dy &lt;- diff(y)
+  xn &lt;- x[2:length(x)]
+
+  n &lt;- length(dy)
+  t1 &lt;- dy[1:n-1];
+  t2 &lt;- dy[2:n];
+  tt &lt;- t1*t2
+  ind &lt;- which(tt &lt; 0)
+
+  # print out the maxima
+  d2y &lt;- diff(dy)
+  vind &lt;- c()
+  vy &lt;- c()
+  vxn &lt;- c()
+
+  xnn &lt;- x[3:length(x)]
+  for(i in 1:length(ind)){
+    ### cat( ind[i], d2y[ind[i]], '\n')
+    if( d2y[ind[i]] &lt; 0 ){
+      ### cat(ind[i], y[ind[i]], xn[ind[i]], "\n" )
+      vind &lt;- c(vind, ind[i])
+      vy &lt;- c(vy, y[ind[i]])
+      vxn &lt;- c(vxn, xn[ind[i]])
+    }
+  }
+
+  maxima &lt;- data.frame(ind=vind, y=vy, x=vxn)
+  maxima &lt;- maxima[order(maxima\$y),]
+
+  #print( maxima )
+
+  results &lt;- list()
+  results\$dy &lt;- dy
+  results\$xn &lt;- xn
+  results\$maxima &lt;- maxima[ order(maxima\$y, decreasing=T), ]
+
+  return(results)
+}
+
+res1 &lt;- zero_crossings_maxima(h1\$x,h1\$ratio)
+res2 &lt;- zero_crossings_maxima(h2\$x,h2\$ratio)
+res3 &lt;- zero_crossings_maxima(h3\$x,h3\$ratio)
+res4 &lt;- zero_crossings_maxima(h4\$x,h4\$ratio)
+res5 &lt;- zero_crossings_maxima(h5\$x,h5\$ratio)
+res6 &lt;- zero_crossings_maxima(h6\$x,h6\$ratio)
+res7 &lt;- zero_crossings_maxima(h7\$x,h7\$ratio)
+res8 &lt;- zero_crossings_maxima(h8\$x,h8\$ratio)
+res9 &lt;- zero_crossings_maxima(h9\$x,h9\$ratio)
+res10 &lt;- zero_crossings_maxima(h10\$x,h10\$ratio)
+res11 &lt;- zero_crossings_maxima(h11\$x,h11\$ratio)
+res12 &lt;- zero_crossings_maxima(h12\$x,h12\$ratio)
+res13 &lt;- zero_crossings_maxima(h13\$x,h13\$ratio)
+res14 &lt;- zero_crossings_maxima(h14\$x,h14\$ratio)
+res15 &lt;- zero_crossings_maxima(h15\$x,h15\$ratio)
+res16 &lt;- zero_crossings_maxima(h16\$x,h16\$ratio)
+res17 &lt;- zero_crossings_maxima(h17\$x,h17\$ratio)
+res18 &lt;- zero_crossings_maxima(h18\$x,h18\$ratio)
+res19 &lt;- zero_crossings_maxima(h19\$x,h19\$ratio)
+res20 &lt;- zero_crossings_maxima(h20\$x,h20\$ratio)
+res21 &lt;- zero_crossings_maxima(h21\$x,h21\$ratio)
+res22 &lt;- zero_crossings_maxima(h22\$x,h22\$ratio)
+res23 &lt;- zero_crossings_maxima(h23\$x,h23\$ratio)
+res24 &lt;- zero_crossings_maxima(h24\$x,h24\$ratio)
+res25 &lt;- 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 &lt;- 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()
+</configfile>
+</configfiles>
+
+
+        <tests>
+        <test>
+        <param name="input1" value="cloudmap_zf_hom_mapping_test1_input.xls" ftype="vcf"/>
+        <param name="job_name" value="test1"/>
+        <param name="runMe" value="$runMe"/>
+        <output name="html_file" file="cloudmap_zf_hom_mapping_test1_output.html" ftype="html" lines_diff="5"/>
+        </test>
+        </tests>
+        
+
+<help>
+
+
+**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 &lt;- Q200var[grep("AF=0.5",Q200var\$V8), ]
+ hom &lt;- Q200var[grep("AF=1.0",Q200var\$V8), ]
+ 
+ get_hom1 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr1")
+   het &lt;- subset(x2,x2\$V1=="chr1")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom2 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr2")
+   het &lt;- subset(x2,x2\$V1=="chr2")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom3 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr3")
+   het &lt;- subset(x2,x2\$V1=="chr3")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom4 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr4")
+   het &lt;- subset(x2,x2\$V1=="chr4")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom5 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr5")
+   het &lt;- subset(x2,x2\$V1=="chr5")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom6 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr6")
+   het &lt;- subset(x2,x2\$V1=="chr6")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom7 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr7")
+   het &lt;- subset(x2,x2\$V1=="chr7")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom8 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr8")
+   het &lt;- subset(x2,x2\$V1=="chr8")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom9 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr9")
+   het &lt;- subset(x2,x2\$V1=="chr9")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom10 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr10")
+   het &lt;- subset(x2,x2\$V1=="chr10")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom11 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr11")
+   het &lt;- subset(x2,x2\$V1=="chr11")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom12 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr12")
+   het &lt;- subset(x2,x2\$V1=="chr12")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom13 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr13")
+   het &lt;- subset(x2,x2\$V1=="chr13")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom14 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr14")
+   het &lt;- subset(x2,x2\$V1=="chr14")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom15 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr15")
+   het &lt;- subset(x2,x2\$V1=="chr15")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom16 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr16")
+   het &lt;- subset(x2,x2\$V1=="chr16")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom17 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr17")
+   het &lt;- subset(x2,x2\$V1=="chr17")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom18 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr18")
+   het &lt;- subset(x2,x2\$V1=="chr18")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom19 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr19")
+   het &lt;- subset(x2,x2\$V1=="chr19")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom20 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr20")
+   het &lt;- subset(x2,x2\$V1=="chr20")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom21 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr21")
+   het &lt;- subset(x2,x2\$V1=="chr21")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom22 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr22")
+   het &lt;- subset(x2,x2\$V1=="chr22")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom23 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr23")
+   het &lt;- subset(x2,x2\$V1=="chr23")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom24 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr24")
+   het &lt;- subset(x2,x2\$V1=="chr24")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ get_hom25 &lt;- function(x1,x2){
+   hom &lt;- subset(x1,x1\$V1=="chr25")
+   het &lt;- subset(x2,x2\$V1=="chr25")
+   d1 &lt;- density(hom\$V2)
+   d2 &lt;- density(het\$V2)
+   return(list( x=d1\$x, ratio=(d1\$y/d2\$y)*(nrow(hom)/nrow(het))))
+ }
+ 
+ h1 &lt;- get_hom1(hom,het)
+ h2 &lt;- get_hom2(hom,het)
+ h3 &lt;- get_hom3(hom,het)
+ h4 &lt;- get_hom4(hom,het)
+ h5 &lt;- get_hom5(hom,het)
+ h6 &lt;- get_hom6(hom,het)
+ h7 &lt;- get_hom7(hom,het)
+ h8 &lt;- get_hom8(hom,het)
+ h9 &lt;- get_hom9(hom,het)
+ h10 &lt;- get_hom10(hom,het)
+ h11 &lt;- get_hom11(hom,het)
+ h12 &lt;- get_hom12(hom,het)
+ h13 &lt;- get_hom13(hom,het)
+ h14 &lt;- get_hom14(hom,het)
+ h15 &lt;- get_hom15(hom,het)
+ h16 &lt;- get_hom16(hom,het)
+ h17 &lt;- get_hom17(hom,het)
+ h18 &lt;- get_hom18(hom,het)
+ h19 &lt;- get_hom19(hom,het)
+ h20 &lt;- get_hom20(hom,het)
+ h21 &lt;- get_hom21(hom,het)
+ h22 &lt;- get_hom22(hom,het)
+ h23 &lt;- get_hom23(hom,het)
+ h24 &lt;- get_hom24(hom,het)
+ h25 &lt;- get_hom25(hom,het)
+ 
+ 
+ 
+ #Extracting maxima (homozygosity-mapping)
+ 
+ zero_crossings_maxima &lt;- function(x,y){
+ 
+   dy &lt;- diff(y)
+   xn &lt;- x[2:length(x)]
+ 
+   n &lt;- length(dy)
+   t1 &lt;- dy[1:n-1];
+   t2 &lt;- dy[2:n];
+   tt &lt;- t1*t2
+   ind &lt;- which(tt &lt; 0)
+ 
+   # print out the maxima
+   d2y &lt;- diff(dy)
+   vind &lt;- c()
+   vy &lt;- c()
+   vxn &lt;- c()
+ 
+   xnn &lt;- x[3:length(x)]
+   for(i in 1:length(ind)){
+     ### cat( ind[i], d2y[ind[i]], '\n')
+     if( d2y[ind[i]] &lt; 0 ){
+       ### cat(ind[i], y[ind[i]], xn[ind[i]], "\n" )
+       vind &lt;- c(vind, ind[i])
+       vy &lt;- c(vy, y[ind[i]])
+       vxn &lt;- c(vxn, xn[ind[i]])
+     }
+   }
+ 
+   maxima &lt;- data.frame(ind=vind, y=vy, x=vxn)
+   maxima &lt;- maxima[order(maxima\$y),]
+ 
+   #print( maxima )
+ 
+   results &lt;- list()
+   results\$dy &lt;- dy
+   results\$xn &lt;- xn
+   results\$maxima &lt;- maxima[ order(maxima\$y, decreasing=T), ]
+ 
+   return(results)
+ }
+ 
+ res1 &lt;- zero_crossings_maxima(h1\$x,h1\$ratio)
+ res2 &lt;- zero_crossings_maxima(h2\$x,h2\$ratio)
+ res3 &lt;- zero_crossings_maxima(h3\$x,h3\$ratio)
+ res4 &lt;- zero_crossings_maxima(h4\$x,h4\$ratio)
+ res5 &lt;- zero_crossings_maxima(h5\$x,h5\$ratio)
+ res6 &lt;- zero_crossings_maxima(h6\$x,h6\$ratio)
+ res7 &lt;- zero_crossings_maxima(h7\$x,h7\$ratio)
+ res8 &lt;- zero_crossings_maxima(h8\$x,h8\$ratio)
+ res9 &lt;- zero_crossings_maxima(h9\$x,h9\$ratio)
+ res10 &lt;- zero_crossings_maxima(h10\$x,h10\$ratio)
+ res11 &lt;- zero_crossings_maxima(h11\$x,h11\$ratio)
+ res12 &lt;- zero_crossings_maxima(h12\$x,h12\$ratio)
+ res13 &lt;- zero_crossings_maxima(h13\$x,h13\$ratio)
+ res14 &lt;- zero_crossings_maxima(h14\$x,h14\$ratio)
+ res15 &lt;- zero_crossings_maxima(h15\$x,h15\$ratio)
+ res16 &lt;- zero_crossings_maxima(h16\$x,h16\$ratio)
+ res17 &lt;- zero_crossings_maxima(h17\$x,h17\$ratio)
+ res18 &lt;- zero_crossings_maxima(h18\$x,h18\$ratio)
+ res19 &lt;- zero_crossings_maxima(h19\$x,h19\$ratio)
+ res20 &lt;- zero_crossings_maxima(h20\$x,h20\$ratio)
+ res21 &lt;- zero_crossings_maxima(h21\$x,h21\$ratio)
+ res22 &lt;- zero_crossings_maxima(h22\$x,h22\$ratio)
+ res23 &lt;- zero_crossings_maxima(h23\$x,h23\$ratio)
+ res24 &lt;- zero_crossings_maxima(h24\$x,h24\$ratio)
+ res25 &lt;- 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 &lt;- 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:41:11
+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
+
+
+</help>
+<citations>
+    <citation type="doi">10.1534/genetics.112.144204</citation>
+    <citation type="doi">10.1093/bioinformatics/bts573</citation>
+</citations>
+</tool>