Mercurial > repos > idot > coverage_correlation

comparison corr.R @ 11:74bfa2464411

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
the bioconductor GRanges stuff blows
author Ido Tamir <ido.tamir@imp.ac.at>
date Sun, 18 Aug 2013 19:39:52 +0200
parents ecad2cd5038f
children f3e037496c18
comparison
equal deleted inserted replaced
10:6df7b07f81c1 11:74bfa2464411
131 }else{ 131 }else{
132 import(mappable, format="bed", asRangedData = FALSE) 132 import(mappable, format="bed", asRangedData = FALSE)
133 } 133 }
134 } 134 }
135 135
136 #I hate this bioconductor ranges stuff more and more
137 #this is a super s**t api
138 mergeSeqInfo <- function(coverages){
139 seqinfos <- do.call("rbind", lapply(coverages, function(cov){ d <- as.data.frame(cov@seqinfo); d$seqnames <- rownames(d); d}))
140 seqinfosMax <- do.call("rbind", by(seqinfos, seqinfos$seqnames, function(byChr){ byChr[order(byChr$seqlengths, decreasing=TRUE),][1,] }))
141 data.frame(seqnames=seqinfosMax$seqnames,seqlengths=seqinfosMax$seqlengths, isCircular=seqinfosMax$isCircular, genome=seqinfosMax$genome)
142 }
143
144
145 changeSeqInfo <- function(coverages){
146 seqinfo <- mergeSeqInfo(coverages)
147 lapply(coverages, function(cov){
148 si <- as.data.frame(cov@seqinfo)
149 si <- data.frame(chr=rownames(si), si)
150 sis <- subset(seqinfo, seqnames %in% si$chr)
151 si <- Seqinfo(seqnames=as.character(sis$seqnames), seqlengths=sis$seqlengths, isCircular=sis$isCircular, genome=sis$genome)
152 cov@seqinfo <- si
153 cov
154 })
155 }
136 156
137 args <- commandArgs(TRUE) 157 args <- commandArgs(TRUE)
138 158
139 infiles <- unlist(strsplit(args[1], ",")) 159 infiles <- unlist(strsplit(args[1], ","))
140 outnames <- unlist(strsplit(args[2], ",")) 160 outnames <- unlist(strsplit(args[2], ","))
149 print("parsed input parameters") 169 print("parsed input parameters")
150 print(lpp) 170 print(lpp)
151 171
152 172
153 coverages <- getCoverages(infiles, formats) 173 coverages <- getCoverages(infiles, formats)
174 coverages <- changeSeqInfo(coverages)
175
154 mappability <- getMappable(coverages, mappable) 176 mappability <- getMappable(coverages, mappable)
177
155 178
156 calcAndPlot(coverages, outnames, mappability, outnamePDF, outnameMat,title) 179 calcAndPlot(coverages, outnames, mappability, outnamePDF, outnameMat,title)
157 180
158 181
159 182