Hi ,

I have a question about averaging microarray probes of a given dataset. It seems like avereps of limma calculates the average for replicate probes. As soon as the missing values (NA values) are introduced the calculations of average are not accurate (see example below).

Thank you.

Toufiq

> Data
    Features Col1_Counts CSC1_Counts BC_Counts
1  Feature_1    6.643856    7.228819  5.643856
2  Feature_1    3.321928    5.643856  6.228819
3  Feature_1   10.965784    7.228819  6.643856
4  Feature_4          NA    1.584963  3.321928
5  Feature_5   10.965784    5.643856  6.228819
6  Feature_6          NA          NA        NA
7  Feature_6    3.459432    3.584963  3.700440
8  Feature_8    3.906891    4.000000  4.087463
9  Feature_9    4.247928    4.321928  4.392317
10 Feature_9          NA    4.523562        NA

> avereps(Data, ID=Data$Features)

     Features    Col1_Counts CSC1_Counts BC_Counts 
[1,] "Feature_1" " 6.643856" "7.228819"  "5.643856"
[2,] "Feature_4" NA          "1.584963"  "3.321928"
[3,] "Feature_5" "10.965784" "5.643856"  "6.228819"
[4,] "Feature_6" NA          NA          NA        
[5,] "Feature_8" " 3.906891" "4.000000"  "4.087463"
[6,] "Feature_9" " 4.247928" "4.321928"  "4.392317"

avereps is working correctly. avereps works on matrices or EList objects or any object that can be coerced to be matrix. If you give it a numeric matrix, then it will compute averages over replicate probes. If you give it a character matrix, then it will simply return the first row for each unique feature ID.

Your Data object is a data.frame with a non-numeric column. When coerced to a matrix, the data.frame will become a character matrix (not a numeric matrix) and hence no averaging will take place. You can see that the output from avereps is a character matrix in this case from the quote marks around all the numbers.

The issue you are seeing has nothing to do with missing values. You simply have to give avereps a proper expression matrix. avereps will work on any data object that would work for any other limma pipeline. But you naturally can't input a data.frame to avereps and expect the function to guess which column of your matrix is supposed to be the gene IDs, just as you can't do that for any other limma functions like lmFit etc. Instead you need either to set the feature IDs to be row.names or to store the feature IDs as a gene annotation column.

avereps handles NAs automatically as can be seen:

> Data
    Features Col1_Counts CSC1_Counts BC_Counts
1  Feature_1    6.643856    7.228819  5.643856
2  Feature_1    3.321928    5.643856  6.228819
3  Feature_1   10.965784    7.228819  6.643856
4  Feature_4          NA    1.584963  3.321928
5  Feature_5   10.965784    5.643856  6.228819
6  Feature_6          NA          NA        NA
7  Feature_6    3.459432    3.584963  3.700440
8  Feature_8    3.906891    4.000000  4.087463
9  Feature_9    4.247928    4.321928  4.392317
10 Feature_9          NA    4.523562        NA
> ExprMatrix <- as.matrix(Data[,-1])
> ExprMatrix
      Col1_Counts CSC1_Counts BC_Counts
 [1,]    6.643856    7.228819  5.643856
 [2,]    3.321928    5.643856  6.228819
 [3,]   10.965784    7.228819  6.643856
 [4,]          NA    1.584963  3.321928
 [5,]   10.965784    5.643856  6.228819
 [6,]          NA          NA        NA
 [7,]    3.459432    3.584963  3.700440
 [8,]    3.906891    4.000000  4.087463
 [9,]    4.247928    4.321928  4.392317
[10,]          NA    4.523562        NA
> avereps(ExprMatrix, ID=Data[,1])
          Col1_Counts CSC1_Counts BC_Counts
Feature_1    6.977189    6.700498  6.172177
Feature_4         NaN    1.584963  3.321928
Feature_5   10.965784    5.643856  6.228819
Feature_6    3.459432    3.584963  3.700440
Feature_8    3.906891    4.000000  4.087463
Feature_9    4.247928    4.422745  4.392317