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,
genenames=as.vector(rownames(geneVals)))
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
summary(myModel)
...output...
Call:
lm(formula = toModel[c("rVals", sigGenes[, c("Gene ID")])])
Residuals:
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
...
