I have a gene expression data-set with log2-transformed expression values (no NAs) for 495 genes for 59 samples for which values of a continuous response variable (r) are also known (no NAs). I want to use leave-one-out cross validation to test if r of the test sample can be predicted from the sample's gene expression.

For this, I intend to use the samr R package for Significance Analysis of Microarrays to identify significant genes associated with r in the training set of samples. Then, I want to generate a linear model using the significant genes as variables, which will then be used to predict r of the test sample. I have tried the following code to begin with, but when I generate the model and examine it, I see many NAs in the model summary, which makes me suspect that I am doing something wrong.

Can someone tell me what I might be doing wrong?

Secondly, I will appreciate any comment on the use of nperms (in SAM) with a value of 100. Is it too low for an expression data-set for 495 genes.

# rVals with the r values is read as a vector from a row of a table for phenotypic data read from a tab-delimited file with sample-names as column names and phenotype features as row-names
# geneVals is the log2-transformed gene expression data-set read as a matrix from a tab-delimited file with sample-names as column names and gene-names as row-names

# Perform SAM with FDR of 5% and obtain list of significant genes

sam <- SAM(x=geneVals, y=rVals, resp.type=c("Quantitative"),
testStatistic=c("standard"), regression.method=c("standard"), logged2=TRUE, 
fdr.output=0.05, eigengene.number=1, knn.neighbors=10, nperms=100, 

sigGenes <- rbind(sam$siggenes.table$genes.up, sam$siggenes.table$genes.lo)

# Generate linear model
toModel <- data.frame(t(rbind(rVals, geneVals)), check.names=FALSE)
myModel <- lm(toModel[c('r', sigGenes[,c("Gene ID")])])

# Examine model


lm(formula = toModel[c("rVals", sigGenes[, c("Gene ID")])])

ALL 59 residuals are 0: no residual degrees of freedom!

Coefficients: (58 not defined because of singularities)
           Estimate Std. Error t value Pr(>|t|)
(Intercept)   -18.29363         NA      NA       NA
`let-7e`       -1.70545         NA      NA       NA
`miR-125a-5p`   2.43177         NA      NA       NA
`miR-151-5p`    2.67439         NA      NA       NA

I can now see that the 'NA' issue arises if during model generation the number of variables is more than the number of samples.

| cite | improve this answer | |
  • 1
    $\begingroup$ If this reasonably solves your problem you can consider accepting the answer to your own question. $\endgroup$ – Andy W Nov 3 '11 at 22:15
  • 1
    $\begingroup$ This appears to be the problem in your case. NAs can also arise when the design matrix is exactly singular (i.e. if you have exact collinearity). $\endgroup$ – Macro Apr 23 '12 at 23:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.