Fixes: Motif selection & feature matrix construction

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: TFBlearner
 Title: Functionality for training TF-specific classifiers to predict TF bindings based on ATAC-seq data.
-Version: 0.0.1.0001
+Version: 0.0.1.1000
 Authors@R: 
     person("Emanuel", "Sonder", , "emanuel.sonder@hest.ethz.ch", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4788-9508"))

R/contextFeatures.R

@@ -119,7 +119,7 @@
                                               saveHdf5, outDir){
 
                                        data.table::setDTthreads(threads)
-                                       x <- as.integer(round(1000*cor(atacMat, assay(cvSe, NORMDEVASSAY)[motif,])))
+                                       x <- as.integer(round(scaleFactAct*cor(atacMat, assay(cvSe, NORMDEVASSAY)[motif,])))
                                        q <- as.integer(round(quantile(x, prob=c(0,0.1,0.2,0.8,0.9,1), na.rm=TRUE)))
                                        x[abs(x)<sparTh] <- 0L
                                        asc <- Matrix::Matrix(x, ncol=1)

R/getFeatureMatrix.R

@@ -71,10 +71,22 @@
   mats <- lapply(mats, as, "TsparseMatrix")
   js <- cumsum(nCols)
   indDt <- lapply(1:length(mats), function(i){
-      j <- js[[i]]
+      if(i==1){
+        s <- 1
+      }else{
+        s <- js[[i-1]]+1
+      }
+      jNew <- s:js[[i]]
+      jOld <- mats[[i]]@j
+      if(length(jOld)==0){
+        repTimes <- 0}
+      else{
+        repTimes <- table(jOld)
+      }
+      jFull <- rep(jNew, times=repTimes)
       indDt <- data.table(x=mats[[i]]@x,
                           i=mats[[i]]@i+1,
-                          j=rep(j, length(mats[[i]]@x)))})
+                          j=jFull)})
   indDt <- rbindlist(indDt)
 
   # Get matrix
@@ -217,13 +229,13 @@ getFeatureMatrix <- function(mae,
             length(assays(seTfContext))+
             length(assays(seAtac))+
             length(intersect(colnames(colData(seAtac)),
-                   paste(MDSDIMFEATNAME, 1:2, sep="_")))
+                   paste(MDSDIMFEATNAME, 1:2, sep="_")))+
+           sum(addLabels)
 
   if(MAXATACCOLNAME %in% colnames(nonContextTfFeat)){
    nFeats <- nFeats+sum(grepl(INSERTFEATNAME, names(assays(seTfContext))))+
                     sum(TOTALOVERLAPSFEATNAME %in% names(assays(seAtac)))
   }
-  if(addLabels) nFeats <- nFeats+1 # for context-label column
 
   if(saveHdf5)
   {
@@ -326,13 +338,13 @@ getFeatureMatrix <- function(mae,
                                                  colnames(featsNormedMat))]
 
     # normalize by maximum ATAC-signal
-    if(MAXATACCOLNAME %in% colnames(nonContextTfFeat)){
+    if(MAXATACCOLNAME %in% colnames(otherFeatMat)){
       whichCol <- grepl(paste(CONTEXTTFFEAT, INSERTFEATNAME, sep="_"),
                         colnames(featsContextMat))
       countCols <- c(colnames(featsContextMat)[whichCol],
                    paste(CONTEXTFEAT, TOTALOVERLAPSFEATNAME, sep="_"))
       scaledSig <- .minMaxNormalization(featsContextMat[,countCols, drop=FALSE])
-      maxSig <- nonContextTfFeat[,MAXATACCOLNAME, drop=TRUE]
+      maxSig <- otherFeatMat[,MAXATACCOLNAME, drop=TRUE]
       maxScaledMat <- scaledSig / pmax(maxSig, 1e-4)
       colnames(maxScaledMat) <- paste(colnames(maxScaledMat),
                                       NORMEDMAXAFFIX, sep="_")
@@ -380,12 +392,12 @@ getFeatureMatrix <- function(mae,
     else{
       return(featsMat)
     }
-  }, seAtac, seTfContext,
-     tfName, tfCofactors,
-     nonContextTfFeat, norm,
-     saveChunk, hdf5FileName,
-     annoCol,
-     addLabels, convertInteger)
+  }, seAtac=seAtac, seTfContext=seTfContext,
+     tfName=tfName, tfCofactors=tfCofactors,
+     otherFeatMat=nonContextTfFeat, norm=norm,
+     saveChunk=saveChunk, hdf5FileName=hdf5FileName,
+     annoCol=annoCol, addLabels=addLabels,
+     convertInteger=convertInteger)
 
   featMats <- Reduce("rbind", featMats[-1], featMats[[1]])
   featMats <- suppressWarnings({Matrix::Matrix(featMats)})

R/scalingFactors.R

@@ -1,2 +1,3 @@
 scalFactMotif <- 1e4
-scalFactPred <- 1e3
+scaleFactAct <- 1e3
+scalFactPred <- 1e3

R/tfFeatures.R

@@ -141,7 +141,6 @@
 .selectMotifs <- function(matchScores,
                           maxScores,
                           labels,
-                          addThr=4,
                           nMotifs=10,
                           subSample=10000)
 {
@@ -156,10 +155,9 @@
 
   # top motif scores
   matchCoScores <- matchSubScores
-  matchCoScores@x[matchCoScores@x < thr[matchCoScores@j + 1] &
-                    matchCoScores@x<addThr*scalFactMotif] <- 0
-  matchCoScores@x[matchCoScores@x >= thr[matchCoScores@j + 1] &
-                    matchCoScores@x>=addThr*scalFactMotif] <- 1
+  matchCoScores@x[matchCoScores@x < thr[matchCoScores@j + 1]] <- 0
+  matchCoScores@x[matchCoScores@x >= thr[matchCoScores@j + 1]] <- 1
+  matchCoScores <- drop0(matchCoScores)
 
   # get mutually exclusive motif scores
   zeroInd <- which(matchSubScores==0, arr.ind = TRUE)
@@ -168,28 +166,34 @@
                                 x=rep(1, nrow(zeroInd)),
                                 dims=c(nrow(matchSubScores),
                                        ncol(matchSubScores)))
+  colnames(matchExScores) <- colnames(matchSubScores)
 
   # jaccard index of mutually exclusive and top co-occuring motifs
   labels <- matrix(labels, nrow=length(labels), ncol=1)
   matchCo <- .jaccard(matchCoScores, labels)
-  matchCo[,motif_id:=1:.N]
   setorder(matchCo, -cont)
-  topCoMotif <- matchCo$motif_id[1:nMotifs]
+  topCoMotif <- matchCo$set1_col[1:nMotifs]
 
   matchEx <- .jaccard(matchExScores, labels)
-  matchEx[,motif_id:=1:.N]
   setorder(matchEx, -cont)
-  topExMotif <- matchEx$motif_id[1:nMotifs]
+  topExMotif <- matchEx$set1_col[1:nMotifs]
 
-  topExMotif <- intersect(topExMotif, colnames(matchScores))
-  topCoMotif <- intersect(topCoMotif, colnames(matchScores))
   selectedMotifs <- c(topCoMotif, topExMotif)
   if(length(selectedMotifs)>0){
-  names(selectedMotifs) <- c(paste(COMOTIFAFFIX, 1:length(topCoMotif), sep="_"),
-                             paste(EXMOTIFAFFIX, 1:length(topExMotif), sep="_"))}
-
-  # can happen if the motif-matches matrix has less columns than motifs to select
-  selectedMotifs <- unique(selectedMotifs[!is.na(selectedMotifs)])
+     if(length(topCoMotif)>0){
+       namesCo <- paste(COMOTIFAFFIX, 1:length(topCoMotif), sep="_")
+     }
+     else{
+       namesCo <- NULL
+     }
+     if(length(topExMotif)>0){
+       namesEx <- paste(EXMOTIFAFFIX, 1:length(topExMotif), sep="_")
+     }
+    else{
+      namesEx <- NULL
+    }
+    names(selectedMotifs) <- c(namesCo, namesEx)
+  }
 
   return(selectedMotifs)
 }
@@ -470,8 +474,7 @@ tfFeatures <- function(mae,
     names(tfSimMotifCols) <- paste(PRIORMOTIFPREFIX, 1:length(tfSimMotifCols),
                                 sep="_")}
 
-  tfCofactorCols <- unique(grep(paste(tfCofactors,collapse="|"),
-                                motifNames, value=TRUE))
+  tfCofactorCols <- intersect(tfCofactors, motifNames)
   if(length(tfCofactorCols)>0){
     names(tfCofactorCols) <- paste(TFCOFACTORMOTIFPREFIX,
                                    1:length(tfCofactorCols), sep="_")}
@@ -511,7 +514,7 @@ tfFeatures <- function(mae,
     maxScores <- colDataMotifs[[MAXSCORECOL]]
     selMotifs <- .selectMotifs(matchScores, maxScores, labels, nMotifs=nMotifs)
     if(length(selMotifs)>0){
-      names(selMotifs) <- paste0(SELMOTIFPREFIX, names(selMotifs))}
+      names(selMotifs) <- paste(SELMOTIFPREFIX, names(selMotifs), sep=".")}
   }
   else{
     selMotifs <- NULL
@@ -526,8 +529,7 @@ tfFeatures <- function(mae,
     names(tfSimMotifCols) <- paste(PRIORMOTIFPREFIX, 1:length(tfSimMotifCols),
                                    sep="_")}
 
-  tfCofactorCols <- unique(grep(paste(tfCofactors,collapse="|"),
-                                actMotifNames, value=TRUE))
+  tfCofactorCols <- intersect(tfCofactors, actMotifNames)
   if(length(tfCofactorCols)>0){
     names(tfCofactorCols) <- paste(TFCOFACTORMOTIFPREFIX,
                                    1:length(tfCofactorCols), sep="_")}
@@ -558,8 +560,8 @@ tfFeatures <- function(mae,
     # select motifs co-occuring around ChIP-peaks or motif matches of TF of interest
     actAssoc <- assays(mae[[ASSOCEXP]])[[ASSOCASSAY]]
     actAssoc <- actAssoc[,!c(colnames(actAssoc) %in% priorMotifCols), drop=FALSE]
-    selActMotifs <- .selectMotifs(actAssoc, rep(1, ncol(actAssoc)), labels,
-                                  addThr=0, nMotifs=nMotifs)
+    selActMotifs <- .selectMotifs(actAssoc, rep(1*scaleFactAct, ncol(actAssoc)),
+                                  labels, nMotifs=nMotifs)
     if(length(selActMotifs)>0){
       names(selActMotifs) <- paste0(SELMOTIFPREFIX, names(selActMotifs))}
   }