Browse code

move log_median and median_log code from affyPLM to preprocessCore

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/preprocessCore@27036 bc3139a8-67e5-0310-9ffc-ced21a209358

Ben Bolstad authored on 11/09/2007 13:53:41
Showing9 changed files

... ...
@@ -1,5 +1,5 @@
1 1
 Package: preprocessCore
2
-Version: 0.99.14
2
+Version: 0.99.15
3 3
 Title: A collection of pre-processing functions
4 4
 Author: Benjamin Milo Bolstad <bmb@bmbolstad.com>
5 5
 Maintainer: Benjamin Milo Bolstad <bmb@bmbolstad.com>
6 6
new file mode 100644
... ...
@@ -0,0 +1,11 @@
1
+#ifndef LOG_MEDIAN_H
2
+#define LOG_MEDIAN_H 1
3
+
4
+
5
+
6
+void LogMedian(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE);
7
+
8
+void logmedian(double *data, int rows, int cols, double *results, double *resultsSE);
9
+void logmedian_no_copy(double *data, int rows, int cols, double *results, double *resultsSE);
10
+
11
+#endif
0 12
new file mode 100644
... ...
@@ -0,0 +1,15 @@
1
+#ifndef MEDIAN_LOG_H
2
+#define MEDIAN_LOG_H 1
3
+
4
+
5
+
6
+void MedianLog(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE);
7
+void MedianLog_noSE(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes);
8
+
9
+void medianlog(double *data, int rows, int cols, double *results, double *resultsSE);
10
+void medianlog_no_copy(double *data, int rows, int cols, double *results, double *resultsSE);
11
+
12
+
13
+
14
+
15
+#endif
... ...
@@ -421,5 +421,87 @@ void rlm_compute_se_anova(double *Y, int y_rows,int y_cols, double *beta, double
421 421
 
422 422
 
423 423
 
424
+void MedianLog(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE){
425
+
426
+  static void(*fun)(double *, int, int, int *, double *, int, double *) = NULL;
427
+
428
+  if (fun == NULL)
429
+    fun = (void(*)(double *, int, int, int *, double *, int, double *))R_GetCCallable("preprocessCore","MedianLog");
430
+
431
+  return fun(data, rows, cols, cur_rows, results, nprobes, resultsSE);
432
+
433
+}
434
+
435
+
436
+void MedianLog_noSE(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes){
437
+
438
+  static void(*fun)(double *, int, int, int *, double *, int) = NULL;
439
+
440
+  if (fun == NULL)
441
+    fun = (void(*)(double *, int, int, int *, double *, int))R_GetCCallable("preprocessCore","MedianLog_noSE");
442
+
443
+  return fun(data, rows, cols, cur_rows, results, nprobes);
444
+
445
+}
446
+
447
+
448
+
449
+void medianlog(double *data, int rows, int cols, double *results, double *resultsSE){
450
+
451
+  static void(*fun)(double *, int, int, double *, double *) = NULL;
452
+  
453
+  if (fun == NULL)
454
+    fun = (void(*)(double *, int, int, double *, double *))R_GetCCallable("preprocessCore","medianlog");
455
+
456
+  return fun(data, rows, cols, results, resultsSE);
457
+
458
+}
459
+
460
+
461
+void medianlog_no_copy(double *data, int rows, int cols, double *results, double *resultsSE){
462
+
463
+  static void(*fun)(double *, int, int, double *, double *) = NULL;
464
+  
465
+  if (fun == NULL)
466
+    fun = (void(*)(double *, int, int, double *, double *))R_GetCCallable("preprocessCore","medianlog_no_copy");
467
+
468
+  return fun(data, rows, cols, results, resultsSE);
469
+
470
+}
471
+
472
+
473
+
474
+void LogMedian(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE){
475
+
476
+  static void(*fun)(double *, int, int, int *, double *, int, double *) = NULL;
477
+
478
+  if (fun == NULL)
479
+    fun = (void(*)(double *, int, int, int *, double *, int, double *))R_GetCCallable("preprocessCore","LogMedian");
480
+
481
+  return fun(data, rows, cols, cur_rows, results, nprobes, resultsSE);
482
+
483
+}
484
+
485
+
486
+void logmedian(double *data, int rows, int cols, double *results, double *resultsSE){
487
+
488
+  static void(*fun)(double *, int, int, double *, double *) = NULL;
489
+  
490
+  if (fun == NULL)
491
+    fun = (void(*)(double *, int, int, double *, double *))R_GetCCallable("preprocessCore","logmedian");
492
+
493
+  return fun(data, rows, cols, results, resultsSE);
494
+
495
+}
424 496
 
425 497
 
498
+void logmedian_no_copy(double *data, int rows, int cols, double *results, double *resultsSE){
499
+
500
+  static void(*fun)(double *, int, int, double *, double *) = NULL;
501
+  
502
+  if (fun == NULL)
503
+    fun = (void(*)(double *, int, int, double *, double *))R_GetCCallable("preprocessCore","logmedian_no_copy");
504
+
505
+  return fun(data, rows, cols, results, resultsSE);
506
+
507
+}
... ...
@@ -10,6 +10,7 @@
10 10
  ** May 20, 2007 - Initial version
11 11
  ** May 24-27, 2007 - add in additional registered functions
12 12
  ** Sep 9, 2007 - add the R_rlm_rma_default and R_wrlm_rma_default_model as registered functions
13
+ ** Sep 10, 2007 - add logmedian medianlog dunctions
13 14
  **
14 15
  *****************************************************/
15 16
 
... ...
@@ -17,6 +18,10 @@
17 18
 #include "medianpolish.h"
18 19
 #include "log_avg.h"
19 20
 #include "avg_log.h"
21
+
22
+#include "median_log.h"
23
+#include "log_median.h"
24
+
20 25
 #include "biweight.h"
21 26
 #include "lm.h"
22 27
 #include "rlm.h"
... ...
@@ -75,6 +80,15 @@ void R_init_preprocessCore(DllInfo *info){
75 80
   R_RegisterCCallable("preprocessCore","averagelog", (DL_FUNC)&averagelog);
76 81
   R_RegisterCCallable("preprocessCore","AverageLog_noSE", (DL_FUNC)&AverageLog_noSE);
77 82
 
83
+  R_RegisterCCallable("preprocessCore","MedianLog", (DL_FUNC)&MedianLog);
84
+  R_RegisterCCallable("preprocessCore","medianlog_no_copy", (DL_FUNC)&medianlog_no_copy);
85
+  R_RegisterCCallable("preprocessCore","medianlog", (DL_FUNC)&medianlog);
86
+  R_RegisterCCallable("preprocessCore","MedianLog_noSE", (DL_FUNC)&MedianLog_noSE);
87
+  
88
+  R_RegisterCCallable("preprocessCore","LogMedian", (DL_FUNC)&LogMedian);
89
+  R_RegisterCCallable("preprocessCore","logmedian_no_copy", (DL_FUNC)&logmedian_no_copy);
90
+  R_RegisterCCallable("preprocessCore","logmedian", (DL_FUNC)&logmedian);
91
+ 
78 92
   R_RegisterCCallable("preprocessCore","logaverage", (DL_FUNC)&logaverage);
79 93
   R_RegisterCCallable("preprocessCore","LogAverage", (DL_FUNC)&LogAverage);
80 94
 
81 95
new file mode 100644
... ...
@@ -0,0 +1,126 @@
1
+/************************************************************************
2
+ **
3
+ ** log_median.c (was previously medianPM.c)
4
+ **
5
+ ** Copyright (C) 2002-2003 Ben Bolstad
6
+ **
7
+ ** created by: B. M. Bolstad   <bolstad@stat.berkeley.edu>
8
+ ** created on: Feb 6, 2003  (but based on earlier work from Nov 2002)
9
+ **
10
+ ** last modified: Feb 6, 2003
11
+ **
12
+ ** License: GPL V2 or later (same as the rest of the Affy package)
13
+ **
14
+ ** General discussion
15
+ **
16
+ ** Implement log2 median pm summarization.
17
+ **
18
+ ** Feb 6, 2003 - Initial version of this summarization method
19
+ ** Feb 24, 2003 - remove unused int i from LogMedian()
20
+ ** Jul 23, 2003 - add a parameter for storing SE (not yet implemented)
21
+ ** Oct 10, 2003 - PLM version of threestep
22
+ ** Sep 10, 2007 - move functionality out of affyPLM (and into preprocessCore)
23
+ **
24
+ ************************************************************************/
25
+
26
+
27
+#include "log_median.h"
28
+#include "rma_common.h"
29
+
30
+#include <R.h> 
31
+#include <Rdefines.h>
32
+#include <Rmath.h>
33
+#include <Rinternals.h>
34
+
35
+#include <stdio.h>
36
+#include <stdlib.h>
37
+#include <math.h>
38
+
39
+
40
+/***************************************************************************
41
+ **
42
+ ** double LogMedian(double *x, int length)
43
+ **
44
+ ** double *x - a vector of PM intensities 
45
+ ** int length - length of *x
46
+ **
47
+ ** take the log2 of the median of PM intensities.
48
+ **
49
+ ***************************************************************************/
50
+
51
+static double log_median(double *x, int length){
52
+
53
+  double med = 0.0;
54
+  
55
+  med = median(x,length);
56
+  med = log(med)/log(2.0);
57
+
58
+  return (med);    
59
+}
60
+
61
+/***************************************************************************
62
+ **
63
+ ** double LogMedianPM(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes)
64
+ **
65
+ ** aim: given a data matrix of probe intensities, and a list of rows in the matrix 
66
+ **      corresponding to a single probeset, compute log2 Median expression measure. 
67
+ **      Note that data is a probes by chips matrix.
68
+ **
69
+ ** double *data - Probe intensity matrix
70
+ ** int rows - number of rows in matrix *data (probes)
71
+ ** int cols - number of cols in matrix *data (chips)
72
+ ** int *cur_rows - indicies of rows corresponding to current probeset
73
+ ** double *results - already allocated location to store expression measures (cols length)
74
+ ** int nprobes - number of probes in current probeset.
75
+ **
76
+ ***************************************************************************/
77
+
78
+void LogMedian(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE){
79
+  int i,j;
80
+  double *z = Calloc(nprobes*cols,double);
81
+
82
+  for (j = 0; j < cols; j++){
83
+    for (i =0; i < nprobes; i++){
84
+      z[j*nprobes + i] = data[j*rows + cur_rows[i]];  
85
+    }
86
+  } 
87
+  
88
+  for (j=0; j < cols; j++){
89
+    results[j] = log_median(&z[j*nprobes],nprobes);
90
+    resultsSE[j] = R_NaReal;
91
+  }
92
+  Free(z);
93
+}
94
+
95
+
96
+
97
+
98
+
99
+void logmedian(double *data, int rows, int cols, double *results, double *resultsSE){
100
+  int i,j;
101
+  double *buffer = Calloc(rows, double);
102
+  
103
+  for (j=0; j < cols; j++){
104
+    for (i = 0; i < rows; i++){
105
+      buffer[i] = data[j*rows + i];
106
+    }
107
+    results[j] = log_median(buffer,rows); 
108
+    resultsSE[j] = R_NaReal;
109
+  }
110
+
111
+  Free(buffer);
112
+
113
+}
114
+
115
+
116
+void logmedian_no_copy(double *data, int rows, int cols, double *results, double *resultsSE){
117
+  int i,j;
118
+  
119
+  for (j=0; j < cols; j++){
120
+    results[j] = log_median(&data[j*rows],rows); 
121
+    resultsSE[j] = R_NaReal;
122
+  }
123
+
124
+  //  Free(buffer);
125
+
126
+}
0 127
new file mode 100644
... ...
@@ -0,0 +1,11 @@
1
+#ifndef LOG_MEDIAN_H
2
+#define LOG_MEDIAN_H 1
3
+
4
+
5
+
6
+void LogMedian(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE);
7
+
8
+void logmedian(double *data, int rows, int cols, double *results, double *resultsSE);
9
+void logmedian_no_copy(double *data, int rows, int cols, double *results, double *resultsSE);
10
+
11
+#endif
0 12
new file mode 100644
... ...
@@ -0,0 +1,160 @@
1
+/************************************************************************
2
+ **
3
+ ** median_logPM.c
4
+ **
5
+ ** created by: B. M. Bolstad   <bmb@bmbolstad.com>
6
+ ** created on: Feb 6, 2003  (but based on earlier work from Nov 2002)
7
+ **
8
+ ** Copyright (C) 2003-2007   Ben Bolstad
9
+ **
10
+ ** last modified: Feb 6, 2003
11
+ **
12
+ ** License: GPL V2 or later (same as the rest of the Affy package)
13
+ **
14
+ ** General discussion
15
+ **
16
+ ** Implement median log2 pm summarization.
17
+ **
18
+ ** Feb 6, 2003 - Initial version of this summarization method
19
+ ** Feb 24, 2003 - Remove unused variable in i from MedianLog
20
+ ** Jul 23, 2003 - add SE parameter (but not yet implemented)
21
+ ** Oct 10, 2003 - added PLM version
22
+ ** Sept 9, 2007 - branch out of affyPLM into a new package preprocessCore
23
+ **
24
+ ************************************************************************/
25
+
26
+#include "rma_common.h"
27
+#include "median_log.h"
28
+
29
+#include <R.h> 
30
+#include <Rdefines.h>
31
+#include <Rmath.h>
32
+#include <Rinternals.h>
33
+
34
+#include <stdio.h>
35
+#include <stdlib.h>
36
+#include <math.h>
37
+
38
+
39
+/***************************************************************************
40
+ **
41
+ ** double MedianLog(double *x, int length)
42
+ **
43
+ ** double *x - a vector of PM intensities 
44
+ ** int length - length of *x
45
+ **
46
+ ** take the log2 of the median of PM intensities.
47
+ **
48
+ ***************************************************************************/
49
+
50
+static double median_log(double *x, int length){
51
+
52
+  double med = 0.0;
53
+  
54
+  med = median_nocopy(x,length);
55
+ 
56
+  return (med);    
57
+}
58
+
59
+/***************************************************************************
60
+ **
61
+ ** double MedianLogPM(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes)
62
+ **
63
+ ** aim: given a data matrix of probe intensities, and a list of rows in the matrix 
64
+ **      corresponding to a single probeset, compute log2 Median expression measure. 
65
+ **      Note that data is a probes by chips matrix.
66
+ **
67
+ ** double *data - Probe intensity matrix
68
+ ** int rows - number of rows in matrix *data (probes)
69
+ ** int cols - number of cols in matrix *data (chips)
70
+ ** int *cur_rows - indicies of rows corresponding to current probeset
71
+ ** double *results - already allocated location to store expression measures (cols length)
72
+ ** int nprobes - number of probes in current probeset.
73
+ **
74
+ ***************************************************************************/
75
+
76
+void MedianLog(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE){
77
+  int i,j;
78
+  double *z = Calloc(nprobes*cols,double);
79
+
80
+  for (j = 0; j < cols; j++){
81
+    for (i =0; i < nprobes; i++){
82
+      z[j*nprobes + i] = log(data[j*rows + cur_rows[i]])/log(2.0);  
83
+    }
84
+  } 
85
+  
86
+  for (j=0; j < cols; j++){
87
+    results[j] = median_log(&z[j*nprobes],nprobes); 
88
+    resultsSE[j] = R_NaReal;
89
+  }
90
+  Free(z);
91
+}
92
+
93
+
94
+
95
+/***************************************************************************
96
+ **
97
+ ** double MedianLogPM_noSE(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes)
98
+ **
99
+ ** aim: given a data matrix of probe intensities, and a list of rows in the matrix 
100
+ **      corresponding to a single probeset, compute log2 Median expression measure. 
101
+ **      Note that data is a probes by chips matrix.
102
+ **
103
+ ** double *data - Probe intensity matrix
104
+ ** int rows - number of rows in matrix *data (probes)
105
+ ** int cols - number of cols in matrix *data (chips)
106
+ ** int *cur_rows - indicies of rows corresponding to current probeset
107
+ ** double *results - already allocated location to store expression measures (cols length)
108
+ ** int nprobes - number of probes in current probeset.
109
+ **
110
+ ***************************************************************************/
111
+
112
+void MedianLog_noSE(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes){
113
+  int i,j;
114
+  double *z = Calloc(nprobes*cols,double);
115
+
116
+  for (j = 0; j < cols; j++){
117
+    for (i =0; i < nprobes; i++){
118
+      z[j*nprobes + i] = log(data[j*rows + cur_rows[i]])/log(2.0);  
119
+    }
120
+  } 
121
+  
122
+  for (j=0; j < cols; j++){
123
+    results[j] = median_log(&z[j*nprobes],nprobes);
124
+
125
+  }
126
+  Free(z);
127
+}
128
+
129
+
130
+
131
+void medianlog(double *data, int rows, int cols, double *results, double *resultsSE){
132
+  int i,j;
133
+  double *buffer = Calloc(rows, double);
134
+  
135
+  for (j=0; j < cols; j++){
136
+    for (i = 0; i < rows; i++){
137
+      buffer[i] = log2(data[j*rows + i]);
138
+    }
139
+    results[j] = median_log(buffer,rows); 
140
+    resultsSE[j] = R_NaReal;
141
+  }
142
+
143
+  Free(buffer);
144
+}
145
+
146
+
147
+
148
+
149
+void medianlog_no_copy(double *data, int rows, int cols, double *results, double *resultsSE){
150
+  int i,j;
151
+    
152
+  for (j=0; j < cols; j++){
153
+    for (i = 0; i < rows; i++){
154
+      data[j*rows + i]= log2(data[j*rows + i]);
155
+    }
156
+    results[j] = median_log(&data[j*rows],rows); 
157
+    resultsSE[j] = R_NaReal;
158
+  }
159
+
160
+}
0 161
new file mode 100644
... ...
@@ -0,0 +1,15 @@
1
+#ifndef MEDIAN_LOG_H
2
+#define MEDIAN_LOG_H 1
3
+
4
+
5
+
6
+void MedianLog(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes, double *resultsSE);
7
+void MedianLog_noSE(double *data, int rows, int cols, int *cur_rows, double *results, int nprobes);
8
+
9
+void medianlog(double *data, int rows, int cols, double *results, double *resultsSE);
10
+void medianlog_no_copy(double *data, int rows, int cols, double *results, double *resultsSE);
11
+
12
+
13
+
14
+
15
+#endif