spikein norm in csaw working in allele-specific mode

Author: katarzyna.otylia.sikora@gmail.com

snakePipes/shared/rscripts/CSAW.R

@@ -39,8 +39,10 @@ message(paste("FDR:", fdr, "\n"))
 message(paste("LFC:", lfc, "\n"))
 message(paste("paired-end? :", pairedEnd, "\n"))
 message(paste("allele-specific? :", allelic_info, "\n"))
+message(paste("spikein normalization requested :", useSpikeInForNorm, "\n"))
 message(paste("External bed? ;", external_bed, "\n"))
 
+
 ## sampleInfo (setup of the experiment)
 sampleInfo <- read.table(sampleInfoFilePath, header = TRUE, colClasses = c("character", "character"))
 rownames(sampleInfo)<-sampleInfo$name
@@ -128,6 +130,9 @@ if (! external_bed) {
         })
     }
     # merge
+    all_levels<-unique(unlist(lapply(allpeaks,function(X)seqlevels(X))))
+    allpeaks<-lapply(allpeaks,function(X){seqlevels(X)<-all_levels
+                                          return(X)})
     allpeaks <- Reduce(function(x,y) GenomicRanges::union(x,y), allpeaks)
 
     } else {

snakePipes/shared/rscripts/DB_functions.R

@@ -121,7 +121,6 @@ makeQCplots_chip_PE <- function(bam.file, outplot, pe.param){
     }
     collected <- plotwc(bam.file)
     xranged <- as.integer(names(collected))
-
     ## plot
     message("Plotting")
     pdf(outplot)
@@ -136,12 +135,13 @@ makeQCplots_chip_PE <- function(bam.file, outplot, pe.param){
 
     # cross correlation
     plot(0:max.delay, CCF, type = "l", ylab = "CCF", xlab = "Delay (bp)", main = "PE-Cross-correlation")
-
+        print(str(collected))
+    if( any(is.null(collected)) | any(is.na(collected))){message("Skipping plotting cross correlation.")}else{
     # coverage in windows
     plot(xranged, collected, type = "l", col = "blue", xlim = c(-1000, 1000), lwd = 2,
-         xlab = "Distance (bp)", ylab = "Relative coverage per base")
+           xlab = "Distance (bp)", ylab = "Relative coverage per base")
     abline(v = c(-150,200), col = "dodgerblue", lty = 2)
-    legend("topright", col = "dodgerblue", legend = "specified window size")
+    legend("topright", col = "dodgerblue", legend = "specified window size")}
 
     dev.off()
 
@@ -204,7 +204,11 @@ tmmNormalize_chip <- function(chipCountObject, binsize, plotfile){
     wider <- csaw::windowCounts(bam.files, bin = TRUE, width = binsize, param = chipCountObject$pe.param)
     if(useSpikeInForNorm){
         tab<-read.table(scale_factors,sep="\t",header=TRUE,as.is=TRUE,quote="")
-        normfacs<-1/(tab$scalingFactor[match(colnames(chipCountObject$windowCounts),tab$sample)]) }else{
+        if(! allelic_info){
+            normfacs<-1/(tab$scalingFactor[match(colnames(chipCountObject$windowCounts),tab$sample)])}else{
+            normfacs<-1/(tab$scalingFactor[match(gsub(".genome[1-2]","",colnames(chipCountObject$windowCounts)),tab$sample)])
+            }
+     }else{
         normfacs <- csaw::normFactors(wider, se.out=FALSE)}
 
     chipCountObject$normFactors <- normfacs