
Dear all,

We are currently analyzing Illumina data using the lumi and affy packages.
The results we are getting are suspicious: we get adjusted p-values in the 10^-40 to 10^-98 range.
Obviously, we are doing something wrong. The experiment consists of 20 patients of whom we have positive and negative samples, in other words, we have paired observations in each patient. We did both a paired and an unpaired analysis, but this did not matter much. We also tried SAM analysis using siggenes, but again with the same results. As a test, we also used the non-normalized data, but this did not change the results. Has anybody an idea what is going wrong?

We used the following code (abbreviated where appropriate):

library(lumi);
library(limma);
library(Matrix);
library(GOstats);
library(siggenes);
fileName <- 'Group Probe Profile TAB paired.csv';
wouter.lumi <- lumiR(fileName, lib='lumiHumanAll.db');
wouter.lumi;
summary(wouter.lumi, 'QC');

wouter.lumi.N.Q <- lumiExpresso(wouter.lumi, QC.evaluation=TRUE);
summary(wouter.lumi.N.Q);
summary(wouter.lumi.N.Q, 'QC');
## Differential expression over all data sets and gene indentification
## wouter.lumi.N <- lumiN(wouter.lumi.N.Q);
dataMatrix <- exprs(wouter.lumi.N.Q);
presentCount <- detectionCall(wouter.lumi.N.Q);
selDataMatrix <- dataMatrix[presentCount > 0,];
selProbe <- rownames(selDataMatrix);
## Paired Comparission
## specify the groups and conditions
PlaqueType <- c('A', 'B', 'A', 'B', 'A', 'B', ... 'A', 'B', 'A', 'B', 'A') ;
PatientType <- c("p01", "p01", "p02", "p02", "p03", "p03", ... "p21", "p21", "p22", "p22", "p16", "p18"); 
## top 10 genes
if(require(limma)) 
{
factor.PlaqueType <- factor(PlaqueType, levels=c("A","B"));
factor.PatientType <- factor(PatientType);
design <- model.matrix(~ 0+factor.PatientType+factor.PlaqueType);
fit <- lmFit(selDataMatrix, design);
## contrast.matrix <-makeContrasts(A-B, levels=design);
## fit2<-contrasts.fit(fit,contrast.matrix)
fit <-eBayes(fit);
if (require(lumiHumanAll.db) & require(annotate))
  {
  if (nrow(fit$genes) == 0) fit$genes <- data.frame(ID=rownames(selDataMatrix))
    geneSymbol <- getSYMBOL(fit$genes$ID, 'lumiHumanAll.db')
    fit$genes <- data.frame(fit$genes, geneSymbol=geneSymbol)
  }
## Print the top 10 genes
print(topTable(fit, number=10))

## get significant gene list with FDR adjusted p values less than 0.01
p.adj <- p.adjust(fit$p.value[,2])
sigGene.adj <- selProbe[ p.adj < 0.01]
## write.table(sigGene.adj, file="significant_genes.csv", sep=";"); 
## without FDR adjustment
sigGene <- selProbe[ fit$p.value[,2] < 0.01]
}
                      ID geneSymbol factor.PatientTypep01 factor.PatientTypep02
1820  ZdpxBOa0B_.YqnniiE   C22orf33              7.324658              7.314660
7974  oeBgiYqwKqi5.p7lTc      MLLT4              7.338124              7.305449
1678  HYFYQOpk4Y6IkqA6Hk   C1orf183              7.310553              7.350842
27    Nonl3upr2ojFAogAsU      ABCA5              7.301086              7.296429
12509 BkA6hCdSI1BR4pUQio     TAS2R3              7.311262              7.324070
11087 oTZyp7gBIx7iV4FdSg       RXRG              7.324615              7.318605
2657  oBKSggJ.exxFIq_qe4      CDY1B              7.394242              7.332206
10176 fSIikgsewouigigijo    PSMC3IP              7.372082              7.354313
12106 Wqi4jiw9KgiCoIiCLc    SPINLW1              7.339931              7.297919
2429  6lWX99edJ1RohUS3dc      CCR10              7.393743              7.383870
      factor.PatientTypep03 factor.PatientTypep04 factor.PatientTypep05
1820               7.335805              7.307950              7.314493
7974               7.344440              7.327994              7.343561
1678               7.314656              7.342585              7.320585
27                 7.275228              7.344058              7.287226
12509              7.328430              7.317047              7.325980
11087              7.362075              7.319144              7.307189
2657               7.392088              7.380562              7.326872
10176              7.418103              7.346580              7.331416
12106              7.355262              7.334333              7.312377
2429               7.402435              7.382996              7.398166

etc...

      factor.PlaqueTypeB  AveExpr        F      P.Value    adj.P.Val
1820        5.812695e-03 7.334350 79437.93 1.524058e-46 6.165472e-43
7974        9.565944e-03 7.351467 76510.25 2.253042e-46 6.165472e-43
1678        7.535491e-03 7.351725 76065.06 2.394156e-46 6.165472e-43
27          5.886626e-05 7.312479 73115.55 3.613654e-46 6.165472e-43
12509      -2.588097e-05 7.331387 72769.81 3.796430e-46 6.165472e-43
11087      -5.466541e-04 7.341570 71478.61 4.574186e-46 6.165472e-43
2657        6.219771e-03 7.359755 69690.49 5.954555e-46 6.165472e-43
10176       1.360613e-03 7.367441 69121.32 6.485202e-46 6.165472e-43
12106      -3.237106e-03 7.327304 68823.54 6.783319e-46 6.165472e-43
2429       -1.786489e-02 7.380408 68733.39 6.876512e-46 6.165472e-43

Thanks you very much in advance,

Mathijs

	[[alternative HTML version deleted]]

