Browse code

fix thresholding

Davin Brownell authored on 06/01/2022 21:01:42
Showing1 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,
Browse code

Merge pull request #57 from FertigLab/plot-pattern-markers-fix

plotPatternMarkers function fix

Jeanette Johnson authored on 19/10/2021 18:33:03 • GitHub committed on 19/10/2021 18:33:03
Showing0 changed files
Browse code

oops left in a hardcoded number

jeanettejohnson authored on 15/10/2021 20:44:13
Showing1 changed files
... ...
@@ -77,7 +77,7 @@ plot.CogapsResult <- function(x, groups=NULL, ...)
77 77
     ylimits <- c(0, max(x@sampleFactors) * 1.1)
78 78
     plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude", xlab="", axes=FALSE) 
79 79
     axis(1, labels=FALSE, at=1:length(grpset))
80
-    text(x = seq(1, 12, by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)
80
+    text(x = seq(1, length(grpset), by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)
81 81
     axis(2)
82 82
     
83 83
     for (i in 1:nFactors)
Browse code

basic plot function takes optional groups param for large datasets

jeanettejohnson authored on 15/10/2021 20:08:28
Showing1 changed files
... ...
@@ -64,25 +64,50 @@ function(object)
64 64
 #' @export
65 65
 #' @importFrom graphics plot legend lines points
66 66
 #' @importFrom grDevices rainbow
67
-plot.CogapsResult <- function(x, ...)
67
+plot.CogapsResult <- function(x, groups=NULL, ...)
68 68
 {
69
+  
70
+  if(!is.null(groups)) {
71
+    grpset <- unique(groups)
72
+    groups <- table(groups)
73
+    ngroups <- length(grpset)
74
+    nFactors <- ncol(x@sampleFactors)
75
+    colors <- rainbow(nFactors)
76
+    xlimits <- c(0, ngroups + 1)
77
+    ylimits <- c(0, max(x@sampleFactors) * 1.1)
78
+    plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude", xlab="", axes=FALSE) 
79
+    axis(1, labels=FALSE, at=1:length(grpset))
80
+    text(x = seq(1, 12, by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)
81
+    axis(2)
82
+    
83
+    for (i in 1:nFactors)
84
+    {
85
+      lines(x=1:length(groups), y=x@sampleFactors[groups,i], col=colors[i])
86
+      points(x=1:length(groups), y=x@sampleFactors[groups,i], col=colors[i], pch=i)
87
+    }
88
+    
89
+    legend("top", paste("Pattern", 1:nFactors, sep = ""), pch = 1:nFactors,
90
+           lty=1, cex=0.8, col=colors, bty="y", ncol=5)
91
+  }
92
+  else{
69 93
     nSamples <- nrow(x@sampleFactors)
70 94
     nFactors <- ncol(x@sampleFactors)
71 95
     colors <- rainbow(nFactors)
72 96
     xlimits <- c(0, nSamples + 1)
73 97
     ylimits <- c(0, max(x@sampleFactors) * 1.1)
74
-
98
+    
75 99
     plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude",
76
-        xlab="Samples")
77
-
100
+         xlab="Samples")
101
+    
78 102
     for (i in 1:nFactors)
79 103
     {
80
-        lines(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i])
81
-        points(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i], pch=i)
104
+      lines(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i])
105
+      points(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i], pch=i)
82 106
     }
83
-
107
+    
84 108
     legend("top", paste("Pattern", 1:nFactors, sep = ""), pch = 1:nFactors,
85
-        lty=1, cex=0.8, col=colors, bty="y", ncol=5)
109
+           lty=1, cex=0.8, col=colors, bty="y", ncol=5)
110
+  }
86 111
 }
87 112
 
88 113
 #' @rdname getFeatureLoadings-methods
Browse code

fixed plotPatternMarkers function?

Jeanette Johnson authored on 21/07/2021 17:10:33
Showing1 changed files
... ...
@@ -437,10 +437,12 @@ function(object, GStoGenes, numPerm, Pw, PwNull)
437 437
 #' @param data the original data as a matrix
438 438
 #' @param patternPalette a vector indicating what color should be used
439 439
 #' for each pattern
440
+#' @param patternMarkers pattern markers to be plotted, as generated by the 
441
+#' patternMarkers function
440 442
 #' @param sampleNames names of the samples to use for labeling 
441 443
 #' @param samplePalette  a vector indicating what color should be used
442 444
 #' for each sample
443
-#' @param colDenogram logical indicating whether to display sample denogram
445
+#' @param colDendrogram logical indicating whether to display sample dendrogram
444 446
 #' @param heatmapCol pallelet giving color scheme for heatmap
445 447
 #' @param scale character indicating if the values should be centered and scaled
446 448
 #' in either the row direction or the column direction, or none. The
... ...
@@ -451,8 +453,8 @@ function(object, GStoGenes, numPerm, Pw, PwNull)
451 453
 #' @importFrom gplots bluered
452 454
 #' @importFrom gplots heatmap.2
453 455
 #' @importFrom stats hclust as.dist cor
454
-plotPatternMarkers <- function(object, data, patternPalette, sampleNames,
455
-samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
456
+plotPatternMarkers <- function(object, data, patternMarkers, patternPalette, sampleNames,
457
+samplePalette=NULL, heatmapCol=bluered, colDendrogram=TRUE, scale="row", ...)
456 458
 {
457 459
     if (is.null(samplePalette))
458 460
         samplePalette <- rep("black", ncol(data))
... ...
@@ -462,10 +464,10 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
462 464
     names(patternCols) <- unlist(patternMarkers)
463 465
     for (i in 1:length(patternMarkers))
464 466
     {
465
-        patternCols[patternMarkers[[i]]] <- patternPalette[i]
467
+        patternCols[unlist(patternMarkers$PatternMarkers[[i]])] <- patternPalette[i]
466 468
     }
467 469
 
468
-    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),
470
+    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers$PatternMarkers),],),
469 471
         symkey=FALSE,
470 472
         symm=FALSE, 
471 473
         scale=scale,
... ...
@@ -473,9 +475,9 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
473 475
         distfun=function(x) as.dist((1-cor(t(x)))/2),
474 476
         hclustfun=function(x) stats::hclust(x,method="average"),
475 477
         Rowv=FALSE,
476
-        Colv=colDenogram,
478
+        Colv=colDendrogram,
477 479
         trace='none',
478
-        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),
480
+        RowSideColors = as.character(patternCols[unlist(patternMarkers$PatternMarkers)]),
479 481
         labCol= sampleNames,
480 482
         cexCol=0.8,
481 483
         ColSideColors=as.character(samplePalette),
Browse code

Fix some bugs in the documentation

sherman5 authored on 30/06/2020 05:53:52
Showing1 changed files
... ...
@@ -97,8 +97,6 @@ function(object)
97 97
 
98 98
 #' @rdname getAmplitudeMatrix-methods
99 99
 #' @aliases getAmplitudeMatrix
100
-#' The method is same as getFeatureLoadings
101
-#' Alternate name for getFeatureLoadings
102 100
 setMethod("getAmplitudeMatrix", signature(object="CogapsResult"),
103 101
 function(object)
104 102
 {
... ...
@@ -117,8 +115,6 @@ function(object)
117 115
 
118 116
 #' @rdname getPatternMatrix-methods
119 117
 #' @aliases getPatternMatrix
120
-#' The method is same as getSampleFactors
121
-#' Alternate name for getSampleFactors
122 118
 setMethod("getPatternMatrix", signature(object="CogapsResult"),
123 119
 function(object)
124 120
 {
Browse code

Added methods getAmplitudeMatrix and getPatternMatrix

Archana Balan authored on 12/11/2019 17:10:54
Showing1 changed files
... ...
@@ -93,6 +93,20 @@ function(object)
93 93
     object@featureLoadings
94 94
 })
95 95
 
96
+
97
+
98
+#' @rdname getAmplitudeMatrix-methods
99
+#' @aliases getAmplitudeMatrix
100
+#' The method is same as getFeatureLoadings
101
+#' Alternate name for getFeatureLoadings
102
+setMethod("getAmplitudeMatrix", signature(object="CogapsResult"),
103
+function(object)
104
+{
105
+    object@featureLoadings
106
+})
107
+
108
+
109
+
96 110
 #' @rdname getSampleFactors-methods
97 111
 #' @aliases getSampleFactors
98 112
 setMethod("getSampleFactors", signature(object="CogapsResult"),
... ...
@@ -101,6 +115,17 @@ function(object)
101 115
     object@sampleFactors
102 116
 })
103 117
 
118
+#' @rdname getPatternMatrix-methods
119
+#' @aliases getPatternMatrix
120
+#' The method is same as getSampleFactors
121
+#' Alternate name for getSampleFactors
122
+setMethod("getPatternMatrix", signature(object="CogapsResult"),
123
+function(object)
124
+{
125
+    object@sampleFactors
126
+})
127
+
128
+
104 129
 #' @rdname getMeanChiSq-methods
105 130
 #' @aliases getMeanChiSq
106 131
 setMethod("getMeanChiSq", signature(object="CogapsResult"),
Browse code

Fix up remaining errors during package check

Tom Sherman authored on 09/09/2019 23:42:33
Showing1 changed files
... ...
@@ -315,7 +315,7 @@ function(object, sets, whichMatrix, numPerm, ...)
315 315
 {
316 316
     ## do this for backwards compatibility
317 317
     if (!is.null(list(...)$GStoGenes))
318
-        sets <- GStoGenes
318
+        sets <- list(...)$GStoGenes
319 319
     ## check input arguments
320 320
     if (!is(sets, "list"))
321 321
         stop("sets must be a list of either measurements or samples")
Browse code

Change std dev slots in CogapsResult

In order to be more consistent with terminology, featureStdDev has
become loadingStdDev, and sampleStdDev has become factorStdDev.

Tom Sherman authored on 05/09/2019 20:59:43
Showing1 changed files
... ...
@@ -56,7 +56,6 @@ function(object)
56 56
     nFeatures <- nrow(object@featureLoadings)
57 57
     nSamples <- nrow(object@sampleFactors)
58 58
     nPatterns <- ncol(object@featureLoadings)
59
-
60 59
     print(paste("CogapsResult object with", nFeatures, "features and", nSamples,
61 60
         "samples"))
62 61
     print(paste(nPatterns, "patterns were learned"))
... ...
@@ -133,7 +132,6 @@ function(object)
133 132
 {
134 133
     if (!is.null(object@metadata$unmatchedPatterns))
135 134
         return(object@metadata$unmatchedPatterns)
136
-
137 135
     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
138 136
     return(NULL)
139 137
 })
... ...
@@ -145,7 +143,6 @@ function(object)
145 143
 {
146 144
     if (!is.null(object@metadata$clusteredPatterns))
147 145
         return(object@metadata$clusteredPatterns)
148
-
149 146
     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
150 147
     return(NULL)
151 148
 })
... ...
@@ -157,7 +154,6 @@ function(object)
157 154
 {
158 155
     if (!is.null(object@metadata$CorrToMeanPattern))
159 156
         return(object@metadata$CorrToMeanPattern)
160
-
161 157
     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
162 158
     return(NULL)
163 159
 })
... ...
@@ -169,7 +165,6 @@ function(object)
169 165
 {
170 166
     if (!is.null(object@metadata$subsets))
171 167
         return(object@metadata$subsets)
172
-
173 168
     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
174 169
     return(NULL)
175 170
 })
... ...
@@ -179,28 +174,14 @@ function(object)
179 174
 setMethod("calcZ", signature(object="CogapsResult"),
180 175
 function(object, whichMatrix)
181 176
 {
182
-    if (whichMatrix=="featureLoadings")
183
-    {
184
-        if (sum(object@featureStdDev==0) > 0)
185
-        {
186
-            warning("zeros detected in the standard deviation matrix")
187
-            object@featureStdDev[object@featureStdDev==0] <- 1e-6
188
-        }
189
-        return(object@featureLoadings / object@featureStdDev)
190
-    }
191
-    else if (whichMatrix=="sampleFactors")
192
-    {
193
-        if (sum(object@sampleStdDev==0) > 0)
194
-        {
195
-            warning("zeros detected in the standard deviation matrix")
196
-            object@sampleStdDev[object@sampleStdDev==0] <- 1e-6
197
-        }
198
-        return(object@sampleFactors / object@sampleStdDev)
199
-    }
200
-    else
201
-    {
177
+    if (!(whichMatrix %in% c("featureLoadings", "sampleFactors")))
202 178
         stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")
203
-    }
179
+    mean <- if (whichMatrix == "sampleFactors") object@sampleFactors else object@featureLoadings
180
+    stddev <- if (whichMatrix == "sampleFactors") object@factorStdDev else object@loadingStdDev
181
+    if (sum(stddev == 0) > 0)
182
+        warning("zeros detected in the standard deviation matrix")
183
+    stddev[stddev == 0] <- 1e-6
184
+    return(mean / stddev)    
204 185
 })
205 186
 
206 187
 #' @rdname reconstructGene-methods
... ...
@@ -384,11 +365,11 @@ function(object, GStoGenes, numPerm, Pw, nullGenes)
384 365
     if (nullGenes)
385 366
     {
386 367
         ZD <- object@featureLoadings[setdiff(row.names(object@featureLoadings), GStoGenes),] /
387
-            object@featureStdDev[setdiff(row.names(object@featureLoadings), GStoGenes),]
368
+            object@loadingStdDev[setdiff(row.names(object@featureLoadings), GStoGenes),]
388 369
     }
389 370
     else
390 371
     {
391
-        ZD <- object@featureLoadings[GStoGenes,]/object@featureStdDev[GStoGenes,]
372
+        ZD <- object@featureLoadings[GStoGenes,]/object@loadingStdDev[GStoGenes,]
392 373
     }
393 374
     outStats <- apply(sweep(ZD,2,gsStat,FUN="*"),1,sum) / (sum(gsStat))
394 375
     outStats <- outStats / apply(ZD,1,sum)
Browse code

tidied up notation around pattern markers

Tom Sherman authored on 05/09/2019 18:14:53
Showing1 changed files
... ...
@@ -258,127 +258,73 @@ function(object, data, uncertainty)
258 258
         lwid=c(1,10), lhei=c(1, 4, 1.2 ), main="Heatmap of Residuals")
259 259
 })
260 260
 
261
-patternMarkersCalculation <- function(mat, threshold, lp)
261
+unitVector <- function(n, length)
262 262
 {
263
+    vec <- rep(0, length)
264
+    vec[n] <- 1
265
+    return(vec)
266
+}
267
+
268
+#' @rdname patternMarkers-methods
269
+#' @aliases patternMarkers
270
+setMethod("patternMarkers", signature(object="CogapsResult"),
271
+function(object, threshold, lp, axis)
272
+{
273
+    ## check inputs to the function
263 274
     if (!(threshold %in% c("cut", "all")))
264
-        stop("invalid thresold argument")
265
-    if (!is.na(lp) & length(lp) != ncol(mat))
275
+        stop("threshold must be either 'cut' or 'all'")
276
+    if (!is.na(lp) & length(lp) != ncol(object@featureLoadings))
266 277
         stop("lp length must equal the number of patterns")
267
-    
268
-    nMeasurements <- nrow(mat)
269
-    nPatterns <- ncol(mat)
270
-    rowMax <- t(apply(mat, 1, function(x) x / max(x)))
278
+    if (!(axis %in% 1:2))
279
+        stop("axis must be either 1 or 2")
280
+    ## need to scale each row of the matrix of interest so that the maximum is 1
281
+    resultMatrix <- if (axis == 1) object@featureLoadings else object@sampleFactors
282
+    normedMatrix <- t(apply(resultMatrix, 1, function(row) row / max(row)))
283
+    ## handle the case where a custom linear combination of patterns was passed in
271 284
     if (!is.na(lp))
272 285
     {
273
-        sstat <- apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp)))
274
-        ssranks <- rank(sstat)
275
-        ssgenes.th <- names(sort(sstat,decreasing=FALSE,na.last=TRUE))
286
+        markerScores <- apply(normedMatrix, 1, function(row) sqrt(sum((row-lp)^2)))
287
+        markersByPattern <- names(sort(markerScores, decreasing=FALSE, na.last=TRUE))
288
+        return(list(
289
+            "PatternMarkers"=markersByPattern,
290
+            "PatternMarkerRanks"=rank(markerScores),
291
+            "PatternMarkerScores"=markerScores
292
+        ))
276 293
     }
277
-    else
278
-    {
279
-        # determine which measurements are most associated with each pattern
280
-        sstat <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))
281
-        ssranks <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))
282
-        ssgenes <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=NULL)
283
-        for (i in 1:nPatterns)
294
+    ## default pattern marker calculation, each pattern has unit weight
295
+    markerScores <- sapply(1:ncol(normedMatrix), function(patternIndex)
296
+        apply(normedMatrix, 1, function(row)
284 297
         {
285
-            lp <- rep(0,nPatterns)
286
-            lp[i] <- 1
287
-            sstat[,i] <- unlist(apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp))))
288
-            ssranks[,i] <- rank(sstat[,i])
289
-            ssgenes[,i] <- names(sort(sstat[,i], decreasing=FALSE, na.last=TRUE))
290
-        }
291
-        if (threshold=="cut")
292
-        {
293
-            geneThresh <- sapply(1:nPatterns, function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))
294
-            ssgenes.th <- sapply(1:nPatterns, function(x) ssgenes[1:geneThresh[x],x])
295
-        }
296
-        else if (threshold=="all")
298
+            lp <- unitVector(patternIndex, ncol(normedMatrix))
299
+            return(sqrt(sum((row-lp)^2)))
300
+        })
301
+    )
302
+    markerRanks <- apply(markerScores, 2, rank)
303
+    colnames(markerScores) <- colnames(markerRanks) <- colnames(normedMatrix)
304
+    ## keep only a subset of markers for each pattern depending on the type of threshold
305
+    if (threshold == "cut") # all markers which achieve minimum rank
306
+    {
307
+        rankCutoff <- sapply(1:ncol(markerRanks), function(patternIndex)
297 308
         {
298
-            pIndx <- unlist(apply(sstat, 1, which.min))
299
-            gBYp <- list()
300
-            for (i in sort(unique(pIndx)))
301
-            {
302
-                gBYp[[i]] <- sapply(strsplit(names(pIndx[pIndx==i]),"[.]"),function(x) x[[1]][1])
303
-            }
304
-            ssgenes.th <- lapply(1:max(sort(unique(pIndx))), function(x)
305
-                ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])
306
-        }
309
+            patternRank <- markerRanks[,patternIndex]
310
+            return(max(patternRank[patternRank == apply(markerRanks, 1, min)]))
311
+        })
312
+        markersByPattern <- lapply(1:ncol(markerRanks), function(patternIndex)
313
+            rownames(markerRanks)[markerRanks[,patternIndex] <= rankCutoff[patternIndex]])
314
+        names(markersByPattern) <- colnames(markerRanks)
307 315
     }
308
-    return(list("PatternMarkers"=ssgenes.th, "PatternRanks"=ssranks,
309
-        "PatternMarkerScores"=sstat))       
310
-}
311
-
312
-#' @rdname patternMarkers-methods
313
-#' @aliases patternMarkers
314
-setMethod("patternMarkers", signature(object="CogapsResult"),
315
-function(object, threshold, lp)
316
-{
317
-    patternMarkersCalculation(object@featureLoadings, threshold, lp)
318
-})
319
-
320
-#' @rdname amplitudeMarkers-methods
321
-#' @aliases amplitudeMarkers
322
-setMethod("amplitudeMarkers", signature(object="CogapsResult"),
323
-function(object, threshold, lp)
324
-{
325
-    patternMarkersCalculation(object@sampleFactors, threshold, lp)
326
-})
327
-
328
-#' heatmap of original data clustered by pattern markers statistic
329
-#' @export
330
-#' @docType methods
331
-#' @rdname plotPatternMarkers-methods
332
-#'
333
-#' @param object an object of type CogapsResult
334
-#' @param data the original data as a matrix
335
-#' @param patternPalette a vector indicating what color should be used
336
-#' for each pattern
337
-#' @param sampleNames names of the samples to use for labeling 
338
-#' @param samplePalette  a vector indicating what color should be used
339
-#' for each sample
340
-#' @param colDenogram logical indicating whether to display sample denogram
341
-#' @param heatmapCol pallelet giving color scheme for heatmap
342
-#' @param scale character indicating if the values should be centered and scaled
343
-#' in either the row direction or the column direction, or none. The
344
-#' default is "row". 
345
-#' @param ... additional graphical parameters to be passed to \code{heatmap.2}
346
-#' @return heatmap of the \code{data} values for the \code{patternMarkers}
347
-#' @seealso  \code{\link{heatmap.2}}
348
-#' @importFrom gplots bluered
349
-#' @importFrom gplots heatmap.2
350
-#' @importFrom stats hclust as.dist cor
351
-plotPatternMarkers <- function(object, data, patternPalette, sampleNames,
352
-samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
353
-{
354
-    if (is.null(samplePalette))
355
-        samplePalette <- rep("black", ncol(data))
356
-
357
-    ### coloring of genes
358
-    patternCols <- rep(NA, length(unlist(patternMarkers)))
359
-    names(patternCols) <- unlist(patternMarkers)
360
-    for (i in 1:length(patternMarkers))
316
+    else if (threshold == "all") # only the markers with the lowest scores
361 317
     {
362
-        patternCols[patternMarkers[[i]]] <- patternPalette[i]
318
+        patternsByMarker <- colnames(markerScores)[apply(markerScores, 1, which.min)]
319
+        markersByPattern <- sapply(colnames(markerScores), USE.NAMES=TRUE, simplify=FALSE,
320
+            function(pattern) rownames(markerScores)[which(patternsByMarker==pattern)])
363 321
     }
364
-
365
-    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),
366
-        symkey=FALSE,
367
-        symm=FALSE, 
368
-        scale=scale,
369
-        col=heatmapCol,
370
-        distfun=function(x) as.dist((1-cor(t(x)))/2),
371
-        hclustfun=function(x) stats::hclust(x,method="average"),
372
-        Rowv=FALSE,
373
-        Colv=colDenogram,
374
-        trace='none',
375
-        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),
376
-        labCol= sampleNames,
377
-        cexCol=0.8,
378
-        ColSideColors=as.character(samplePalette),
379
-        rowsep=cumsum(sapply(patternMarkers,length))
380
-    )
381
-}
322
+    return(list(
323
+        "PatternMarkers"=markersByPattern,
324
+        "PatternMarkerRanks"=markerRanks,
325
+        "PatternMarkerScores"=markerScores
326
+    ))
327
+})
382 328
 
383 329
 #' @rdname calcCoGAPSStat-methods
384 330
 #' @aliases calcCoGAPSStat
... ...
@@ -394,20 +340,18 @@ function(object, sets, whichMatrix, numPerm, ...)
394 340
         stop("sets must be a list of either measurements or samples")
395 341
     ## do a permutation test
396 342
     zMatrix <- calcZ(object, whichMatrix)
397
-    pvalUpReg <- sapply(sets,
398
-        function(thisSet)
343
+    pvalUpReg <- sapply(sets, function(thisSet)
344
+    {
345
+        lessThanCount <- rep(0, ncol(zMatrix))
346
+        actualZScore <- colMeans(zMatrix[rownames(zMatrix) %in% thisSet,])
347
+        for (n in 1:numPerm)
399 348
         {
400
-            numValues <- length(thisSet)
401
-            lessThanCount <- rep(0, ncol(zMatrix))
402
-            actualZScore <- colMeans(zMatrix[rownames(zMatrix) %in% thisSet,])
403
-            for (n in 1:numPerm)
404
-            {
405
-                permutedIndices <- sample(1:nrow(zMatrix), size=numValues, replace=FALSE)
406
-                permutedZScore <- colMeans(zMatrix[permutedIndices,])
407
-                lessThanCount <- lessThanCount + (actualZScore < permutedZScore)
408
-            }
409
-            return(lessThanCount / numPerm)
410
-        })
349
+            permutedIndices <- sample(1:nrow(zMatrix), size=length(thisSet), replace=FALSE)
350
+            permutedZScore <- colMeans(zMatrix[permutedIndices,])
351
+            lessThanCount <- lessThanCount + (actualZScore < permutedZScore)
352
+        }
353
+        return(lessThanCount / numPerm)
354
+    })
411 355
     ## calculate other quantities of interest, return a list
412 356
     pvalDownReg <- 1 - pvalUpReg # this is techinically >=, but == should rarely happen
413 357
     activityEstimate <- 1 - 2 * pvalUpReg
... ...
@@ -482,4 +426,57 @@ function(object, GStoGenes, numPerm, Pw, PwNull)
482 426
     return(finalStats)
483 427
 })
484 428
 
429
+#' heatmap of original data clustered by pattern markers statistic
430
+#' @export
431
+#' @docType methods
432
+#' @rdname plotPatternMarkers-methods
433
+#'
434
+#' @param object an object of type CogapsResult
435
+#' @param data the original data as a matrix
436
+#' @param patternPalette a vector indicating what color should be used
437
+#' for each pattern
438
+#' @param sampleNames names of the samples to use for labeling 
439
+#' @param samplePalette  a vector indicating what color should be used
440
+#' for each sample
441
+#' @param colDenogram logical indicating whether to display sample denogram
442
+#' @param heatmapCol pallelet giving color scheme for heatmap
443
+#' @param scale character indicating if the values should be centered and scaled
444
+#' in either the row direction or the column direction, or none. The
445
+#' default is "row". 
446
+#' @param ... additional graphical parameters to be passed to \code{heatmap.2}
447
+#' @return heatmap of the \code{data} values for the \code{patternMarkers}
448
+#' @seealso  \code{\link{heatmap.2}}
449
+#' @importFrom gplots bluered
450
+#' @importFrom gplots heatmap.2
451
+#' @importFrom stats hclust as.dist cor
452
+plotPatternMarkers <- function(object, data, patternPalette, sampleNames,
453
+samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
454
+{
455
+    if (is.null(samplePalette))
456
+        samplePalette <- rep("black", ncol(data))
457
+
458
+    ### coloring of genes
459
+    patternCols <- rep(NA, length(unlist(patternMarkers)))
460
+    names(patternCols) <- unlist(patternMarkers)
461
+    for (i in 1:length(patternMarkers))
462
+    {
463
+        patternCols[patternMarkers[[i]]] <- patternPalette[i]
464
+    }
485 465
 
466
+    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),
467
+        symkey=FALSE,
468
+        symm=FALSE, 
469
+        scale=scale,
470
+        col=heatmapCol,
471
+        distfun=function(x) as.dist((1-cor(t(x)))/2),
472
+        hclustfun=function(x) stats::hclust(x,method="average"),
473
+        Rowv=FALSE,
474
+        Colv=colDenogram,
475
+        trace='none',
476
+        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),
477
+        labCol= sampleNames,
478
+        cexCol=0.8,
479
+        ColSideColors=as.character(samplePalette),
480
+        rowsep=cumsum(sapply(patternMarkers,length))
481
+    )
482
+}
Browse code

Snapshots working for both phases

Tom Sherman authored on 05/09/2019 00:20:30
Showing1 changed files
... ...
@@ -18,6 +18,34 @@ createCogapsResult <- function(returnList, allParams)
18 18
         diagnostics = append(returnList$diagnostics,
19 19
             list("params"=allParams$gaps, "version"=utils::packageVersion("CoGAPS")))
20 20
     )
21
+
22
+    # label snapshots
23
+    if (allParams$nSnapshots > 0)
24
+    {
25
+        freq <- floor(allParams$gaps@nIterations / allParams$nSnapshots)
26
+        labels <- seq(freq, allParams$gaps@nIterations, freq)
27
+        if (allParams$snapshotPhase == 'equilibration' | allParams$snapshotPhase == 'all')
28
+        {
29
+            for (n in 1:allParams$nSnapshots)
30
+            {
31
+                dimnames(res@metadata$equilibrationSnapshotsA[[n]]) <- dimnames(res@featureLoadings)
32
+                dimnames(res@metadata$equilibrationSnapshotsP[[n]]) <- dimnames(res@sampleFactors)
33
+            }
34
+            names(res@metadata$equilibrationSnapshotsA) <- labels
35
+            names(res@metadata$equilibrationSnapshotsP) <- labels
36
+        }
37
+        if (allParams$snapshotPhase == 'sampling' | allParams$snapshotPhase == 'all')
38
+        {
39
+            for (n in 1:allParams$nSnapshots)
40
+            {
41
+                dimnames(res@metadata$samplingSnapshotsA[[n]]) <- dimnames(res@featureLoadings)
42
+                dimnames(res@metadata$samplingSnapshotsP[[n]]) <- dimnames(res@sampleFactors)
43
+            }
44
+            names(res@metadata$samplingSnapshotsA) <- labels
45
+            names(res@metadata$samplingSnapshotsP) <- labels
46
+        }
47
+    }        
48
+    
21 49
     validObject(res)
22 50
     return(res)
23 51
 }
... ...
@@ -232,16 +260,16 @@ function(object, data, uncertainty)
232 260
 
233 261
 patternMarkersCalculation <- function(mat, threshold, lp)
234 262
 {
263
+    if (!(threshold %in% c("cut", "all")))
264
+        stop("invalid thresold argument")
265
+    if (!is.na(lp) & length(lp) != ncol(mat))
266
+        stop("lp length must equal the number of patterns")
267
+    
235 268
     nMeasurements <- nrow(mat)
236 269
     nPatterns <- ncol(mat)
237 270
     rowMax <- t(apply(mat, 1, function(x) x / max(x)))
238
-    
239 271
     if (!is.na(lp))
240 272
     {
241
-        if (length(lp) != nPatterns)
242
-        {
243
-            warning("lp length must equal the number of columns of the Amatrix")
244
-        }
245 273
         sstat <- apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp)))
246 274
         ssranks <- rank(sstat)
247 275
         ssgenes.th <- names(sort(sstat,decreasing=FALSE,na.last=TRUE))
... ...
@@ -276,10 +304,6 @@ patternMarkersCalculation <- function(mat, threshold, lp)
276 304
             ssgenes.th <- lapply(1:max(sort(unique(pIndx))), function(x)
277 305
                 ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])
278 306
         }
279
-        else
280
-        {
281
-            stop("invalid thresold argument")
282
-        }
283 307
     }
284 308
     return(list("PatternMarkers"=ssgenes.th, "PatternRanks"=ssranks,
285 309
         "PatternMarkerScores"=sstat))       
Browse code

Change zeros in StdDev Matrix to a warning

Tom Sherman authored on 30/08/2019 17:40:16
Showing1 changed files
... ...
@@ -151,17 +151,22 @@ function(object)
151 151
 setMethod("calcZ", signature(object="CogapsResult"),
152 152
 function(object, whichMatrix)
153 153
 {
154
-
155 154
     if (whichMatrix=="featureLoadings")
156 155
     {
157 156
         if (sum(object@featureStdDev==0) > 0)
158
-            stop("zeros detected in the standard deviation matrix")
157
+        {
158
+            warning("zeros detected in the standard deviation matrix")
159
+            object@featureStdDev[object@featureStdDev==0] <- 1e-6
160
+        }
159 161
         return(object@featureLoadings / object@featureStdDev)
160 162
     }
161 163
     else if (whichMatrix=="sampleFactors")
162 164
     {
163 165
         if (sum(object@sampleStdDev==0) > 0)
164
-            stop("zeros detected in the standard deviation matrix")
166
+        {
167
+            warning("zeros detected in the standard deviation matrix")
168
+            object@sampleStdDev[object@sampleStdDev==0] <- 1e-6
169
+        }
165 170
         return(object@sampleFactors / object@sampleStdDev)
166 171
     }
167 172
     else
Browse code

Set minimum p-value in permutation test

Tom Sherman authored on 30/08/2019 17:38:53
Showing1 changed files
... ...
@@ -383,7 +383,7 @@ function(object, sets, whichMatrix, numPerm, ...)
383 383
     pvalDownReg <- 1 - pvalUpReg # this is techinically >=, but == should rarely happen
384 384
     activityEstimate <- 1 - 2 * pvalUpReg
385 385
     return(list(
386
-        'twoSidedPValue'=pmin(pvalDownReg, pvalUpReg),
386
+        'twoSidedPValue'=pmax(pmin(pvalDownReg, pvalUpReg), 1 / numPerm),
387 387
         'GSUpreg'=pvalUpReg,
388 388
         'GSDownreg'=pvalDownReg,
389 389
         'GSActEst'=activityEstimate
Browse code

Fix bogus error check

Tom Sherman authored on 30/08/2019 17:34:26
Showing1 changed files
... ...
@@ -151,12 +151,23 @@ function(object)
151 151
 setMethod("calcZ", signature(object="CogapsResult"),
152 152
 function(object, whichMatrix)
153 153
 {
154
+
154 155
     if (whichMatrix=="featureLoadings")
156
+    {
157
+        if (sum(object@featureStdDev==0) > 0)
158
+            stop("zeros detected in the standard deviation matrix")
155 159
         return(object@featureLoadings / object@featureStdDev)
160
+    }
156 161
     else if (whichMatrix=="sampleFactors")
162
+    {
163
+        if (sum(object@sampleStdDev==0) > 0)
164
+            stop("zeros detected in the standard deviation matrix")
157 165
         return(object@sampleFactors / object@sampleStdDev)
166
+    }
158 167
     else
168
+    {
159 169
         stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")
170
+    }
160 171
 })
161 172
 
162 173
 #' @rdname reconstructGene-methods
... ...
@@ -352,8 +363,6 @@ function(object, sets, whichMatrix, numPerm, ...)
352 363
     ## check input arguments
353 364
     if (!is(sets, "list"))
354 365
         stop("sets must be a list of either measurements or samples")
355
-    if (sum(object@featureStdDev==0) > 0)
356
-        stop("zeros detected in the standard deviation matrix")
357 366
     ## do a permutation test
358 367
     zMatrix <- calcZ(object, whichMatrix)
359 368
     pvalUpReg <- sapply(sets,
Browse code

Enable calcCoGAPSStat for P Matrix

Added the option to pass either 'featureLoadings' or 'sampleFactors' to the calcCoGAPSStat function so that the statistic can be used for cell type analysis as well as gene set analysis.

sherman5 authored on 29/08/2019 23:34:55
Showing1 changed files
... ...
@@ -149,20 +149,14 @@ function(object)
149 149
 #' @rdname calcZ-methods
150 150
 #' @aliases calcZ
151 151
 setMethod("calcZ", signature(object="CogapsResult"),
152
-function(object, which)
152
+function(object, whichMatrix)
153 153
 {
154
-    if (which == "feature")
155
-    {
156
-        object@featureLoadings / object@featureStdDev
157
-    }
158
-    else if (which == "sample")
159
-    {
160
-        object@sampleFactors / object@sampleStdDev
161
-    }
154
+    if (whichMatrix=="featureLoadings")
155
+        return(object@featureLoadings / object@featureStdDev)
156
+    else if (whichMatrix=="sampleFactors")
157
+        return(object@sampleFactors / object@sampleStdDev)
162 158
     else
163
-    {
164
-        stop("\"which\" must be either \"feature\" or \"sample\"")
165
-    }
159
+        stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")
166 160
 })
167 161
 
168 162
 #' @rdname reconstructGene-methods
... ...
@@ -348,118 +342,43 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
348 342
 
349 343
 #' @rdname calcCoGAPSStat-methods
350 344
 #' @aliases calcCoGAPSStat
351
-#' @importFrom stats runif
352 345
 #' @importFrom methods is
353 346
 setMethod("calcCoGAPSStat", signature(object="CogapsResult"),
354
-function(object, GStoGenes, numPerm)
347
+function(object, sets, whichMatrix, numPerm, ...)
355 348
 {
356
-    # test for std dev of zero, possible mostly in simple simulations
349
+    ## do this for backwards compatibility
350
+    if (!is.null(list(...)$GStoGenes))
351
+        sets <- GStoGenes
352
+    ## check input arguments
353
+    if (!is(sets, "list"))
354
+        stop("sets must be a list of either measurements or samples")
357 355
     if (sum(object@featureStdDev==0) > 0)
358
-    {
359
-        object@featureStdDev[object@featureStdDev==0] <- 1e-6
360
-    }
361
-
362
-    # calculate Z scores
363
-    zMatrix <- calcZ(object)
364
-
365
-    # check input arguments
366
-    if (!is(GStoGenes, "data.frame") && !is(GStoGenes, "list") && !is(GStoGenes,"GSA.genesets"))
367
-    {
368
-        stop("GStoGenes must be a data.frame,GSA.genesets, or list with format specified in the users manual.")
369
-    }
370
-
371
-    if (is(GStoGenes, "GSA.genesets"))
372
-    {
373
-        names(GStoGenes$genesets) <- GStoGenes$geneset.names
374
-        GStoGenes <- GStoGenes$genesets
375
-    }
376
-
377
-    if (is(GStoGenes, "list"))
378
-    {
379
-        GStoGenesList <- GStoGenes
380
-    }
381
-    else
382
-    {
383
-        GStoGenesList <- list()
384
-        for (i in 1:dim(GStoGenes)[2])
356
+        stop("zeros detected in the standard deviation matrix")
357
+    ## do a permutation test
358
+    zMatrix <- calcZ(object, whichMatrix)
359
+    pvalUpReg <- sapply(sets,
360
+        function(thisSet)
385 361
         {
386
-            GStoGenesList[[as.character(colnames(GStoGenes)[i])]] <- as.character(unique(GStoGenes[,i]))
387
-        }
388
-    }
389
-
390
-    # get dimensions
391
-    numGS   <- length(names(GStoGenesList))
392
-    numPatt <- dim(zMatrix)[2]
393
-    numG    <- dim(zMatrix)[1]+0.9999
394
-
395
-    # initialize matrices
396
-    statsUp   <- matrix(ncol = numGS, nrow = numPatt)
397
-    statsDown <- matrix(ncol = numGS, nrow = numPatt)
398
-    actEst    <- matrix(ncol = numGS, nrow = numPatt)
399
-    results   <- list()
400
-    zPerm     <- matrix(ncol=numPerm,nrow=numPatt)
401
-
402
-    # do permutation test
403
-    for (gs in 1:numGS)
404
-    {
405
-        genes <- GStoGenesList[[names(GStoGenesList)[gs]]]
406
-        index <- rownames(zMatrix) %in% genes
407
-        zValues <- zMatrix[index,1]
408
-        numGenes <- length(zValues)
409
-        label <- as.character(numGenes)
410
-
411
-        if (!any(names(results)==label))
412
-        {
413
-            for (p in 1:numPatt)
362
+            numValues <- length(thisSet)
363
+            lessThanCount <- rep(0, ncol(zMatrix))
364
+            actualZScore <- colMeans(zMatrix[rownames(zMatrix) %in% thisSet,])
365
+            for (n in 1:numPerm)
414 366
             {
415
-                for (j in 1:numPerm)
416
-                {
417
-                    temp <- floor(runif(numGenes,1,numG))
418
-                    temp2 <- zMatrix[temp,p]
419
-                    zPerm[p,j] <- mean(temp2)
420
-                }
367
+                permutedIndices <- sample(1:nrow(zMatrix), size=numValues, replace=FALSE)
368
+                permutedZScore <- colMeans(zMatrix[permutedIndices,])
369
+                lessThanCount <- lessThanCount + (actualZScore < permutedZScore)
421 370
             }
422
-            results[[label]] <- zPerm
423
-        }
424
-    }
425
-
426
-    # get p-value
427
-    for (p in 1:numPatt)
428
-    {
429
-        for (gs in 1:numGS)
430
-        {
431
-            genes <- GStoGenesList[[names(GStoGenesList)[gs]]]
432
-            index <- rownames(zMatrix) %in% genes
433
-            zValues <- zMatrix[index,p]
434
-            zScore <- mean(zValues)
435
-
436
-            numGenes <- length(zValues)
437
-            label <- as.character(numGenes)
438
-
439
-            permzValues <- results[[label]][p,]
440
-            ordering <- order(permzValues)
441
-            permzValues <- permzValues[ordering]
442
-
443
-            statUpReg <- 1 - sum(zScore > permzValues) / length(permzValues)
444
-            statsUp[p,gs] <- max(statUpReg, 1 / numPerm)
445
-
446
-            statDownReg <- 1 - sum(zScore < permzValues) / length(permzValues)
447
-            statsDown[p,gs] <- max(statDownReg, 1 / numPerm)
448
-
449
-            activity <- 1 - 2*max(statUpReg, 1/numPerm)
450
-            actEst[p,gs] <- activity
451
-        }
452
-    }
453
-
454
-    # format output
455
-    rownames(statsUp) <- rownames(statsDown) <- rownames(actEst) <- colnames(zMatrix)
456
-    colnames(statsUp) <- colnames(statsDown) <- colnames(actEst) <- names(GStoGenesList)
457
-
458
-    results[['GSUpreg']] <- statsUp
459
-    results[['GSDownreg']] <- statsDown
460
-    results[['GSActEst']] <- actEst
461
-
462
-    return(results)
371
+            return(lessThanCount / numPerm)
372
+        })
373
+    ## calculate other quantities of interest, return a list
374
+    pvalDownReg <- 1 - pvalUpReg # this is techinically >=, but == should rarely happen
375
+    activityEstimate <- 1 - 2 * pvalUpReg
376
+    return(list(
377
+        'twoSidedPValue'=pmin(pvalDownReg, pvalUpReg),
378
+        'GSUpreg'=pvalUpReg,
379
+        'GSDownreg'=pvalDownReg,
380
+        'GSActEst'=activityEstimate
381
+    ))
463 382
 })
464 383
 
465 384
 #' @rdname calcGeneGSStat-methods
Browse code

fixed typo

Tom Sherman authored on 18/06/2019 16:00:07
Showing1 changed files
... ...
@@ -280,7 +280,7 @@ patternMarkersCalculation <- function(mat, threshold, lp)
280 280
 setMethod("patternMarkers", signature(object="CogapsResult"),
281 281
 function(object, threshold, lp)
282 282
 {
283
-    patternMarkersCalculation(object@feataureLoadings, threshold, lp)
283
+    patternMarkersCalculation(object@featureLoadings, threshold, lp)
284 284
 })
285 285
 
286 286
 #' @rdname amplitudeMarkers-methods
Browse code

added amplitudeMarkers function

Tom Sherman authored on 18/06/2019 15:09:16
Showing1 changed files
... ...
@@ -220,66 +220,75 @@ function(object, data, uncertainty)
220 220
         lwid=c(1,10), lhei=c(1, 4, 1.2 ), main="Heatmap of Residuals")
221 221
 })
222 222
 
223
-#' @rdname patternMarkers-methods
224
-#' @aliases patternMarkers
225
-setMethod("patternMarkers", signature(object="CogapsResult"),
226
-function(object, threshold, lp)
223
+patternMarkersCalculation <- function(mat, threshold, lp)
227 224
 {
228
-    nGenes <- nrow(object@featureLoadings)
229
-    nPatterns <- ncol(object@featureLoadings)
230
-
231
-    # find the A with the highest magnitude
232
-    Arowmax <- t(apply(object@featureLoadings, 1, function(x) x/max(x)))
233
-
225
+    nMeasurements <- nrow(mat)
226
+    nPatterns <- ncol(mat)
227
+    rowMax <- t(apply(mat, 1, function(x) x / max(x)))
228
+    
234 229
     if (!is.na(lp))
235 230
     {
236 231
         if (length(lp) != nPatterns)
237 232
         {
238 233
             warning("lp length must equal the number of columns of the Amatrix")
239 234
         }
240
-        sstat <- apply(Arowmax, 1, function(x) sqrt(t(x-lp)%*%(x-lp)))
241
-        ssranks <-rank(sstat)
235
+        sstat <- apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp)))
236
+        ssranks <- rank(sstat)
242 237
         ssgenes.th <- names(sort(sstat,decreasing=FALSE,na.last=TRUE))
243 238
     }
244 239
     else
245 240
     {
246
-        # determine which genes are most associated with each pattern
247
-        sstat <- matrix(NA, nrow=nGenes, ncol=nPatterns,
248
-            dimnames=dimnames(object@featureLoadings))
249
-        ssranks <- matrix(NA, nrow=nGenes, ncol=nPatterns,
250
-            dimnames=dimnames(object@featureLoadings))
251
-        ssgenes <- matrix(NA, nrow=nGenes, ncol=nPatterns, dimnames=NULL)
241
+        # determine which measurements are most associated with each pattern
242
+        sstat <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))
243
+        ssranks <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))
244
+        ssgenes <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=NULL)
252 245
         for (i in 1:nPatterns)
253 246
         {
254 247
             lp <- rep(0,nPatterns)
255 248
             lp[i] <- 1
256
-            sstat[,i] <- unlist(apply(Arowmax, 1, function(x) sqrt(t(x-lp)%*%(x-lp))))
257
-            ssranks[,i]<-rank(sstat[,i])
258
-            ssgenes[,i]<-names(sort(sstat[,i],decreasing=FALSE,na.last=TRUE))
249
+            sstat[,i] <- unlist(apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp))))
250
+            ssranks[,i] <- rank(sstat[,i])
251
+            ssgenes[,i] <- names(sort(sstat[,i], decreasing=FALSE, na.last=TRUE))
259 252
         }
260 253
         if (threshold=="cut")
261 254
         {
262
-            geneThresh <- sapply(1:nPatterns,function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))
263
-            ssgenes.th <- sapply(1:nPatterns,function(x) ssgenes[1:geneThresh[x],x])
255
+            geneThresh <- sapply(1:nPatterns, function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))
256
+            ssgenes.th <- sapply(1:nPatterns, function(x) ssgenes[1:geneThresh[x],x])
264 257
         }
265 258
         else if (threshold=="all")
266 259
         {
267
-            pIndx<-unlist(apply(sstat,1,which.min))
260
+            pIndx <- unlist(apply(sstat, 1, which.min))
268 261
             gBYp <- list()
269
-            for(i in sort(unique(pIndx)))
262
+            for (i in sort(unique(pIndx)))
270 263
             {
271
-                gBYp[[i]]<-sapply(strsplit(names(pIndx[pIndx==i]),"[.]"),function(x) x[[1]][1])
264
+                gBYp[[i]] <- sapply(strsplit(names(pIndx[pIndx==i]),"[.]"),function(x) x[[1]][1])
272 265
             }
273 266
             ssgenes.th <- lapply(1:max(sort(unique(pIndx))), function(x)
274 267
                 ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])
275 268
         }
276 269
         else
277 270
         {
278
-            stop("Threshold arguement not viable option")
271
+            stop("invalid thresold argument")
279 272
         }
280 273
     }
281 274
     return(list("PatternMarkers"=ssgenes.th, "PatternRanks"=ssranks,
282
-        "PatternMarkerScores"=sstat))
275
+        "PatternMarkerScores"=sstat))       
276
+}
277
+
278
+#' @rdname patternMarkers-methods
279
+#' @aliases patternMarkers
280
+setMethod("patternMarkers", signature(object="CogapsResult"),
281
+function(object, threshold, lp)
282
+{
283
+    patternMarkersCalculation(object@feataureLoadings, threshold, lp)
284
+})
285
+
286
+#' @rdname amplitudeMarkers-methods
287
+#' @aliases amplitudeMarkers
288
+setMethod("amplitudeMarkers", signature(object="CogapsResult"),
289
+function(object, threshold, lp)
290
+{
291
+    patternMarkersCalculation(object@sampleFactors, threshold, lp)
283 292
 })
284 293
 
285 294
 #' heatmap of original data clustered by pattern markers statistic
Browse code

allow rds for parameters; move all critical parametrers in CogapsParams

Tom Sherman authored on 18/02/2019 19:33:02
Showing1 changed files
... ...
@@ -4,6 +4,7 @@
4 4
 #' @param returnList list from cogaps_cpp
5 5
 #' @param allParams list of all parameters
6 6
 #' @return CogapsResult object
7
+#' @importFrom utils packageVersion
7 8
 createCogapsResult <- function(returnList, allParams)
8 9
 {
9 10
     res <- new("CogapsResult",
... ...
@@ -57,6 +58,22 @@ plot.CogapsResult <- function(x, ...)
57 58
         lty=1, cex=0.8, col=colors, bty="y", ncol=5)
58 59
 }
59 60
 
61
+#' @rdname getFeatureLoadings-methods
62
+#' @aliases getFeatureLoadings
63
+setMethod("getFeatureLoadings", signature(object="CogapsResult"),
64
+function(object)
65
+{
66
+    object@featureLoadings
67
+})
68
+
69
+#' @rdname getSampleFactors-methods
70
+#' @aliases getSampleFactors
71
+setMethod("getSampleFactors", signature(object="CogapsResult"),
72
+function(object)
73
+{
74
+    object@sampleFactors
75
+})
76
+
60 77
 #' @rdname getMeanChiSq-methods
61 78
 #' @aliases getMeanChiSq
62 79
 setMethod("getMeanChiSq", signature(object="CogapsResult"),
Browse code

added fixedPatterns option for normal CoGAPS

Tom Sherman authored on 19/12/2018 14:59:41
Showing1 changed files
... ...
@@ -1,3 +1,26 @@
1
+#' convert list output from c++ code to a CogapsResult object
2
+#' @keywords internal
3
+#'
4
+#' @param returnList list from cogaps_cpp
5
+#' @param allParams list of all parameters
6
+#' @return CogapsResult object
7
+createCogapsResult <- function(returnList, allParams)
8
+{
9
+    res <- new("CogapsResult",
10
+        Amean       = returnList$Amean,
11
+        Asd         = returnList$Asd,
12
+        Pmean       = returnList$Pmean,
13
+        Psd         = returnList$Psd,
14
+        meanChiSq   = returnList$meanChiSq,
15
+        geneNames   = returnList$geneNames,
16
+        sampleNames = returnList$sampleNames,
17
+        diagnostics = append(returnList$diagnostics,
18
+            list("params"=allParams$gaps, "version"=utils::packageVersion("CoGAPS")))
19
+    )
20
+    validObject(res)
21
+    return(res)
22
+}
23
+
1 24
 setMethod("show", signature("CogapsResult"),
2 25
 function(object)
3 26
 {
Browse code

updated config to commit file permissions

Tom Sherman authored on 29/10/2018 19:56:14
Showing1 changed files
1 1
old mode 100644
2 2
new mode 100755
Browse code

added more features to GWCoGAPS and scCoGAPS

Tom Sherman authored on 08/08/2018 22:34:56
Showing1 changed files
... ...
@@ -47,7 +47,7 @@ function(object)
47 47
 setMethod("getVersion", signature(object="CogapsResult"),
48 48
 function(object)
49 49
 {
50
-    object@metadata$diagnostics$version
50
+    object@metadata$version
51 51
 })
52 52
 
53 53
 #' @rdname getOriginalParameters-methods
... ...
@@ -55,7 +55,55 @@ function(object)
55 55
 setMethod("getOriginalParameters", signature(object="CogapsResult"),
56 56
 function(object)
57 57
 {
58
-    object@metadata$diagnostics$params
58
+    object@metadata$params
59
+})
60
+
61
+#' @rdname getUnmatchedPatterns-methods
62
+#' @aliases getUnmatchedPatterns
63
+setMethod("getUnmatchedPatterns", signature(object="CogapsResult"),
64
+function(object)
65
+{
66
+    if (!is.null(object@metadata$unmatchedPatterns))
67
+        return(object@metadata$unmatchedPatterns)
68
+
69
+    message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
70
+    return(NULL)
71
+})
72
+
73
+#' @rdname getClusteredPatterns-methods
74
+#' @aliases getClusteredPatterns
75
+setMethod("getClusteredPatterns", signature(object="CogapsResult"),
76
+function(object)
77
+{
78
+    if (!is.null(object@metadata$clusteredPatterns))
79
+        return(object@metadata$clusteredPatterns)
80
+
81
+    message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
82
+    return(NULL)
83
+})
84
+
85
+#' @rdname getCorrelationToMeanPattern-methods
86
+#' @aliases getCorrelationToMeanPattern
87
+setMethod("getCorrelationToMeanPattern", signature(object="CogapsResult"),
88
+function(object)
89
+{
90
+    if (!is.null(object@metadata$CorrToMeanPattern))
91
+        return(object@metadata$CorrToMeanPattern)
92
+
93
+    message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
94
+    return(NULL)
95
+})
96
+
97
+#' @rdname getSubsets-methods
98
+#' @aliases getSubsets
99
+setMethod("getSubsets", signature(object="CogapsResult"),
100
+function(object)
101
+{
102
+    if (!is.null(object@metadata$subsets))
103
+        return(object@metadata$subsets)
104
+
105
+    message("this result was not generated with a call to GWCoGAPS or scCoGAPS")
106
+    return(NULL)
59 107
 })
60 108
 
61 109
 #' @rdname calcZ-methods
Browse code

finished up docs/vignette

Tom Sherman authored on 02/08/2018 18:25:49
Showing1 changed files
... ...
@@ -42,6 +42,22 @@ function(object)
42 42
     object@metadata$meanChiSq
43 43
 })
44 44
 
45
+#' @rdname getVersion-methods
46
+#' @aliases getVersion
47
+setMethod("getVersion", signature(object="CogapsResult"),
48
+function(object)
49
+{
50
+    object@metadata$diagnostics$version
51
+})
52
+
53
+#' @rdname getOriginalParameters-methods
54
+#' @aliases getOriginalParameters
55
+setMethod("getOriginalParameters", signature(object="CogapsResult"),
56
+function(object)
57
+{
58
+    object@metadata$diagnostics$params
59
+})
60
+
45 61
 #' @rdname calcZ-methods
46 62
 #' @aliases calcZ
47 63
 setMethod("calcZ", signature(object="CogapsResult"),
Browse code

clean up output

Tom Sherman authored on 02/08/2018 16:52:31
Showing1 changed files
... ...
@@ -21,7 +21,8 @@ plot.CogapsResult <- function(x, ...)
21 21
     xlimits <- c(0, nSamples + 1)
22 22
     ylimits <- c(0, max(x@sampleFactors) * 1.1)
23 23
 
24
-    plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude")
24
+    plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude",
25
+        xlab="Samples")
25 26
 
26 27
     for (i in 1:nFactors)
27 28
     {
... ...
@@ -77,6 +78,7 @@ function(object, genes)
77 78
 #' @aliases binaryA
78 79
 #' @importFrom gplots heatmap.2
79 80
 #' @importFrom graphics mtext
81
+#' @importFrom RColorBrewer brewer.pal
80 82
 setMethod("binaryA", signature(object="CogapsResult"),
81 83
 function(object, threshold)
82 84
 {
... ...
@@ -95,6 +97,7 @@ function(object, threshold)
95 97
 #' @aliases plotResiduals
96 98
 #' @importFrom gplots heatmap.2
97 99
 #' @importFrom grDevices colorRampPalette
100
+#' @importFrom RColorBrewer brewer.pal
98 101
 setMethod("plotResiduals", signature(object="CogapsResult"),
99 102
 function(object, data, uncertainty)
100 103
 {
... ...
@@ -142,9 +145,9 @@ function(object, threshold, lp)
142 145
         ssranks <- matrix(NA, nrow=nGenes, ncol=nPatterns,
143 146
             dimnames=dimnames(object@featureLoadings))
144 147
         ssgenes <- matrix(NA, nrow=nGenes, ncol=nPatterns, dimnames=NULL)
145
-        for (i in 1:nP)
148
+        for (i in 1:nPatterns)
146 149
         {
147
-            lp <- rep(0,dim(Amatrix)[2])
150
+            lp <- rep(0,nPatterns)
148 151
             lp[i] <- 1
149 152
             sstat[,i] <- unlist(apply(Arowmax, 1, function(x) sqrt(t(x-lp)%*%(x-lp))))
150 153
             ssranks[,i]<-rank(sstat[,i])
... ...
@@ -152,8 +155,8 @@ function(object, threshold, lp)
152 155
         }
153 156
         if (threshold=="cut")
154 157
         {
155
-            geneThresh <- sapply(1:nP,function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))
156
-            ssgenes.th <- sapply(1:nP,function(x) ssgenes[1:geneThresh[x],x])
158
+            geneThresh <- sapply(1:nPatterns,function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))
159
+            ssgenes.th <- sapply(1:nPatterns,function(x) ssgenes[1:geneThresh[x],x])
157 160
         }
158 161
         else if (threshold=="all")
159 162
         {
Browse code

vignette coming together

Tom Sherman authored on 01/08/2018 20:53:40
Showing1 changed files
... ...
@@ -15,18 +15,18 @@ function(object)
15 15
 #' @importFrom grDevices rainbow
16 16
 plot.CogapsResult <- function(x, ...)
17 17
 {
18
-    nSamples <- nrow(object@sampleFactors)
19
-    nFactors <- ncol(object@sampleFactors)
18
+    nSamples <- nrow(x@sampleFactors)
19
+    nFactors <- ncol(x@sampleFactors)
20 20
     colors <- rainbow(nFactors)
21 21
     xlimits <- c(0, nSamples + 1)
22
-    ylimits <- c(0, max(object@sampleFactors) * 1.1)
22
+    ylimits <- c(0, max(x@sampleFactors) * 1.1)
23 23
 
24 24
     plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude")
25 25
 
26 26
     for (i in 1:nFactors)
27 27
     {
28
-        lines(x=1:nSamples, y=object@sampleFactors[,i], col=colors[i])
29
-        points(x=1:nSamples, y=object@sampleFactors[,i], col=colors[i], pch=i)
28
+        lines(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i])
29
+        points(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i], pch=i)
30 30
     }
31 31
 
32 32
     legend("top", paste("Pattern", 1:nFactors, sep = ""), pch = 1:nFactors,
Browse code

moved file writers to FileParser

Tom Sherman authored on 30/07/2018 02:34:58
Showing1 changed files
... ...
@@ -11,7 +11,8 @@ function(object)
11 11
 })
12 12
 
13 13
 #' @export
14
-#' @importFrom graphics plot
14
+#' @importFrom graphics plot legend lines points
15
+#' @importFrom grDevices rainbow
15 16
 plot.CogapsResult <- function(x, ...)
16 17
 {
17 18
     nSamples <- nrow(object@sampleFactors)
... ...
@@ -75,6 +76,7 @@ function(object, genes)
75 76
 #' @rdname binaryA-methods
76 77
 #' @aliases binaryA
77 78
 #' @importFrom gplots heatmap.2
79
+#' @importFrom graphics mtext
78 80
 setMethod("binaryA", signature(object="CogapsResult"),
79 81
 function(object, threshold)
80 82
 {
... ...
@@ -92,6 +94,7 @@ function(object, threshold)
92 94
 #' @rdname plotResiduals-methods
93 95
 #' @aliases plotResiduals
94 96
 #' @importFrom gplots heatmap.2
97
+#' @importFrom grDevices colorRampPalette
95 98
 setMethod("plotResiduals", signature(object="CogapsResult"),
96 99
 function(object, data, uncertainty)
97 100
 {
... ...
@@ -194,7 +197,7 @@ function(object, threshold, lp)
194 197
 #' @seealso  \code{\link{heatmap.2}}
195 198
 #' @importFrom gplots bluered
196 199
 #' @importFrom gplots heatmap.2
197
-#' @importFrom stats hclust
200
+#' @importFrom stats hclust as.dist cor
198 201
 plotPatternMarkers <- function(object, data, patternPalette, sampleNames,
199 202
 samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
200 203
 {
... ...
@@ -229,6 +232,8 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
229 232
 
230 233
 #' @rdname calcCoGAPSStat-methods
231 234
 #' @aliases calcCoGAPSStat
235
+#' @importFrom stats runif
236
+#' @importFrom methods is
232 237
 setMethod("calcCoGAPSStat", signature(object="CogapsResult"),
233 238
 function(object, GStoGenes, numPerm)
234 239
 {
... ...
@@ -346,7 +351,7 @@ function(object, GStoGenes, numPerm)
346 351
 setMethod("calcGeneGSStat", signature(object="CogapsResult"),
347 352
 function(object, GStoGenes, numPerm, Pw, nullGenes)
348 353
 {
349
-    gsStat <- calcCoGAPSStat(object, data.frame(GSGenes), numPerm=numPerm)
354
+    gsStat <- calcCoGAPSStat(object, data.frame(GStoGenes), numPerm=numPerm)
350 355
     gsStat <- gsStat$GSUpreg
351 356
     gsStat <- -log(gsStat)
352 357
 
... ...
@@ -361,12 +366,12 @@ function(object, GStoGenes, numPerm, Pw, nullGenes)
361 366
 
362 367
     if (nullGenes)
363 368
     {
364
-        ZD <- object@featureLoadings[setdiff(row.names(object@featureLoadings), GSGenes),] /
365
-            object@featureStdDev[setdiff(row.names(object@featureLoadings), GSGenes),]
369
+        ZD <- object@featureLoadings[setdiff(row.names(object@featureLoadings), GStoGenes),] /
370
+            object@featureStdDev[setdiff(row.names(object@featureLoadings), GStoGenes),]
366 371
     }
367 372
     else
368 373
     {
369
-        ZD <- object@featureLoadings[GSGenes,]/object@featureStdDev[GSGenes,]
374
+        ZD <- object@featureLoadings[GStoGenes,]/object@featureStdDev[GStoGenes,]
370 375
     }
371 376
     outStats <- apply(sweep(ZD,2,gsStat,FUN="*"),1,sum) / (sum(gsStat))
372 377
     outStats <- outStats / apply(ZD,1,sum)
... ...
@@ -384,21 +389,21 @@ function(object, GStoGenes, numPerm, Pw, nullGenes)
384 389
 setMethod("computeGeneGSProb", signature(object="CogapsResult"),
385 390
 function(object, GStoGenes, numPerm, Pw, PwNull)
386 391
 {
387
-    geneGSStat <- calcGeneGSStat(object, Pw=Pw, GSGenes=GSGenes,
392
+    geneGSStat <- calcGeneGSStat(object, Pw=Pw, GStoGenes=GStoGenes,
388 393
         numPerm=numPerm)
389 394
 
390 395
     if (PwNull)
391 396
     {
392
-        permGSStat <- calcGeneGSStat(object, GSGenes=GSGenes, numPerm=numPerm,
397
+        permGSStat <- calcGeneGSStat(object, GStoGenes=GStoGenes, numPerm=numPerm,
393 398
             Pw=Pw, nullGenes=TRUE)
394 399
     }
395 400
     else
396 401
     {
397
-        permGSStat <- calcGeneGSStat(object, GSGenes=GSGenes, numPerm=numPerm,
402
+        permGSStat <- calcGeneGSStat(object, GStoGenes=GStoGenes, numPerm=numPerm,
398 403
             nullGenes=TRUE)
399 404
     }
400 405
 
401
-    finalStats <- sapply(GSGenes, function(x)
406
+    finalStats <- sapply(GStoGenes, function(x)
402 407
         length(which(permGSStat > geneGSStat[x])) / length(permGSStat))
403 408
 
404 409
     return(finalStats)
Tom Sherman authored on 26/07/2018 19:21:36
Showing1 changed files
... ...
@@ -225,7 +225,7 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)
225 225
         ColSideColors=as.character(samplePalette),
226 226
         rowsep=cumsum(sapply(patternMarkers,length))
227 227
     )
228
-})
228
+}
229 229
 
230 230
 #' @rdname calcCoGAPSStat-methods
231 231
 #' @aliases calcCoGAPSStat
Browse code

cleaned up directory

Tom Sherman authored on 26/07/2018 18:15:21
Showing1 changed files