Mercurial > repos > idot > coverage_correlation

comparison corr.R @ 13:de452b96da8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added bam input
author Ido Tamir <ido.tamir@imp.ac.at>
date Mon, 23 Sep 2013 13:46:24 +0200
parents f3e037496c18
children
comparison
equal deleted inserted replaced
12:f3e037496c18 13:de452b96da8e
15 #args[6] <- title 15 #args[6] <- title
16 #args[7] <- optional mappability bed file #not implemented because we have to change the coverage function (covWith0) to 1. intersection then diff because mappability can be wrong 16 #args[7] <- optional mappability bed file #not implemented because we have to change the coverage function (covWith0) to 1. intersection then diff because mappability can be wrong
17 17
18 options("useFancyQuotes" = FALSE) 18 options("useFancyQuotes" = FALSE)
19 19
20 TEST=FALSE
21
20 library(rtracklayer) 22 library(rtracklayer)
21 library(lattice) 23 library(lattice)
22 library(latticeExtra) 24 library(latticeExtra)
23 library(hexbin) 25 library(hexbin)
26 library(Rsamtools)
24 27
25 ## creates the union of coverages i.e. mappable regions 28 ## creates the union of coverages i.e. mappable regions
26 ## todo: add mappability argument 29 ## todo: add mappability argument
27 createMappable <- function(coverages){ 30 createMappable <- function(coverages){
28 suppressWarnings(Reduce(union, coverages)) #supress because missing chromosomes spit warinings 31 suppressWarnings(Reduce(union, coverages)) #supress because missing chromosomes spit warinings
107 110
108 } 111 }
109 112
110 getCoverage <- function(infile, format){ 113 getCoverage <- function(infile, format){
111 print(paste("reading", infile, format)) 114 print(paste("reading", infile, format))
112 values <- import(infile, format=format, asRangedData = FALSE) 115 values <- if(format == "bam"){
116 import(infile, format=format)
117 }else{
118 import(infile, format=format, asRangedData = FALSE)
119 }
113 if(tolower(format) %in% c("bigwig","wig")){ 120 if(tolower(format) %in% c("bigwig","wig")){
114 values 121 values
115 }else{ 122 }else{
116 cov <- coverage(values) 123 cov <- coverage(values)
117 gcov <- as(cov, "GRanges") 124 gcov <- as(cov, "GRanges")
155 }) 162 })
156 } 163 }
157 164
158 args <- commandArgs(TRUE) 165 args <- commandArgs(TRUE)
159 166
167 if(TEST){
168 args[1] <- c("file1.bed,file2.bed,file3.bed,file4.bed")
169 args[1] <- c("file1.bam,file2.bam,file3.bam,file4.bam")
170 args[2] <- c("file1,file2,file3,file4")
171 args[3] <- c("bam,bam,bam,bam")
172 args[4] <- "out.pdf"
173 args[5] <- "mat.mat"
174 args[6] <- "title"
175 args[7] <- NA
176 }
177
160 infiles <- unlist(strsplit(args[1], ",")) 178 infiles <- unlist(strsplit(args[1], ","))
161 outnames <- unlist(strsplit(args[2], ",")) 179 outnames <- unlist(strsplit(args[2], ","))
162 formats <- unlist(strsplit(args[3], ",")) 180 formats <- unlist(strsplit(args[3], ","))
163 outnamePDF <- args[4] 181 outnamePDF <- args[4]
164 outnameMat <- args[5] 182 outnameMat <- args[5]