Merge pull request #54 from c-mertes/fraser2

Fraser2

Author: ischeller

R/AllGenerics.R

@@ -788,7 +788,8 @@ FRASER.results <- function(object, sampleIDs, fdrCutoff,
             deltaPsiVals <- deltaPsiValue(tmp_x, type)
             rho          <- rho(tmp_x, type)
             aberrant     <- aberrant.FRASER(tmp_x, type=type, 
-                                            padjCutoff=fdrCutoff, 
+                                            padjCutoff=ifelse(isTRUE(aggregate), 
+                                                        NA, fdrCutoff), 
                                             deltaPsiCutoff=dPsiCutoff, 
                                             minCount=minCount, 
                                             rhoCutoff=rhoCutoff,
@@ -1080,6 +1081,10 @@ aberrant.FRASER <- function(object, type=fitMetrics(object),
     if(is.na(padjCutoff)){
         padjCutoff <- 1
     }
+    if(isTRUE(aggregate)){
+        padjCutoffGene <- padjCutoff
+        padjCutoff <- 1
+    }
 
     if(isTRUE(all)){
         aberrantEvents <- matrix(TRUE, nrow=nrow(object), ncol=ncol(object))
@@ -1121,7 +1126,7 @@ aberrant.FRASER <- function(object, type=fitMetrics(object),
         if(isFALSE(all)){
             aberrantEvents <- aberrantEvents & as.matrix(
                     padj_gene[rownames(aberrantEvents),colnames(aberrantEvents)]
-                ) <= padjCutoff
+                ) <= padjCutoffGene
         }
     }
 

R/getNSetterFuns.R

@@ -748,16 +748,18 @@ getPlottingDT <- function(fds, axis=c("row", "col"), type=currentType(fds),
     )
     dt[, deltaPsi:=obsPsi - predPsi]
 
-    # add aberrant information to it
-    aberrantVec <- aberrant(fds, ..., padjVals=dt[,.(padj)],
-        dPsi=dt[,.(deltaPsi)], n=dt[,.(n)], 
-        rhoVals=dt[,.(rho)], aggregate=FALSE)
-    dt[,aberrant:=aberrantVec]
-
     # if requested return gene p values
     if(isTRUE(aggregate)){
-        dt <- dt[!is.na(featureID)]
+        # get gene-level aberrant status
+        aberrantGeneLevel <- aberrant(fds[, idxcol], ..., aggregate=TRUE)
+        aberrantGeneLevel <- melt(
+                data.table(featureID=rownames(aberrantGeneLevel), 
+                           aberrantGeneLevel), 
+                value.name="aberrant", id.vars="featureID", 
+                variable.name="sampleID")
+        
         # split featureID into several rows if more than one
+        dt <- dt[!is.na(featureID)]
         dt[, dt_idx:=seq_len(.N)]
         dt_tmp <- dt[, splitGenes(featureID), by="dt_idx"]
         dt <- dt[dt_tmp$dt_idx,]
@@ -780,6 +782,10 @@ getPlottingDT <- function(fds, axis=c("row", "col"), type=currentType(fds),
         pvalsGene <- merge(pvalsGene[[1]], pvalsGene[[2]], 
                             by=c("featureID", "sampleID"))
 
+        # merge with gene level aberrant status
+        pvalsGene <- merge(pvalsGene, aberrantGeneLevel, 
+                           by=c("featureID", "sampleID"))
+        
         # merge with gene pval matrix
         dt <- merge(dt, pvalsGene, by=c("featureID", "sampleID"))
         dt[,`:=`(pval=gene_pval, padj=gene_padj, 
@@ -789,6 +795,12 @@ getPlottingDT <- function(fds, axis=c("row", "col"), type=currentType(fds),
         dt <- dt[order(sampleID, featureID, type, -aberrant,
                         padj, -abs(deltaPsi))][
                 !duplicated(data.table(sampleID, featureID, type))]
+    } else{
+        # add aberrant information to it
+        aberrantVec <- aberrant(fds, ..., padjVals=dt[,.(padj)],
+                                dPsi=dt[,.(deltaPsi)], n=dt[,.(n)], 
+                                rhoVals=dt[,.(rho)], aggregate=FALSE)
+        dt[,aberrant:=aberrantVec]
     }
 
     # return object

R/mergeExternalData.R

@@ -100,8 +100,8 @@ mergeExternalData <- function(fds, countFiles, sampleIDs, annotation=NULL){
     # merge psi5/psi3 data
     #
     extractExtData <- function(fds, countFun, type, ov, extData, extName){
-        ctsOri <- as.matrix(countFun(fds, type=type)[from(ov),])
-        ctsExt <- as.matrix(mcols(extData[[extName]])[to(ov),])
+        ctsOri <- as.matrix(countFun(fds, type=type)[from(ov),,drop=FALSE])
+        ctsExt <- as.matrix(mcols(extData[[extName]])[to(ov),,drop=FALSE])
         ans <- cbind(ctsOri, ctsExt)
         mode(ans) <- "integer"
         ans

R/plotMethods.R

@@ -337,25 +337,6 @@ plotVolcano.FRASER <- function(object, sampleID,
         theme(legend.position="none") +
         scale_color_manual(values=c("gray40", "firebrick"))
 
-    if(!is.na(deltaPsiCutoff)){
-        g <- g + 
-            geom_vline(xintercept=c(-deltaPsiCutoff, deltaPsiCutoff),
-                    color="firebrick", linetype=2)
-    }
-    
-    if(!is.na(padjCutoff)){
-        if(dt[padj <= padjCutoff, .N] > 0){
-            padj_line <- min(dt[padj <= padjCutoff, -log10(pval)])
-            padj_line <- min(dt[padj <= padjCutoff, -log10(pval)])
-        }
-        if(!"padj_line" %in% ls() || padj_line > 10 || is.na(padj_line)){
-            padj_line <- 6
-        }
-        g <- g + 
-            geom_hline(yintercept=padj_line, color="firebrick", linetype=4)
-    }
-    
-    
     if(isFALSE(basePlot)){
         g <- g + xlab(paste("delta", 
                             ggplotLabelPsi(type, asCharacter=TRUE)[[1]])) +