I use R in bioinformatic analysis of metabolomics and gene expression data. One of my colleagues wrote a (nearly published) R program to bundle the features of many R packages for informatics, including limma and plsgenomics packages. I have come to rely on the complex program primarily for something much simpler than intended, which is extraction of significance testing from limma and plsgenomics. In the case of limma, I take a file which assigns a p-value to each feature (row). In the case of plsgenomics, I take a file which assigns a VIP score to each feature (row). When I read the manual for either package I have a difficult time finding any examples that describe even approximately what I am attempting to do. If it's not already clear, I want to generate a p-value or VIP score from either package for features in a given data matrix.

I have done some basic coding in R which I will show, the most advanced is a "for" loop to get t.test results in the same fashion I am looking to use limma and plsgenomics. If there is a way to do something very similar with limma and plsgenomics to get the results I want I would greatly appreciate the assistance in generating such a code.

t.test Code:

#Read in your file
data <- read.table ("expression-data-for-tstat-2.txt", fill=TRUE, header=TRUE, sep="\t", stringsAsFactors=FALSE)

#Assign rownames
rownames(data) <- data$TCID
data <- data [-1]

#Convert to matrix
matrix_data <- data.matrix(data)

#Group samples
colnames(matrix_data) <- c("Con", "Con", "Con", "Cd", "Cd", "Cd")

#Variable declaration
ttest_result <- NULL
tt_pval <- NULL
ttest <- NULL
wilcox_result <- NULL
wx_pval <- NULL
wxtest <- NULL
foldchange <- NULL


#Initialization
count <- (1:nrow(matrix_data))

#"for" loop: Performs t.test
for (k in count){
    index1 <- grep ("Con", colnames(matrix_data))
    index2 <- grep ("Cd", colnames(matrix_data))

    #check if columns exist
    if(length(index1) == 0) next
    if(length(index2) == 0) next

#Welch Two Sample t-test
    ttest_result <- t.test (as.numeric(matrix_data[k,][index1]), na.rm = TRUE, as.numeric(matrix_data[k,][index2]), na.rm = TRUE, paired = FALSE)
    tt_pval <- rbind (tt_pval, ttest_result$p.value)
    ttest <- rbind (ttest, ttest_result$statistic)

#Wilcoxon Test
    #wilcox_result <- wilcox.test (as.numeric(matrix_data[k,][index1]), na.rm = TRUE, as.numeric(matrix_data[k,][index2]), na.rm = TRUE, paired = FALSE)
    #wx_pval <- rbind (wx_pval, wilcox_result$p.value)
    #wxtest <- rbind (wxtest, wilcox_result$statistic)

#Fold change
    #pair_fold <- mean(as.numeric(matrix_data[k,][index2], na.rm = TRUE))/mean(as.numeric(matrix_data[k,][index1], na.rm = TRUE))
    #foldchange <- rbind (foldchange, pair_fold)
}

#Probably don't need this; puts row names on ttest, will fix with cbind function
#rownames(ttest) <- rownames(data)

#Bind object columns (so results are presented with original data)
results <- cbind (data, ttest, tt_pval)

#Write the new file
write.table(results, file="expressiondata-withtstat.txt", sep="\t", col.names=colnames(results), row.names=TRUE)

Well, at least for limma, you should be able to do something like:

y <- new("EList", list(E=matrix_data))
design <- model.matrix(~factor(colnames(matrix_data)))
fit <- lmFit(y, design)
fit <- eBayes(fit)
res <- topTable(fit, coef=2, sort.by="none", n=Inf)

Then, calling res$P.Value will give you a numeric vector containing a p-value for each row. Note that each p-value is the result of running a moderated t-test on the data from each row. This differs from the standard t-test, as moderation shares information between rows to improve the precision of the variance estimates (and thus, the power of the testing procedure).

In addition, normalization of expression values is often desirable for standard applications of the limma pipeline. For example, loess-based normalization can be performed by adding

matrix_data <- normalizeBetweenArrays(matrix_data, method="cyclicloess")

prior to the construction of the y object.

I think you are somehow assuming that limma computes a univariate test statistic that you will call in a loop. This isn't how it works. limma uses the entire matrix of expression values at once to compute improved t-tests for each row. limma computes all the p-values for you at once -- you don't need a for loop. It also collates the results for you into a data.frame with the log-fold-changes and so on.

limma can compute both empirical Bayes and ordinary t-tests for you in much less time than it would take to run a for loop.

I am amazed by your comment that you couldn't find any advice in the limma User's Guide. There are heaps of examples of how to compute p-values for two group comparison. The "quick start" code on page 13 for example tells you what you need to know.

Duplicate comment deleted