damidseq_polii_gene_call: polii.gene.call annotate

annotate polii.gene.call @ 0:1b5bd3b955ed draft

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty

author	mvdbeek
date	Fri, 20 Apr 2018 10:27:20 -0400
parents
children

rev	line source
0 1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	1 #!/usr/bin/env Rscript
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	2
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	3 # polii.gene.call.r
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	4 # Copyright © 2014-15, Owen Marshall
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	5
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	6 # This program is free software; you can redistribute it and/or modify
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	7 # it under the terms of the GNU General Public License as published by
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	8 # the Free Software Foundation; either version 2 of the License, or (at
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	9 # your option) any later version.
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	10 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	11 # This program is distributed in the hope that it will be useful, but
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	12 # WITHOUT ANY WARRANTY; without even the implied warranty of
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	14 # General Public License for more details.
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	15 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	16 # You should have received a copy of the GNU General Public License
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	17 # along with this program; if not, write to the Free Software
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	18 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	19 # USA
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	20
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	21 ### FDR calcs ###
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	22 # Method based on original perl scripts by Tony Southall (TDS) as published in
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	23 # Southall et al. (2013). Dev Cell, 26(1), 101–12. doi:10.1016/j.devcel.2013.05.020
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	24 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	25 # Significant modifications to the original methodology include:
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	26 # * taking a linear regression of log data rather than trial-and-error curve fitting of non-log data
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	27 # * using a linear regression for the final intercept value rather than using the average intercept value for all conditions
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	28 # -- both of these should increase the accuracy of the final FDR value.
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	29
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	30 version <- "1.0.2"
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	31 cat(paste("polii.gene.call v",version,"\n", sep=""))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	32
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	33 library(tools)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	34
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	35 ### Read CLI options
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	36 input.args <- commandArgs(trailingOnly = TRUE)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	37
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	38 in.files <- vector()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	39 read.ops <- function (x) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	40 for (op in x) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	41 if (any(grepl("^--",op))) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	42 op <- gsub("^--","",op)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	43 y <- unlist(strsplit(op,"="))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	44
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	45 if (y[1] == "help") {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	46 cat(paste("Usage: Rscript polii.gene.call.r --genes.file=[some_genes_file.gff] [list of .gatc.bedgraph or .gatc.gff ratio files to process]\n\n", sep=""))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	47
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	48 cat("Options:\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	49 for (n in names(op.args)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	50 cat(paste(" ",n,"=",op.args[[n]],"\n",sep=""))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	51 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	52 cat("\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	53 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	54 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	55
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	56 if (!is.null(op.args[[ y[1] ]])) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	57 op.args[[ y[1] ]] <<- y[2]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	58 } else {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	59 cat("Error: Option",y[1],"not recognised ...\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	60 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	61 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	62 } else {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	63 in.files <<- c(in.files,op)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	64 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	65 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	66 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	67
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	68 write.ops <- function () {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	69 out.df <- data.frame()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	70 for (n in names(op.args)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	71 v <<- as.character(op.args[[n]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	72 df.line <- data.frame(
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	73 option=n,
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	74 value=v
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	75 )
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	76 out.df <- rbind(out.df, df.line)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	77 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	78 write.table(out.df,"input.args.single.txt",row.names=F)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	79 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	80
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	81 op.args <- list(
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	82 "genes.file" = "/mnt/data/Genomes/dmel_release/DmR6/DmR6.genes.gff",
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	83 "iter" = 50000,
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	84 "fdr" = 0.01
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	85 )
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	86
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	87 read.ops(input.args)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	88
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	89 if (length(in.files) == 0) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	90 cat("Usage: Rscript polii.gene.call.r [list of .gatc.gff ratio files to process]\n\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	91 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	92 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	93
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	94 write.ops()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	95
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	96 ### save random seed for future reproducibility
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	97 dump.random <- runif(1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	98 my.seed <- .Random.seed
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	99 write.table(my.seed,".randomseed")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	100
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	101 ### read genes file
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	102 cat("Reading genes data file ...\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	103 genes.file=op.args[["genes.file"]]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	104 genes <- read.table(genes.file, comment.char="#", sep="\t", quote="", fill=T)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	105 names(genes) <- c('chr','source','type','start','end','score','strand','c','details')
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	106
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	107 # only subset if there is a type termed "gene"
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	108 if (any(genes$type == 'gene')) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	109 genes <- subset(genes, type=='gene')
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	110 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	111
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	112 genes$name <- sapply(genes$details, FUN = function (x) {regmatches(x,gregexpr("(?<=Name=).*?(?=;)", x, perl=T))} )
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	113 genes <- genes[,c('chr','start','end','strand','name')]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	114 genes$chr <- gsub("^chr","",genes$chr,perl=T)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	115
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	116 if (nrow(genes) == 0) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	117 cat("Error: unable to extract gene information from genes file\n\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	118 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	119 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	120
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	121 ### functions
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	122 read.gff <- function (x,name="score") {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	123 fn.ext <- file_ext(x)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	124
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	125 if (grepl("gff",ignore.case=T,fn.ext)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	126 temp.data <- read.table(x,row.names=NULL)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	127 if (ncol(temp.data) > 5) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	128 # GFF
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	129 trim.data <- temp.data[,c(1,4,5,6)]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	130 } else {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	131 cat("Error: file does not appear to be in GFF format\n\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	132 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	133 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	134 } else if (grepl("bed",ignore.case=T,fn.ext)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	135 temp.data <- read.table(x,row.names=NULL,skip=1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	136 if (ncol(temp.data) == 4) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	137 # bedgraph
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	138 trim.data <- temp.data
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	139 } else {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	140 cat("Error: file does not appear to be in bedGraph format\n\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	141 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	142 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	143 } else {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	144 cat("Error: input file does not appear to be in bedGraph or GFF format ...\n\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	145 quit("no",1)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	146 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	147
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	148 names(trim.data) <- c("chr","start","end",name)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	149
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	150 trim.data$chr <- gsub("^chr","",trim.data$chr,perl=T)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	151
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	152 return(trim.data)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	153 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	154
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	155 gene.exp <- function (input.df, buffer=0, iter=50000, debug=F) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	156
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	157 avg.exp <- data.frame(input.df[1,c(4:(length(names(input.df))))])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	158 avg <- vector(length=(length(names(input.df)) - 4))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	159 avg.exp <- avg.exp[0,]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	160
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	161 ### FDR calcs ###
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	162 # Method based off perl scripts by Tony Southall (TDS) as published in
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	163 # Southall et al. (2013). Dev Cell, 26(1), 101–12. doi:10.1016/j.devcel.2013.05.020
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	164 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	165 # Significant modifications to the original methodology include:
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	166 # * taking a linear regression of log data rather than trial-and-error curve fitting of non-log data
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	167 # * using a linear regression for the final intercept value rather than using the average intercept value for all conditions
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	168 # -- both of these should increase the accuracy of the final FDR value.
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	169
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	170 input.len <- length(input.df[,1])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	171 frag.samp <- c(1,2,3,4,6,8,10,12,15)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	172 thres.samp <- c(0.1,0.2,0.3,0.4,0.5,0.65,0.8,1.0,1.5,2.0)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	173
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	174 rand <- list()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	175
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	176 for (thres in thres.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	177
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	178 cat(paste(" Calculating FDR for threshold",thres,"\n",sep=" "))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	179
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	180 # init vars
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	181 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	182 # List names are, e.g. frag 1, thres 0.2: rand[[thres.1.0.2]]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	183 rand[[paste("thres.",f,".",thres,sep="")]] <- 0;
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	184 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	185
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	186 for (i in 1:iter) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	187 if (i %% 200 == 0) {cat(paste(" iter",i,"\r"))}
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	188 # get random sample for different fragment lengths
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	189
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	190 rand.samp <- list()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	191
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	192 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	193 # Using the fourth column as we're only calculating FDR for one sample ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	194 rand.samp[[paste("rand.",f,sep="")]] <- mean(input.df[runif(f,1,input.len),4])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	195 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	196
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	197 # count number of times exp > thres
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	198 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	199 if (rand.samp[[paste("rand.",f,sep="")]] > thres) {rand[[paste("thres.",f,".",thres,sep="")]] <- rand[[paste("thres.",f,".",thres,sep="")]] + 1}
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	200 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	201 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	202 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	203
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	204 rand.fdr <- list()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	205
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	206 for (thres in thres.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	207 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	208 rand.fdr[[paste("thres.",f,".",thres,sep="")]] <- rand[[paste("thres.",f,".",thres,sep="")]]/iter
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	209 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	210 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	211
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	212 cat("Fitting curves ...\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	213
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	214 # curve fit: fdr vs thresholds
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	215 var.thres <- list()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	216
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	217 for (thres in thres.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	218 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	219 var.thres[[paste("frags.",f,sep="")]] <- append(var.thres[[paste("frags.",f,sep="")]], rand.fdr[[paste("thres.",f,".",thres,sep="")]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	220 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	221 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	222
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	223 inf.log.lm <- function (v) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	224 non.inf <- log(v) != -Inf
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	225 ret <- lm(log(v)[non.inf] ~ thres.samp[non.inf])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	226 return(ret)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	227 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	228
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	229 # The relationship is exponential, so we need log data for a linear regression
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	230 # (in R, linear regression is: y = lm$coefficients[[2]]x + lm$coefficients[[1]] ... )
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	231 var.lm <- list()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	232 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	233 var.lm[[paste("frags.",f,sep="")]] <- inf.log.lm(var.thres[[paste("frags.",f,sep="")]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	234 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	235
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	236 # ... and now we do a linear regression on the slopes and intercepts of our previous regressions
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	237 # (This is the clever bit, and it actually seems to work. The correlation of slope to fragment size is linear ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	238 # By doing this on the slope and intercept, we can now predict the FDR for any number of fragments with any expression.)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	239 slope <- vector()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	240 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	241 slope <- append(slope, var.lm[[paste("frags.",f,sep="")]]$coefficients[[2]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	242 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	243
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	244 # slope regression predicts the average slope
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	245 slope.lm <- lm(slope ~ frag.samp)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	246
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	247 # TDS used an average intercept value for the intercept, however ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	248 inter <- vector()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	249 for (f in frag.samp) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	250 inter <- append(inter, var.lm[[paste("frags.",f,sep="")]]$coefficients[[1]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	251 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	252
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	253 # ... there's actually quite a bit of variation of the intercept with real data,
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	254 # so we're going to do a lin regression on the intercept instead.
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	255 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	256 # (I'm not convinced it's a true linear relationship. But it's close to linear,
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	257 # and will certainly perform better than taking the mean intercept ...)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	258 #
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	259 # If you're interested, set the debug flag to TRUE and take a look at the plots generated below ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	260 inter.lm <- lm(inter ~ frag.samp)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	261
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	262 if (debug == T) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	263 # plots for debugging/checking
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	264 plot.debug <- function (y,x,l,name="Debug plot") {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	265 plot(y ~ x)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	266 abline(l)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	267
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	268 lsum <- summary(l)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	269 r2 <- lsum$r.squared
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	270
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	271 legend("topleft",legend=r2,bty='n')
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	272 title(name)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	273
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	274 dev.copy(png,paste(name,".png",sep=""));dev.off()
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	275 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	276
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	277 plot.debug(slope,frag.samp,slope.lm,name="correlation of slope")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	278 plot.debug(inter,frag.samp,inter.lm,name="correlation of intercepts")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	279 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	280
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	281 # ok, so putting that all together ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	282 fdr <- function (frags, expr) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	283 inter.test <- inter.lm$coefficients[[2]] * frags + inter.lm$coefficients[[1]]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	284 slope.test <- slope.lm$coefficients[[2]] * frags + slope.lm$coefficients[[1]]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	285
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	286 fdr.out <- exp(slope.test * expr + inter.test)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	287 return(fdr.out)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	288 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	289
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	290 ### Gene expression values ###
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	291 cat("Calculating gene values ...\n")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	292
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	293 count <- 0
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	294
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	295 # unroll chromosomes for speed:
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	296 for (chromo in unique(genes$chr)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	297 input.chr <- subset(input.df, chr==chromo)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	298 genes.chr <- subset(genes, chr==chromo)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	299 for (i in 1:length(genes.chr$name)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	300 # Roll through each gene
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	301
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	302 # Note: the original script calculated expression values for a table of proteins,
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	303 # here we just limit ourselves to the FDR for the first column past "chr", "start" and "end"
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	304 # so if you're thinking of looking at chromatin, for example, put PolII as column 4 in your table!
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	305
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	306 gene.start <- genes.chr[i,"start"] - buffer
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	307 gene.end <- genes.chr[i,"end"] + buffer
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	308
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	309 gene.start <- ifelse(gene.start < 1, 1, gene.start)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	310
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	311 # Create data frames for all gatc fragments covering current gene
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	312 exp <- data.frame(input.chr[ (input.chr$start <= gene.end)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	313 & (input.chr$end >= gene.start)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	314 ,] )
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	315
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	316 gatc.num <- length(exp[,1])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	317
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	318 # skip if no gatc fragments cover gene :(
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	319 if (gatc.num == 0) {next}
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	320
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	321 # trim to gene boundaries ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	322 exp$start[1] <- gene.start
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	323 exp$end[length(exp[,1])] <- gene.end
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	324
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	325 # gene length covered by gatc fragments
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	326 len <- sum(exp$end-exp$start)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	327
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	328 # calculate weighted score for each column (representing different proteins)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	329 for (j in 4:length(names(input.chr))) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	330 avg[j] <- (sum((exp$end-exp$start)*exp[j]))/len
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	331 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	332
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	333 # make data.frame of averages (to be appended to avg.exp)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	334 df <- cbind(avg[1])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	335 for (k in 2:length(avg)) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	336 df <- cbind(df,avg[k])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	337 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	338 df <- cbind(df,gatc.num)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	339
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	340 # only fdr for first column for now ...
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	341 gene.fdr <- fdr(gatc.num,avg[4])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	342 df <- cbind(df, gene.fdr)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	343
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	344 # append current gene to list
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	345 avg.exp <- rbind(avg.exp,data.frame(name=as.character(genes.chr[i,"name"]), df))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	346 count <- count+1
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	347 if (count %% 50 == 0) {cat(paste(count,"genes averaged ...\r"))}
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	348 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	349 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	350
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	351 avg.exp <- avg.exp[,c(1,5:(length(names(avg.exp))))]
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	352 names(avg.exp) <- c("name",names(input.df)[c(4:(length(names(input.df))))],"gatc.num","FDR")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	353
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	354 return(avg.exp)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	355 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	356
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	357 for (name in in.files) {
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	358 cat(paste("\nNow working on",name,"...\n"))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	359 bname <- basename(name)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	360 fname <- sub("\\..*","",bname,perl=T)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	361
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	362 polii <- read.gff(name,"polii")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	363
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	364 polii.exp <- gene.exp(polii, iter=as.numeric(op.args[["iter"]]))
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	365
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	366 out <- subset(polii.exp,FDR < op.args[["fdr"]])
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	367
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	368 write.table(polii.exp,paste(fname,"genes.details.csv",sep="."),row.names=F,col.names=T,quote=T,sep="\t")
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	369 write.table(out$name, paste(fname,"genes",sep="."),row.names=F,col.names=F,quote=F)
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	370 }
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	371
1b5bd3b955ed planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/damidseq_polII_gene_call commit 2e97fdaeafdafe2c63588de17777cd32a798adb1-dirty mvdbeek parents: diff changeset	372 cat("\nAll done.\n\n")

Mercurial > repos > mvdbeek > damidseq_polii_gene_call

annotate polii.gene.call @ 0:1b5bd3b955ed draft