Browse code

fix thresholding

Davin Brownell authored on 06/01/2022 21:01:42
Showing 1 changed files

... ...
@@ -331,20 +331,60 @@ function(object, threshold, lp, axis)
331 331
     ## keep only a subset of markers for each pattern depending on the type of threshold
332 332
     if (threshold == "cut") # all markers which achieve minimum rank
333 333
     {
334
-        rankCutoff <- sapply(1:ncol(markerRanks), function(patternIndex)
335
-        {
336
-            patternRank <- markerRanks[,patternIndex]
337
-            return(max(patternRank[patternRank == apply(markerRanks, 1, min)]))
338
-        })
339
-        markersByPattern <- lapply(1:ncol(markerRanks), function(patternIndex)
340
-            rownames(markerRanks)[markerRanks[,patternIndex] <= rankCutoff[patternIndex]])
341
-        names(markersByPattern) <- colnames(markerRanks)
334
+      simplicityGENES <- function(As, Ps) {
335
+        # rescale p's to have max 1
336
+        pscale <- apply(Ps,1,max)
337
+        
338
+        # rescale A in accordance with p's having max 1
339
+        As <- sweep(As, 2, pscale, FUN="*")
340
+        
341
+        # find the A with the highest magnitude
342
+        Arowmax <- t(apply(As, 1, function(x) x/max(x)))
343
+        pmax <- apply(As, 1, max)
344
+        
345
+        # determine which genes are most associated with each pattern
346
+        ssl <- matrix(NA, nrow=nrow(As), ncol=ncol(As),
347
+                      dimnames=dimnames(As))
348
+        for (i in 1:ncol(As)) {
349
+          lp <- rep(0, ncol(As))
350
+          lp[i] <- 1
351
+          ssl.stat <- apply(Arowmax, 1, function(x) sqrt(t(x-lp)%*%(x-lp)))
352
+          ssl[order(ssl.stat),i] <- 1:length(ssl.stat)
353
+        }
354
+        
355
+        return(ssl)
356
+        
357
+      }
358
+      simGenes <- simplicityGENES(As=object@featureLoadings,
359
+                                  Ps=object@sampleFactors)
360
+      
361
+      patternMarkers <- list()
362
+      
363
+      nP <- ncol(simGenes)
364
+      
365
+      for (i in 1:nP) {
366
+        
367
+        sortSim <- names(sort(simGenes[,i],decreasing=F))
368
+        
369
+        geneThresh <- min(which(simGenes[sortSim,i] > 
370
+                                  apply(simGenes[sortSim,],1,min)))
371
+        
372
+        markerGenes <- sortSim[1:geneThresh]
373
+        
374
+        markerGenes <- unique(markerGenes)
375
+        
376
+        patternMarkers[[i]] <- markerGenes
377
+        
378
+        markersByPattern <- patternMarkers
379
+        
380
+      }
342 381
     }
343 382
     else if (threshold == "all") # only the markers with the lowest scores
344 383
     {
345 384
         patternsByMarker <- colnames(markerScores)[apply(markerScores, 1, which.min)]
346 385
         markersByPattern <- sapply(colnames(markerScores), USE.NAMES=TRUE, simplify=FALSE,
347 386
             function(pattern) rownames(markerScores)[which(patternsByMarker==pattern)])
387
+        
348 388
     }
349 389
     return(list(
350 390
         "PatternMarkers"=markersByPattern,