diff corr.R @ 11:74bfa2464411

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
line wrap: on
line diff
--- a/corr.R	Sun Aug 18 18:10:48 2013 +0200
+++ b/corr.R	Sun Aug 18 19:39:52 2013 +0200
@@ -133,6 +133,26 @@
     }
 }
 
+#I hate this bioconductor ranges stuff more and more
+#this is a super s**t api
+mergeSeqInfo <- function(coverages){
+  seqinfos <- do.call("rbind", lapply(coverages, function(cov){ d <- as.data.frame(cov@seqinfo); d$seqnames <- rownames(d); d}))
+  seqinfosMax <- do.call("rbind", by(seqinfos, seqinfos$seqnames, function(byChr){ byChr[order(byChr$seqlengths, decreasing=TRUE),][1,] }))
+  data.frame(seqnames=seqinfosMax$seqnames,seqlengths=seqinfosMax$seqlengths, isCircular=seqinfosMax$isCircular, genome=seqinfosMax$genome)
+}
+
+
+changeSeqInfo <- function(coverages){
+  seqinfo <- mergeSeqInfo(coverages)
+  lapply(coverages, function(cov){ 
+	  si <- as.data.frame(cov@seqinfo)
+      si <- data.frame(chr=rownames(si), si)
+	  sis <- subset(seqinfo, seqnames %in% si$chr) 
+	  si <- Seqinfo(seqnames=as.character(sis$seqnames), seqlengths=sis$seqlengths, isCircular=sis$isCircular, genome=sis$genome)   
+      cov@seqinfo <- si
+	  cov
+  })
+}
 
 args <- commandArgs(TRUE)
 
@@ -151,8 +171,11 @@
 
 
 coverages <- getCoverages(infiles, formats)
+coverages <- changeSeqInfo(coverages)
+
 mappability <- getMappable(coverages, mappable)
 
+
 calcAndPlot(coverages, outnames, mappability, outnamePDF, outnameMat,title)