Browse code

Breaks up discordant.R into multiple files

All exported / large functions are given their
own file, with fishersTrans and remaining
non-exported functions placed in utils.R

Max McGrath authored on 22/07/2021 17:38:58
Showing5 changed files

1 1
new file mode 100644
... ...
@@ -0,0 +1,141 @@
1
+################################################################################
2
+# File: SparCC.R
3
+# Aim : SparCC 
4
+# ------------------------------------------------------------------------------
5
+# Author: Fang Huaying (Peking University)
6
+# Email : hyfang@pku.edu.cn
7
+# Date  : 11/12/2014
8
+# ------------------------------------------------------------------------------
9
+# SparCC for counts known
10
+#   function: SparCC.count
11
+#   input:
12
+#     x ------ nxp count data matrix, row is sample, col is variable
13
+#     imax ------ resampling times from posterior distribution. default 20
14
+#     kmax ------ max iteration steps for SparCC. default is 10
15
+#     alpha ------ the threshold for strong correlation. default is 0.1
16
+#     Vmin ------ minimal variance if negative variance appears. default is 1e-4
17
+#   output: a list structure
18
+#     cov.w ------ covariance estimation
19
+#     cor.w ------ correlation estimation
20
+
21
+SparCC.count <- function(x, imax = 20, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
22
+    # dimension for w (latent variables)
23
+    p <- ncol(x);
24
+    n <- nrow(x);
25
+    # posterior distribution (alpha)
26
+    x <- x + 1;
27
+    # store generate data
28
+    y <- matrix(0, n, p);
29
+    # store covariance/correlation matrix
30
+    cov.w <- cor.w <- matrix(0, p, p);
31
+    indLow <- lower.tri(cov.w, diag = TRUE);
32
+    # store covariance/correlation for several posterior samples
33
+    covs <- cors <- matrix(0, p * (p + 1) / 2, imax);
34
+    for(i in 1:imax) {
35
+        # generate fractions from posterior distribution
36
+        y <- t(apply(x, 1, function(x) rdirichlet(n = 1, alpha = x)));
37
+        # estimate covariance/correlation
38
+        cov_cor <- SparCC.frac(x = y, kmax = kmax, alpha = alpha, Vmin = Vmin);
39
+        # store variance/correlation only low triangle 
40
+        covs[, i] <- cov_cor$cov.w[indLow];
41
+        cors[, i] <- cov_cor$cor.w[indLow];
42
+    }
43
+    # calculate median for several posterior samples
44
+    cov.w[indLow] <- apply(covs, 1, median); 
45
+    cor.w[indLow] <- apply(cors, 1, median);
46
+    #
47
+    cov.w <- cov.w + t(cov.w);
48
+    diag(cov.w) <- diag(cov.w) / 2;
49
+    cor.w <- cor.w + t(cor.w);
50
+    diag(cor.w) <- 1;
51
+    #
52
+    return(list(cov.w = cov.w, cor.w = cor.w));
53
+}
54
+
55
+#-------------------------------------------------------------------------------
56
+# SparCC for fractions known
57
+#   function: SparCC.frac
58
+#   input:
59
+#          x ------ nxp fraction data matrix, row is sample, col is variable
60
+#       kmax ------ max iteration steps for SparCC. default is 10
61
+#      alpha ------ the threshold for strong correlation. default is 0.1
62
+#       Vmin ------ minimal variance if negative variance appears. default is 1e-4
63
+#   output: a list structure
64
+#      cov.w ------ covariance estimation
65
+#      cor.w ------ correlation estimation
66
+SparCC.frac <- function(x, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
67
+    # Log transformation
68
+    x <- log(x);
69
+    p <- ncol(x);
70
+    # T0 = var(log(xi/xj)) variation matrix
71
+    TT <- stats::var(x);
72
+    T0 <- diag(TT) + rep(diag(TT), each = p) - 2 * TT;
73
+    # Variance and correlation coefficients for Basic SparCC  
74
+    rowT0 <- rowSums(T0);
75
+    var.w <- (rowT0 - sum(rowT0) / (2 * p - 2))/(p - 2);
76
+    var.w[var.w < Vmin] <- Vmin;
77
+    #cor.w <- (outer(var.w, var.w, "+") - T0 ) / 
78
+    #  sqrt(outer(var.w, var.w, "*")) / 2;
79
+    Is <- sqrt(1/var.w);
80
+    cor.w <- (var.w + rep(var.w, each = p) - T0) * Is * rep(Is, each = p) * 0.5;
81
+    # Truncated correlation in [-1, 1]
82
+    cor.w[cor.w <= - 1] <- - 1; 
83
+    cor.w[cor.w >= 1] <- 1;
84
+    # Left matrix of estimation equation
85
+    Lmat <- diag(rep(p - 2, p)) + 1; 
86
+    # Remove pairs
87
+    rp <- NULL;
88
+    # Left components
89
+    cp <- rep(TRUE, p);
90
+    # Do loops until max iteration or only 3 components left
91
+    k <- 0;  
92
+    while(k < kmax && sum(cp) > 3) {
93
+        # Left T0 = var(log(xi/xj)) after removing pairs
94
+        T02 <- T0;
95
+        # Store current correlation to find the strongest pair
96
+        curr_cor.w <- cor.w;
97
+        # Remove diagonal
98
+        diag(curr_cor.w) <- 0;
99
+        # Remove removed pairs
100
+        if(!is.null(rp)) {
101
+            curr_cor.w[rp] <- 0;
102
+        }
103
+        # Find the strongest pair in vector form
104
+        n_rp <- which.max(abs(curr_cor.w));
105
+        # Remove the pair if geater than alpha
106
+        if(abs(curr_cor.w[n_rp]) >= alpha) {
107
+            # Which pair in matrix form
108
+            t_id <- c(arrayInd(n_rp, .dim = c(p, p)));
109
+            Lmat[t_id, t_id] <- Lmat[t_id, t_id] - 1;
110
+            # Update remove pairs
111
+            n_rp <- c(n_rp, (p + 1) * sum(t_id) - 2 * p - n_rp);
112
+            rp <- c(rp, n_rp);
113
+            # Update T02
114
+            T02[rp] <- 0;
115
+            # Which component left
116
+            cp <- (diag(Lmat) > 0);
117
+            # Update variance and truncated lower by Vmin
118
+            var.w[cp] <- solve(Lmat[cp, cp], rowSums(T02[cp, cp]));
119
+            var.w[var.w <= Vmin] <- Vmin;
120
+            # Update correlation matrix and truncated by [-1, 1]
121
+            #cor.w <- (outer(var.w, var.w, "+") - T0 ) / 
122
+            #  sqrt(outer(var.w, var.w, "*")) / 2;    
123
+            Is <- sqrt(1/var.w);
124
+            cor.w <- (var.w + rep(var.w, each = p) - T0) * 
125
+                Is * rep(Is, each = p) * 0.5;
126
+            # Truncated correlation in [-1, 1]
127
+            cor.w[cor.w <= - 1] <- - 1;
128
+            cor.w[cor.w >= 1] <- 1;
129
+        }
130
+        else {
131
+            break;
132
+        }
133
+        # 
134
+        k <- k + 1;
135
+    }
136
+    # Covariance
137
+    Is <- sqrt(var.w);
138
+    cov.w <- cor.w * Is * rep(Is, each = p);
139
+    #
140
+    return(list(cov.w = cov.w, cor.w = cor.w));
141
+}
0 142
\ No newline at end of file
1 143
new file mode 100644
... ...
@@ -0,0 +1,89 @@
1
+#' @import Biobase
2
+#' @import biwt
3
+#' @import gtools
4
+#' @import MASS
5
+#' @import stats
6
+#' @import tools
7
+#' 
8
+#' @export
9
+createVectors <- function(x, y = NULL, groups, cor.method = c("spearman")) {
10
+    print(x)
11
+    if(checkInputs(x,y,groups)) {
12
+        stop("Please fix inputs.")
13
+    }
14
+    
15
+    index1 <- which(groups == 1)
16
+    index2 <- which(groups == 2)
17
+    
18
+    check <- match(cor.method, c("spearman", "pearson", "bwmc", "sparcc"))
19
+    if(is.na(check)) {
20
+        stop("Please enter spearman, pearson, bwmc or sparcc for correlation 
21
+            metric.")
22
+    }
23
+    
24
+    x <- exprs(x)
25
+    
26
+    if(is.null(y) == FALSE) {
27
+        y <- exprs(y)
28
+        data <- rbind(x, y)
29
+        data1 <- data[,index1]
30
+        data2 <- data[,index2]
31
+        featureSize = dim(x)[1]
32
+        if(cor.method == c("spearman") || cor.method == c("pearson")) {
33
+            statMatrix1 <- cor(t(data1), method = cor.method)
34
+            statMatrix2 <- cor(t(data2), method = cor.method)
35
+        }
36
+        if(cor.method == c("bwmc")) {
37
+            statMatrix1 <- biwt.cor(data1)
38
+            statMatrix2 <- biwt.cor(data2)
39
+        }
40
+        if(cor.method == c("sparcc")) {
41
+            if(min(data) < 0) {
42
+                stop("SparCC can only be applied to sequencing data.")
43
+            }
44
+            statMatrix1 <- SparCC.count(t(data1))
45
+            statMatrix2 <- SparCC.count(t(data2))
46
+        }
47
+        statMatrix1 <- statMatrix1[1:featureSize,
48
+                                   (featureSize + 1):dim(data1)[1]]
49
+        statMatrix2 <- statMatrix2[1:featureSize,
50
+                                   (featureSize + 1):dim(data1)[1]]
51
+        statVector1 <- as.vector(statMatrix1)
52
+        statVector2 <- as.vector(statMatrix2)
53
+        vector_names <- getNames(x,y)
54
+    }
55
+    
56
+    if(is.null(y)) {
57
+        data <- x # NOTE: Line added by MM to fix min(data) error below
58
+        data1 <- x[,index1]
59
+        data2 <- x[,index2]
60
+        if(cor.method == c("spearman") || cor.method == c("pearson")) {
61
+            statMatrix1 <- cor(t(data1), method = cor.method)
62
+            statMatrix2 <- cor(t(data2), method = cor.method)
63
+        }   
64
+        if(cor.method == c("bwmc")) {
65
+            statMatrix1 <- biwt.cor(data1)
66
+            statMatrix2 <- biwt.cor(data2)
67
+        }   
68
+        if(cor.method == c("sparcc")) {
69
+            if(min(data) < 0) {
70
+                stop("SparCC can only be applied to sequencing data.")
71
+            }   
72
+            # NOTE: $cor.w added to each statement below by MM as SparCC.count
73
+            #    returns two vectors, one for correlation and one for covariance
74
+            statMatrix1 <- SparCC.count(t(data1))$cor.w
75
+            statMatrix2 <- SparCC.count(t(data2))$cor.w
76
+        }   
77
+        statVector1 <- as.vector(statMatrix1)
78
+        statVector2 <- as.vector(statMatrix2)
79
+        diag <- lower.tri(statMatrix1, diag = FALSE)
80
+        statVector1 <- statMatrix1[diag]
81
+        statVector2 <- statMatrix2[diag]
82
+        vector_names <- getNames(x)
83
+    }
84
+    
85
+    names(statVector1) <- vector_names
86
+    names(statVector2) <- vector_names
87
+    
88
+    return(list(v1 = statVector1, v2 = statVector2))
89
+}
0 90
\ No newline at end of file
... ...
@@ -1,168 +1,187 @@
1 1
 # Applied from Lai et al 2007 Bioinformatics.
2
-
2
+#
3 3
 # Code written by Charlotte Siska
4
-# R functions fisherTrans, createVectors, discordantRun and makeTable
4
+# R functions fishersTrans, createVectors, discordantRun and makeTable
5 5
 # Email: charlotte.siska@ucdenver.edu
6
-
6
+#
7 7
 # Code written by Yinglei Lai
8 8
 # C code and R functions unmap and em.normal.partial.concordant
9 9
 # E-mail: ylai@gwu.edu
10
-
10
+#
11 11
 # Code written by Fang Huaying
12 12
 # R function for SparCC, adapted from https://bitbucket.org/yonatanf/sparcc
13 13
 # Email: hyfang@pku.edu.cn
14 14
 
15
-################################################################################
16
-# File: SparCC.R
17
-# Aim : SparCC 
18
-#-------------------------------------------------------------------------------
19
-# Author: Fang Huaying (Peking University)
20
-# Email : hyfang@pku.edu.cn
21
-# Date  : 11/12/2014
22
-#-------------------------------------------------------------------------------
23
-# SparCC for counts known
24
-#   function: SparCC.count
25
-#   input:
26
-#          x ------ nxp count data matrix, row is sample, col is variable
27
-#       imax ------ resampling times from posterior distribution. default 20
28
-#       kmax ------ max iteration steps for SparCC. default is 10
29
-#      alpha ------ the threshold for strong correlation. default is 0.1
30
-#       Vmin ------ minimal variance if negative variance appears. default is 1e-4
31
-#   output: a list structure
32
-#      cov.w ------ covariance estimation
33
-#      cor.w ------ correlation estimation
34
-#
35
-
36
-SparCC.count <- function(x, imax = 20, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
37
-  # dimension for w (latent variables)
38
-  p <- ncol(x);
39
-  n <- nrow(x);
40
-  # posterior distribution (alpha)
41
-  x <- x + 1;
42
-  # store generate data
43
-  y <- matrix(0, n, p);
44
-  # store covariance/correlation matrix
45
-  cov.w <- cor.w <- matrix(0, p, p);
46
-  indLow <- lower.tri(cov.w, diag = TRUE);
47
-  # store covariance/correlation for several posterior samples
48
-  covs <- cors <- matrix(0, p * (p + 1) / 2, imax);
49
-  for(i in 1:imax) {
50
-    # generate fractions from posterior distribution
51
-    y <- t(apply(x, 1, function(x) rdirichlet(n = 1, alpha = x)));
52
-    # estimate covariance/correlation
53
-    cov_cor <- SparCC.frac(x = y, kmax = kmax, alpha = alpha, Vmin = Vmin);
54
-    # store variance/correlation only low triangle 
55
-    covs[, i] <- cov_cor$cov.w[indLow];
56
-    cors[, i] <- cov_cor$cor.w[indLow];
15
+discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE, 
16
+                          subsampling = FALSE, subSize = dim(x)[1], iter = 100, 
17
+                          components = 3) {
18
+  
19
+  if(checkInputs(x,y)) {
20
+    stop("Please fix inputs.")
57 21
   }
58
-  # calculate median for several posterior samples
59
-  cov.w[indLow] <- apply(covs, 1, median); 
60
-  cor.w[indLow] <- apply(cors, 1, median);
61
-  #
62
-  cov.w <- cov.w + t(cov.w);
63
-  diag(cov.w) <- diag(cov.w) / 2;
64
-  cor.w <- cor.w + t(cor.w);
65
-  diag(cor.w) <- 1;
66
-  #
67
-  return(list(cov.w = cov.w, cor.w = cor.w));
68
-}
69
-
70
-#-------------------------------------------------------------------------------
71
-# SparCC for fractions known
72
-#   function: SparCC.frac
73
-#   input:
74
-#          x ------ nxp fraction data matrix, row is sample, col is variable
75
-#       kmax ------ max iteration steps for SparCC. default is 10
76
-#      alpha ------ the threshold for strong correlation. default is 0.1
77
-#       Vmin ------ minimal variance if negative variance appears. default is 1e-4
78
-#   output: a list structure
79
-#      cov.w ------ covariance estimation
80
-#      cor.w ------ correlation estimation
81
-SparCC.frac <- function(x, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
82
-  # Log transformation
83
-  x <- log(x);
84
-  p <- ncol(x);
85
-  # T0 = var(log(xi/xj)) variation matrix
86
-  TT <- stats::var(x);
87
-  T0 <- diag(TT) + rep(diag(TT), each = p) - 2 * TT;
88
-  # Variance and correlation coefficients for Basic SparCC  
89
-  rowT0 <- rowSums(T0);
90
-  var.w <- (rowT0 - sum(rowT0) / (2 * p - 2))/(p - 2);
91
-  var.w[var.w < Vmin] <- Vmin;
92
-  #cor.w <- (outer(var.w, var.w, "+") - T0 ) / 
93
-  #  sqrt(outer(var.w, var.w, "*")) / 2;
94
-  Is <- sqrt(1/var.w);
95
-  cor.w <- (var.w + rep(var.w, each = p) - T0) * Is * rep(Is, each = p) * 0.5;
96
-  # Truncated correlation in [-1, 1]
97
-  cor.w[cor.w <= - 1] <- - 1; 
98
-  cor.w[cor.w >= 1] <- 1;
99
-  # Left matrix of estimation equation
100
-  Lmat <- diag(rep(p - 2, p)) + 1; 
101
-  # Remove pairs
102
-  rp <- NULL;
103
-  # Left components
104
-  cp <- rep(TRUE, p);
105
-  # Do loops until max iteration or only 3 components left
106
-  k <- 0;  
107
-  while(k < kmax && sum(cp) > 3) {
108
-    # Left T0 = var(log(xi/xj)) after removing pairs
109
-    T02 <- T0;
110
-    # Store current correlation to find the strongest pair
111
-    curr_cor.w <- cor.w;
112
-    # Remove diagonal
113
-    diag(curr_cor.w) <- 0;
114
-    # Remove removed pairs
115
-    if(!is.null(rp)) {
116
-      curr_cor.w[rp] <- 0;
117
-    }
118
-    # Find the strongest pair in vector form
119
-    n_rp <- which.max(abs(curr_cor.w));
120
-    # Remove the pair if geater than alpha
121
-    if(abs(curr_cor.w[n_rp]) >= alpha) {
122
-      # Which pair in matrix form
123
-      t_id <- c(arrayInd(n_rp, .dim = c(p, p)));
124
-      Lmat[t_id, t_id] <- Lmat[t_id, t_id] - 1;
125
-      # Update remove pairs
126
-      n_rp <- c(n_rp, (p + 1) * sum(t_id) - 2 * p - n_rp);
127
-      rp <- c(rp, n_rp);
128
-      # Update T02
129
-      T02[rp] <- 0;
130
-      # Which component left
131
-      cp <- (diag(Lmat) > 0);
132
-      # Update variance and truncated lower by Vmin
133
-      var.w[cp] <- solve(Lmat[cp, cp], rowSums(T02[cp, cp]));
134
-      var.w[var.w <= Vmin] <- Vmin;
135
-      # Update correlation matrix and truncated by [-1, 1]
136
-      #cor.w <- (outer(var.w, var.w, "+") - T0 ) / 
137
-      #  sqrt(outer(var.w, var.w, "*")) / 2;    
138
-      Is <- sqrt(1/var.w);
139
-      cor.w <- (var.w + rep(var.w, each = p) - T0) * 
140
-        Is * rep(Is, each = p) * 0.5;
141
-      # Truncated correlation in [-1, 1]
142
-      cor.w[cor.w <= - 1] <- - 1;
143
-      cor.w[cor.w >= 1] <- 1;
22
+  
23
+  if(transform == TRUE) {
24
+    if(range(v1)[1] < -1 || range(v1)[2] > 1 || range(v2)[1] < -1 || 
25
+       range(v2)[2] > 1) {
26
+      stop("correlation vectors have values less than -1 and/or greater 
27
+                than 1.")
144 28
     }
145
-    else {
146
-      break;
29
+    v1 <- fishersTrans(v1)
30
+    v2 <- fishersTrans(v2)
31
+  }
32
+  
33
+  check <- match(components, c(3,5))
34
+  if(is.na(check)) {
35
+    stop("components must be 3 or 5")
36
+  }
37
+  
38
+  x <- exprs(x)
39
+  if(is.null(y) == FALSE) { y <- exprs(y) }
40
+  featureSize = dim(x)[1]
41
+  
42
+  pdata <- cbind(v1, v2)
43
+  
44
+  param1 <- sd(v1)
45
+  param2 <- sd(v2)
46
+  
47
+  if(components == 3) {
48
+    class <- cbind(rep(0, dim(pdata)[1]), rep(0, dim(pdata)[1]))
49
+    class[pdata[,1]>0+param1,1] <- 2
50
+    class[pdata[,1]<0-param1,1] <- 1
51
+    class[pdata[,2]>0+param2,2] <- 2
52
+    class[pdata[,2]<0-param2,2] <- 1
53
+    discordClass <- c(2,3,4,6,7,8)
54
+  }
55
+  
56
+  if(components == 5) {
57
+    class <- cbind(rep(0, dim(pdata)[1]), rep(0, dim(pdata)[1]))
58
+    class[pdata[,1]>0+param1,1] <- 2
59
+    class[pdata[,1]<0-param1,1] <- 1
60
+    class[pdata[,1]>0+(2*param1),1] <- 4
61
+    class[pdata[,1]<0-(2*param1),1] <- 3
62
+    class[pdata[,2]>0+param2,2] <- 2
63
+    class[pdata[,2]<0-param2,2] <- 1
64
+    class[pdata[,2]>0+(2*param2),2] <- 4
65
+    class[pdata[,2]<0-(2*param2),2] <- 3
66
+    concordClass <- c(1,7,13,19,25)
67
+    discordClass <- setdiff(1:25,concordClass)
68
+  }
69
+  
70
+  if(subsampling == TRUE) {
71
+    subSize = nrow(x)
72
+    if(is.null(y) == FALSE & nrow(y) < subSize) {
73
+      subSize = nrow(y)
74
+    }
75
+    total_mu <- rep(0,3)
76
+    total_sigma <- rep(0,3)
77
+    total_nu <- rep(0,3)
78
+    total_tau <- rep(0,3)
79
+    total_pi <- rep(0,3)
80
+    for(i in 1:iter) {
81
+      # make sure pairs are independent
82
+      rowIndex <- sample(nrow(x), subSize)
83
+      colIndex <- sample(nrow(y), subSize)
84
+      mat1 <- matrix(v1, nrow = nrow(x), byrow = FALSE)
85
+      mat2 <- matrix(v2, nrow = nrow(x), byrow = FALSE)
86
+      
87
+      subSampV1 <- sapply(1:subSize, function(x) mat1[rowIndex[x], 
88
+                                                      colIndex[x]])
89
+      subSampV2 <- sapply(1:subSize, function(x) mat2[rowIndex[x], 
90
+                                                      colIndex[x]])
91
+      
92
+      sub.pdata <- cbind(subSampV1, subSampV2)
93
+      sub.class <- cbind(rep(0, subSize), rep(0, subSize))
94
+      if(components == 3) {
95
+        sub.class[sub.pdata[,1]>0+param1,1] <- 2
96
+        sub.class[sub.pdata[,1]<0-param1,1] <- 1
97
+        sub.class[sub.pdata[,2]>0+param2,2] <- 2
98
+        sub.class[sub.pdata[,2]<0-param2,2] <- 1
99
+      }
100
+      if(components == 5) {
101
+        sub.class[sub.pdata[,1]>0+param1,1] <- 2
102
+        sub.class[sub.pdata[,1]<0-param1,1] <- 1
103
+        sub.class[sub.pdata[,1]>0+(2*param1),1] <- 4
104
+        sub.class[sub.pdata[,1]<0-(2*param1),1] <- 3
105
+        sub.class[pdata[,2]>0+param2,2] <- 2
106
+        sub.class[pdata[,2]<0-param2,2] <- 1
107
+        sub.class[pdata[,2]>0+(2*param2),2] <- 4
108
+        sub.class[pdata[,2]<0-(2*param2),2] <- 3
109
+      }
110
+      pd <- em.normal.partial.concordant(sub.pdata, sub.class, tol = 0.001, 
111
+                                         restriction = 0, 
112
+                                         constrain = c(0,-sd(pdata),sd(pdata)),
113
+                                         iteration = 1000, 
114
+                                         components = components)
115
+      total_mu <- total_mu + pd$mu_sigma[1,]
116
+      total_sigma <- total_sigma + pd$mu_sigma[2,]
117
+      total_nu <- total_nu + pd$nu_tau[1,]
118
+      total_tau <- total_tau + pd$nu_tau[2,]
119
+      total_pi <- total_pi + pd$pi
147 120
     }
148
-    # 
149
-    k <- k + 1;
121
+    
122
+    mu <- total_mu / iter
123
+    sigma <- total_sigma / iter
124
+    nu <- total_nu / iter
125
+    tau <- total_tau / iter
126
+    pi <- total_pi / iter
127
+    
128
+    finalResult <- subSampleData(pdata, class, mu, sigma, nu, tau, pi, components)
129
+    zTable <- finalResult$z
130
+    classVector <- finalResult$class
131
+  } else {
132
+    pd <- em.normal.partial.concordant(pdata, class, tol = 0.001, 
133
+                                       restriction = 0, 
134
+                                       constrain = c(0, -sd(pdata), sd(pdata)), 
135
+                                       iteration = 1000, 
136
+                                       components = components)
137
+    zTable <- pd$z
138
+    classVector <- pd$class
150 139
   }
151
-  # Covariance
152
-  Is <- sqrt(var.w);
153
-  cov.w <- cor.w * Is * rep(Is, each = p);
154
-  #
155
-  return(list(cov.w = cov.w, cor.w = cor.w));
140
+  
141
+  discordPPV <- apply(zTable, 1, function(x) sum(x[discordClass])/sum(x))
142
+  
143
+  if(is.null(y) == FALSE) {
144
+    discordPPMatrix <- matrix(discordPPV, nrow = featureSize, 
145
+                              byrow = FALSE)
146
+    classMatrix <- matrix(classVector, nrow = featureSize, byrow = FALSE)
147
+    rownames(discordPPMatrix) <- rownames(x)
148
+    colnames(discordPPMatrix) <- rownames(y)
149
+    rownames(classMatrix) <- rownames(x)
150
+    colnames(classMatrix) <- rownames(y)
151
+    
152
+    vector_names <- getNames(x,y)
153
+    names(discordPPV) <- vector_names
154
+    names(classVector) <- vector_names
155
+  }
156
+  
157
+  if(is.null(y)) {
158
+    discordPPMatrix <- matrix(NA,nrow = featureSize, ncol = featureSize)
159
+    classMatrix <- discordPPMatrix
160
+    diag <- lower.tri(discordPPMatrix, diag = FALSE)
161
+    discordPPMatrix[diag] <- discordPPV
162
+    classMatrix[diag] <- classVector
163
+    rownames(discordPPMatrix) <- rownames(x)
164
+    colnames(discordPPMatrix) <- rownames(x)
165
+    rownames(classMatrix) <- rownames(x)
166
+    colnames(classMatrix) <- rownames(x)
167
+    vector_names <- getNames(x)
168
+    names(discordPPV) <- vector_names
169
+    names(classVector) <- vector_names
170
+  }
171
+  
172
+  zTable <- t(apply(zTable, 1, function(x) x/sum(x)))
173
+  rownames(zTable) <- vector_names
174
+  
175
+  return(list(discordPPMatrix = discordPPMatrix, discordPPVector = 
176
+                discordPPV, classMatrix = classMatrix, classVector = classVector, 
177
+              probMatrix = zTable, loglik = pd$loglik))
156 178
 }
157 179
 
158
-
159
-
160 180
 #modified from package mclust
161 181
 unmap <- function(classification, components){
162 182
     n <- length(classification)
163 183
     # u <- sort(unique(classification)) # OG Code
164
-    # u <- 0:max(classification) # Max's potential fix
165
-    u <- 0:(components - 1)
184
+    u <- 0:(components - 1) # Max's potential fix
166 185
     labs <- as.character(u)
167 186
     k <- length(u)
168 187
     z <- matrix(0, n, k)
... ...
@@ -171,7 +190,6 @@ unmap <- function(classification, components){
171 190
     return(z)
172 191
 }
173 192
 
174
-
175 193
 em.normal.partial.concordant <- function(data, class, tol=0.001, 
176 194
     restriction=0, constrain=0, iteration=1000, components){
177 195
     n <- as.integer(dim(data)[1])
... ...
@@ -195,10 +213,6 @@ em.normal.partial.concordant <- function(data, class, tol=0.001,
195 213
     tau <- double(g)
196 214
     loglik <- double(1)
197 215
     convergence <- integer(1)
198
-    # results <- .C("em_normal_partial_concordant", as.double(data[,1]), 
199
-    #     as.double(data[,2]), as.double(t(zxy)), n, pi, mu, sigma, nu, tau, g, 
200
-    #     loglik, as.double(tol), as.integer(restriction), 
201
-    #     as.integer(constrain), as.integer(iteration), convergence)
202 216
     
203 217
     results <- em_normal_partial_concordant_cpp(as.double(data[,1]), 
204 218
                                                 as.double(data[,2]), 
... ...
@@ -215,375 +229,4 @@ em.normal.partial.concordant <- function(data, class, tol=0.001,
215 229
         results[[7]]), nu_tau=rbind(results[[8]], results[[9]]), 
216 230
         loglik=results[[11]], class=apply(array(results[[3]], dim=c(n,g*g)),1,
217 231
         order,decreasing=TRUE)[1,], z=array(results[[3]], dim=c(n,g*g))))
218
-}
219
-
220
-subSampleData <- function(pdata, class, mu, sigma, nu, tau, pi, components) {
221
-    n <- as.integer(dim(pdata)[1])
222
-    g <- as.integer(nlevels(as.factor(class)))
223
-
224
-    yl.outer <- function(k, zx, zy){
225
-        return( c(zx[k,] %o% zy[k,]) )
226
-    }
227
-    yl.diag <- function(k, z){
228
-        return( c(diag(z[k,])) )
229
-    }
230
-
231
-    zx <- unmap(class[,1], components = components)
232
-    zy <- unmap(class[,2], components = components)
233
-    zxy <- sapply(1:dim(zx)[1], yl.outer, zx, zy)
234
-
235
-    # results <- .C("subsampling", as.double(pdata[,1]), as.double(pdata[,2]), 
236
-    #     as.double(t(zxy)), n, as.double(pi), as.double(mu), as.double(sigma), 
237
-    #     as.double(nu), as.double(tau), g)
238
-    
239
-    results <- subsampling_cpp(as.double(pdata[,1]), as.double(pdata[,2]), 
240
-                               as.double(t(zxy)), n, as.double(pi), 
241
-                               as.double(mu), as.double(sigma), as.double(nu), 
242
-                               as.double(tau), g)
243
-    
244
-    return(list(pi=t(array(results[[5]],dim=c(g,g))), 
245
-        mu_sigma=rbind(results[[6]], results[[7]]), nu_tau=rbind(results[[8]],
246
-        results[[9]]), class=apply(array(results[[3]], dim=c(n,g*g)),1,order,
247
-        decreasing=TRUE)[1,], z=array(results[[3]], dim=c(n,g*g))))
248
-} 
249
-
250
-#' @export
251
-fishersTrans <- function(rho) {
252
-    r = (1+rho)/(1-rho)
253
-    z = 0.5*log(r,base = exp(1))
254
-    return(z)
255
-}
256
-
257
-getNames <- function(x, y = NULL) {
258
-    if(is.null(y) == FALSE) {
259
-        namesMatrix <- NULL
260
-        for(i in 1:nrow(x)) {
261
-            tempMatrix <- cbind(rep(rownames(x)[i],nrow(y)), rownames(y))
262
-            namesMatrix <- rbind(namesMatrix, tempMatrix)
263
-        }
264
-    }  
265
-
266
-    if(is.null(y)) {
267
-        temp <- matrix(NA,nrow = nrow(x), ncol = nrow(x))
268
-        diag <- lower.tri(temp, diag = FALSE)
269
-        temp[diag] <- rep(1, sum(diag == TRUE))
270
-
271
-        namesMatrix <- NULL
272
-
273
-        for(i in 1:dim(temp)[1]) {
274
-            outputCol <- temp[,i]
275
-            index <- which(is.na(outputCol) == FALSE)
276
-            if(length(index) > 0) {
277
-                tempMatrix <- cbind(rep(rownames(x)[i],length(index)), 
278
-                    rownames(x)[index])
279
-                namesMatrix <- rbind(namesMatrix, tempMatrix)
280
-            }
281
-        }
282
-    }
283
-
284
-    vector_names <- apply(namesMatrix, 1, function(k) paste(k[1],"_",k[2],
285
-        sep = ""))
286
-    return(vector_names)
287
-}
288
-
289
-checkInputs <- function(x,y,groups = NULL) {
290
-    issue = FALSE
291
-    if(is.null(groups) == FALSE  && unique(unique(groups) != c(1,2))) {
292
-        print("groups vector must consist of 1s and 2s corresponding to first
293
-             and second group.")
294
-        issue = TRUE
295
-    }
296
-    if(mode(x) != "S4" || (is.null(y) == FALSE && mode(y) != "S4")) {
297
-        print("data matrices x and/or y must be type ExpressionSet")
298
-        issue = TRUE
299
-    }
300
-    return(issue)
301
-}
302
-
303
-#' @import Biobase
304
-#' @import biwt
305
-#' @import gtools
306
-#' @import MASS
307
-#' @import stats
308
-#' @import tools
309
-#' 
310
-#' @export
311
-createVectors <- function(x, y = NULL, groups, cor.method = c("spearman")) {
312
-    print(x)
313
-    if(checkInputs(x,y,groups)) {
314
-        stop("Please fix inputs.")
315
-    }
316
-
317
-    index1 <- which(groups == 1)
318
-    index2 <- which(groups == 2)
319
-
320
-    check <- match(cor.method, c("spearman", "pearson", "bwmc", "sparcc"))
321
-    if(is.na(check)) {
322
-        stop("Please enter spearman, pearson, bwmc or sparcc for correlation 
323
-            metric.")
324
-    }
325
-
326
-    x <- exprs(x)
327
-
328
-    if(is.null(y) == FALSE) {
329
-        y <- exprs(y)
330
-        data <- rbind(x, y)
331
-        data1 <- data[,index1]
332
-        data2 <- data[,index2]
333
-        featureSize = dim(x)[1]
334
-        if(cor.method == c("spearman") || cor.method == c("pearson")) {
335
-            statMatrix1 <- cor(t(data1), method = cor.method)
336
-            statMatrix2 <- cor(t(data2), method = cor.method)
337
-        }
338
-        if(cor.method == c("bwmc")) {
339
-            statMatrix1 <- biwt.cor(data1)
340
-            statMatrix2 <- biwt.cor(data2)
341
-        }
342
-        if(cor.method == c("sparcc")) {
343
-            if(min(data) < 0) {
344
-                stop("SparCC can only be applied to sequencing data.")
345
-            }
346
-            statMatrix1 <- SparCC.count(t(data1))
347
-            statMatrix2 <- SparCC.count(t(data2))
348
-        }
349
-        statMatrix1 <- statMatrix1[1:featureSize,
350
-            (featureSize + 1):dim(data1)[1]]
351
-        statMatrix2 <- statMatrix2[1:featureSize,
352
-            (featureSize + 1):dim(data1)[1]]
353
-        statVector1 <- as.vector(statMatrix1)
354
-        statVector2 <- as.vector(statMatrix2)
355
-        vector_names <- getNames(x,y)
356
-    }
357
-
358
-    if(is.null(y)) {
359
-        data <- x # NOTE: Line added by MM to fix min(data) error below
360
-        data1 <- x[,index1]
361
-        data2 <- x[,index2]
362
-        if(cor.method == c("spearman") || cor.method == c("pearson")) {
363
-            statMatrix1 <- cor(t(data1), method = cor.method)
364
-            statMatrix2 <- cor(t(data2), method = cor.method)
365
-        }   
366
-        if(cor.method == c("bwmc")) {
367
-            statMatrix1 <- biwt.cor(data1)
368
-            statMatrix2 <- biwt.cor(data2)
369
-        }   
370
-        if(cor.method == c("sparcc")) {
371
-            if(min(data) < 0) {
372
-                stop("SparCC can only be applied to sequencing data.")
373
-            }   
374
-            # NOTE: $cor.w added to each statement below  by MM as SparCC.count
375
-            #    returns two vectors, one for correlation and one for covariance
376
-            statMatrix1 <- SparCC.count(t(data1))$cor.w
377
-            statMatrix2 <- SparCC.count(t(data2))$cor.w
378
-        }   
379
-        statVector1 <- as.vector(statMatrix1)
380
-        statVector2 <- as.vector(statMatrix2)
381
-        diag <- lower.tri(statMatrix1, diag = FALSE)
382
-        statVector1 <- statMatrix1[diag]
383
-        statVector2 <- statMatrix2[diag]
384
-        vector_names <- getNames(x)
385
-    }
386
-
387
-    names(statVector1) <- vector_names
388
-    names(statVector2) <- vector_names
389
-
390
-    return(list(v1 = statVector1, v2 = statVector2))
391
-}
392
-
393
-#' @export
394
-discordantRun <- function(v1, v2, x, y = NULL, transform = TRUE, 
395
-    subsampling = FALSE, subSize = dim(x)[1], iter = 100, components = 3) {
396
-
397
-    if(checkInputs(x,y)) {
398
-        stop("Please fix inputs.")
399
-    }
400
-
401
-    if(transform == TRUE) {
402
-        if(range(v1)[1] < -1 || range(v1)[2] > 1 || range(v2)[1] < -1 || 
403
-            range(v2)[2] > 1) {
404
-            stop("correlation vectors have values less than -1 and/or greater 
405
-                than 1.")
406
-        }
407
-        v1 <- fishersTrans(v1)
408
-        v2 <- fishersTrans(v2)
409
-    }
410
-
411
-    check <- match(components, c(3,5))
412
-    if(is.na(check)) {
413
-        stop("components must be 3 or 5")
414
-    }
415
-
416
-    x <- exprs(x)
417
-    if(is.null(y) == FALSE) { y <- exprs(y) }
418
-    featureSize = dim(x)[1]
419
-
420
-    pdata <- cbind(v1, v2)
421
-	
422
-    param1 <- sd(v1)
423
-    param2 <- sd(v2)
424
-	
425
-    if(components == 3) {
426
-        class <- cbind(rep(0, dim(pdata)[1]), rep(0, dim(pdata)[1]))
427
-        class[pdata[,1]>0+param1,1] <- 2
428
-        class[pdata[,1]<0-param1,1] <- 1
429
-        class[pdata[,2]>0+param2,2] <- 2
430
-        class[pdata[,2]<0-param2,2] <- 1
431
-        discordClass <- c(2,3,4,6,7,8)
432
-    }
433
-
434
-    if(components == 5) {
435
-        class <- cbind(rep(0, dim(pdata)[1]), rep(0, dim(pdata)[1]))
436
-        class[pdata[,1]>0+param1,1] <- 2
437
-        class[pdata[,1]<0-param1,1] <- 1
438
-        class[pdata[,1]>0+(2*param1),1] <- 4
439
-        class[pdata[,1]<0-(2*param1),1] <- 3
440
-        class[pdata[,2]>0+param2,2] <- 2
441
-        class[pdata[,2]<0-param2,2] <- 1
442
-        class[pdata[,2]>0+(2*param2),2] <- 4
443
-        class[pdata[,2]<0-(2*param2),2] <- 3
444
-        concordClass <- c(1,7,13,19,25)
445
-        discordClass <- setdiff(1:25,concordClass)
446
-    }
447
-
448
-    if(subsampling == TRUE) {
449
-        subSize = nrow(x)
450
-        if(is.null(y) == FALSE & nrow(y) < subSize) {
451
-            subSize = nrow(y)
452
-        }
453
-        total_mu <- rep(0,3)
454
-        total_sigma <- rep(0,3)
455
-        total_nu <- rep(0,3)
456
-        total_tau <- rep(0,3)
457
-        total_pi <- rep(0,3)
458
-        for(i in 1:iter) {
459
-        # make sure pairs are independent
460
-            rowIndex <- sample(nrow(x), subSize)
461
-            colIndex <- sample(nrow(y), subSize)
462
-            mat1 <- matrix(v1, nrow = nrow(x), byrow = FALSE)
463
-            mat2 <- matrix(v2, nrow = nrow(x), byrow = FALSE)
464
-
465
-            subSampV1 <- sapply(1:subSize, function(x) mat1[rowIndex[x], 
466
-                colIndex[x]])
467
-            subSampV2 <- sapply(1:subSize, function(x) mat2[rowIndex[x], 
468
-                colIndex[x]])
469
-
470
-            sub.pdata <- cbind(subSampV1, subSampV2)
471
-            sub.class <- cbind(rep(0, subSize), rep(0, subSize))
472
-            if(components == 3) {
473
-                sub.class[sub.pdata[,1]>0+param1,1] <- 2
474
-                sub.class[sub.pdata[,1]<0-param1,1] <- 1
475
-                sub.class[sub.pdata[,2]>0+param2,2] <- 2
476
-                sub.class[sub.pdata[,2]<0-param2,2] <- 1
477
-            }
478
-            if(components == 5) {
479
-                sub.class[sub.pdata[,1]>0+param1,1] <- 2
480
-                sub.class[sub.pdata[,1]<0-param1,1] <- 1
481
-                sub.class[sub.pdata[,1]>0+(2*param1),1] <- 4
482
-                sub.class[sub.pdata[,1]<0-(2*param1),1] <- 3
483
-                sub.class[pdata[,2]>0+param2,2] <- 2
484
-                sub.class[pdata[,2]<0-param2,2] <- 1
485
-                sub.class[pdata[,2]>0+(2*param2),2] <- 4
486
-                sub.class[pdata[,2]<0-(2*param2),2] <- 3
487
-            }
488
-            pd <- em.normal.partial.concordant(sub.pdata, sub.class, 
489
-                tol=0.001, restriction=0, constrain=c(0,-sd(pdata),sd(pdata)),
490
-                iteration=1000, components=components)
491
-            total_mu <- total_mu + pd$mu_sigma[1,]
492
-            total_sigma <- total_sigma + pd$mu_sigma[2,]
493
-            total_nu <- total_nu + pd$nu_tau[1,]
494
-            total_tau <- total_tau + pd$nu_tau[2,]
495
-            total_pi <- total_pi + pd$pi
496
-        }
497
-
498
-        mu <- total_mu/iter
499
-        sigma <- total_sigma/iter
500
-        nu <- total_nu/iter
501
-        tau <- total_tau/iter
502
-        pi <- total_pi/iter
503
-
504
-        finalResult <- subSampleData(pdata, class, mu, sigma, nu, tau, pi, components)
505
-        zTable <- finalResult$z
506
-        classVector <- finalResult$class
507
-    } else {
508
-        pd <- em.normal.partial.concordant(pdata, class, tol=0.001, 
509
-            restriction=0, constrain=c(0,-sd(pdata),sd(pdata)), 
510
-            iteration=1000, components = components)
511
-        zTable <- pd$z
512
-        classVector <- pd$class
513
-    }
514
-
515
-    discordPPV <- apply(zTable, 1, function(x) sum(x[discordClass])/sum(x))
516
-
517
-    if(is.null(y) == FALSE) {
518
-        discordPPMatrix <- matrix(discordPPV, nrow = featureSize, 
519
-            byrow = FALSE)
520
-        classMatrix <- matrix(classVector, nrow = featureSize, byrow = FALSE)
521
-        rownames(discordPPMatrix) <- rownames(x)
522
-        colnames(discordPPMatrix) <- rownames(y)
523
-        rownames(classMatrix) <- rownames(x)
524
-        colnames(classMatrix) <- rownames(y)
525
-
526
-        vector_names <- getNames(x,y)
527
-        names(discordPPV) <- vector_names
528
-        names(classVector) <- vector_names
529
-    }
530
-
531
-    if(is.null(y)) {
532
-        discordPPMatrix <- matrix(NA,nrow = featureSize, ncol = featureSize)
533
-        classMatrix <- discordPPMatrix
534
-        diag <- lower.tri(discordPPMatrix, diag = FALSE)
535
-        discordPPMatrix[diag] <- discordPPV
536
-        classMatrix[diag] <- classVector
537
-        rownames(discordPPMatrix) <- rownames(x)
538
-        colnames(discordPPMatrix) <- rownames(x)
539
-        rownames(classMatrix) <- rownames(x)
540
-        colnames(classMatrix) <- rownames(x)
541
-        vector_names <- getNames(x)
542
-        names(discordPPV) <- vector_names
543
-        names(classVector) <- vector_names
544
-    }
545
-
546
-    zTable <- t(apply(zTable, 1, function(x) x/sum(x)))
547
-    rownames(zTable) <- vector_names
548
-
549
-    return(list(discordPPMatrix = discordPPMatrix, discordPPVector = 
550
-        discordPPV, classMatrix = classMatrix, classVector = classVector, 
551
-        probMatrix = zTable, loglik = pd$loglik))
552
-}
553
-
554
-#' @export
555
-splitMADOutlier <- function(mat, filter0 = TRUE, threshold = 2) {
556
-    if(mode(mat) != "S4") {
557
-        stop("data matrix mat must be type ExpressionSet")
558
-    }
559
-    mat <- exprs(mat)
560
-    maxMAD <- c()
561
-    mat.filtered <- NULL
562
-    index <- c()
563
-    for(i in 1:nrow(mat)) {
564
-        y <- mat[i,]
565
-        x <- y
566
-        if(length(which(is.infinite(x))) >0 ){
567
-            x <- x[-which(is.infinite(x))]
568
-        }
569
-        if(length(which(is.na(x)))> 0) {
570
-            x <- x[-which(is.na(x))]
571
-        }
572
-        median.x <- median(x)
573
-        left <- x[x <= median.x]
574
-        right <- x[x > median.x]
575
-        left.mad <- mad(left)
576
-        right.mad <- mad(right)
577
-        leftThresh <- median(left) - left.mad*threshold
578
-        rightThresh <- median(right) + right.mad*threshold 
579
-        left.check <- sum(left < leftThresh)
580
-        right.check <- sum(right > rightThresh)
581
-        if(left.check == 0 & right.check == 0) {
582
-            mat.filtered <- rbind(mat.filtered, y)
583
-            rownames(mat.filtered)[nrow(mat.filtered)] <- rownames(mat)[i]
584
-            index <- c(index,i)
585
-        }
586
-    }
587
-    colnames(mat.filtered) = colnames(mat)
588
-    return(list(mat.filtered = mat.filtered, index = index))
589
-}
232
+}
590 233
\ No newline at end of file
591 234
new file mode 100644
... ...
@@ -0,0 +1,36 @@
1
+#' @export
2
+splitMADOutlier <- function(mat, filter0 = TRUE, threshold = 2) {
3
+    if(mode(mat) != "S4") {
4
+        stop("data matrix mat must be type ExpressionSet")
5
+    }
6
+    mat <- exprs(mat)
7
+    maxMAD <- c()
8
+    mat.filtered <- NULL
9
+    index <- c()
10
+    for(i in 1:nrow(mat)) {
11
+        y <- mat[i,]
12
+        x <- y
13
+        if(length(which(is.infinite(x))) >0 ){
14
+            x <- x[-which(is.infinite(x))]
15
+        }
16
+        if(length(which(is.na(x)))> 0) {
17
+            x <- x[-which(is.na(x))]
18
+        }
19
+        median.x <- median(x)
20
+        left <- x[x <= median.x]
21
+        right <- x[x > median.x]
22
+        left.mad <- mad(left)
23
+        right.mad <- mad(right)
24
+        leftThresh <- median(left) - left.mad*threshold
25
+        rightThresh <- median(right) + right.mad*threshold 
26
+        left.check <- sum(left < leftThresh)
27
+        right.check <- sum(right > rightThresh)
28
+        if(left.check == 0 & right.check == 0) {
29
+            mat.filtered <- rbind(mat.filtered, y)
30
+            rownames(mat.filtered)[nrow(mat.filtered)] <- rownames(mat)[i]
31
+            index <- c(index,i)
32
+        }
33
+    }
34
+    colnames(mat.filtered) = colnames(mat)
35
+    return(list(mat.filtered = mat.filtered, index = index))
36
+}
0 37
\ No newline at end of file
1 38
new file mode 100644
... ...
@@ -0,0 +1,78 @@
1
+#' @export
2
+fishersTrans <- function(rho) {
3
+    r = (1 + rho) / (1 - rho)
4
+    z = 0.5 * log(r, base = exp(1))
5
+    return(z)
6
+}
7
+
8
+subSampleData <- function(pdata, class, mu, sigma, nu, tau, pi, components) {
9
+    n <- as.integer(dim(pdata)[1])
10
+    g <- as.integer(nlevels(as.factor(class)))
11
+    
12
+    yl.outer <- function(k, zx, zy){
13
+        return( c(zx[k,] %o% zy[k,]) )
14
+    }
15
+    yl.diag <- function(k, z){
16
+        return( c(diag(z[k,])) )
17
+    }
18
+    
19
+    zx <- unmap(class[,1], components = components)
20
+    zy <- unmap(class[,2], components = components)
21
+    zxy <- sapply(1:dim(zx)[1], yl.outer, zx, zy)
22
+    
23
+    results <- subsampling_cpp(as.double(pdata[,1]), as.double(pdata[,2]), 
24
+                               as.double(t(zxy)), n, as.double(pi), 
25
+                               as.double(mu), as.double(sigma), as.double(nu), 
26
+                               as.double(tau), g)
27
+    
28
+    return(list(pi = t(array(results[[5]], dim=c(g,g))), 
29
+                mu_sigma = rbind(results[[6]], results[[7]]), 
30
+                nu_tau = rbind(results[[8]], results[[9]]), 
31
+                class = apply(array(results[[3]], dim = c(n, g*g)), 1, order,
32
+                              decreasing=TRUE)[1,], 
33
+                z = array(results[[3]], dim=c(n,g*g))))
34
+} 
35
+
36
+getNames <- function(x, y = NULL) {
37
+    if(is.null(y) == FALSE) {
38
+        namesMatrix <- NULL
39
+        for(i in 1:nrow(x)) {
40
+            tempMatrix <- cbind(rep(rownames(x)[i], nrow(y)), rownames(y))
41
+            namesMatrix <- rbind(namesMatrix, tempMatrix)
42
+        }
43
+    } else {
44
+        temp <- matrix(NA,nrow = nrow(x), ncol = nrow(x))
45
+        diag <- lower.tri(temp, diag = FALSE)
46
+        temp[diag] <- rep(1, sum(diag == TRUE))
47
+        
48
+        namesMatrix <- NULL
49
+        
50
+        for(i in 1:dim(temp)[1]) {
51
+            outputCol <- temp[,i]
52
+            index <- which(is.na(outputCol) == FALSE)
53
+            if(length(index) > 0) {
54
+                tempMatrix <- cbind(rep(rownames(x)[i],length(index)), 
55
+                                    rownames(x)[index])
56
+                namesMatrix <- rbind(namesMatrix, tempMatrix)
57
+            }
58
+        }
59
+    }
60
+    
61
+    vector_names <- apply(namesMatrix, 1, function(k) paste(k[1],"_",k[2],
62
+                                                            sep = ""))
63
+    return(vector_names)
64
+}
65
+
66
+checkInputs <- function(x, y, groups = NULL) {
67
+    issue = FALSE
68
+    if(is.null(groups) == FALSE  && unique(unique(groups) != c(1,2))) {
69
+        print("groups vector must consist of 1s and 2s corresponding to first
70
+             and second group.")
71
+        issue = TRUE
72
+    }
73
+    if(mode(x) != "S4" || (is.null(y) == FALSE && mode(y) != "S4")) {
74
+        print("data matrices x and/or y must be type ExpressionSet")
75
+        issue = TRUE
76
+    }
77
+    return(issue)
78
+}
0 79
\ No newline at end of file