R: obtaining p-values from t.test for each gene

Question

I'm using the Bioconductor suite (ALL dataset) and trying to conduct a t.test for each gene. The goal is to view gene expression differences between the sexes. I can get the basic t.test down with the following:

> males <- exprs[, pData(ALL)$sex == "M"]
> females<-exprs[, pData(ALL)$sex == "F"]
> t.test(males, females)

But when I try the apply function to extract a p-value for each gene the command never ends, just keeps going in an endless loop (I think).

pvals=apply(exprs,1,function(x) {t.test(x[males],x[females])$p.value})

here is a sample of the males, there are 12625 rows (i.e., probe IDs).

> males
                                01005     01010     04006     04007     04008
1000_at                      7.597323  7.479445  7.384684  7.905312  7.065914
1001_at                      5.046194  4.932537  4.922627  4.844565  5.147762
1002_f_at                    3.900466  4.208155  4.206798  3.416923  3.945869
1003_s_at                    5.903856  6.169024  6.116890  5.687997  6.208061

Please show some sample data. Are genes individual rows or indicators in one column? Also *it falls apart* is not helpful. Post actual errors or undesired results. — Parfait, Aug 28 '18 at 21:21
Thank you. I have the expression data (in rows), this is noted as `exprs`, and then males and females. I believe my males and females are already in individual data.frames. Could I use the following: 'sapply(exprs[-1], function(x) {t.test(x[males[,1] == 1], x[females[,1] == 2])$p.value})' or should I make a few data.frame? The original data.frame already has both sexes together (I separated them). It's in `pData(ALL)$sex`. — Oars, Aug 28 '18 at 21:57
@Oars You have been amazingly resistant to any advice I have given to your previous & related questions on what to consider when doing a differential gene expression analysis. In a nutshell: 1. Don't use a t-test (or ANOVA as in your previous question) to look for differentially expressed genes. 2. You're working with microarray data where values are expression values per probe **not per gene**. You **first** need to summarise probe values per gene (usually this is done using Tukey's robust median polish), **and then** use e.g. `limma` to look for differentially expressed genes. — Maurits Evers, Aug 28 '18 at 23:33
Maurits - thanks as always. my instructor references the exprs(ALL) call as obtaining "gene expression". I'm sure you are far more accurate (seriously). For whatever reason, we wants us run a loop (using apply) that provides a p-value from a t.test for each...I don't want to say genes but he calls them genes, i.e., 445_at, 40419_at, etc. I seriously wish you were teaching the course! I wish he'd just show us the most appropriate package, why it is used, how it is applied and be done with it. I'm going through torture. — Oars, Aug 28 '18 at 23:57
@Oars Please see my post below, which should get you started. — Maurits Evers, Aug 29 '18 at 04:51

score 1 · Answer 1 · answered Aug 29 '18 at 04:47

Here is something to get you started. (At the risk of repeating myself;-) please note that this is more of a statistics/computational exercise than something you really should be doing; as explained in my comment, there exist sophisticated methods to characterise differential gene expression. A t-test (or ANOVA) is a very crude method in comparison.

We load the ALL library and data.
```
library(ALL)
data(ALL)
```
To characterise differences in mean probe intensities between male and female individuals we perform a two-sample two-sided t-test and store results in a list.
```
lst <- apply(exprs(ALL), 1, function(x)
    t.test(x[which(pData(ALL)$sex == "M")], x[which(pData(ALL)$sex == "F")]))
```

We extract the per-probe t statistic, difference in mean probe intensity and p-value, and store results in a a data.frame.

df <- do.call(rbind, lapply(lst, function(x) c(
    statistic = unname(x$statistic),
    diff = unname(diff(x$estimate)),
    pval = unname(x$p.value))))

We correct p-values for multiple hypothesis testing using the FDR method of Benjamini and Hochberg.
```
df <- transform(df, padj = p.adjust(pval, method = "BH"))
```

We inspect the first 10 rows of df sorted from smallest to largest adjusted p-value.

head(df[order(df$padj), ], n = 10)
#        statistic       diff         pval         padj
#37583_at  18.935092 -1.7717178 1.710570e-36 2.159594e-32
#38355_at  20.542586 -4.9979077 6.129942e-32 3.869526e-28
#41214_at  21.494496 -4.3233221 3.937217e-31 1.656912e-27
#34477_at  14.469711 -1.1639971 2.606867e-28 8.227924e-25
#35885_at  14.417265 -1.4006757 5.806146e-28 1.466052e-24
#38446_at -14.357159  2.3848176 1.956173e-21 4.116115e-18
#38182_at  11.052181 -0.7151076 1.140089e-19 2.056232e-16
#40097_at   9.401626 -0.5798433 8.801566e-16 1.388997e-12
#36321_at   9.208492 -0.6499951 1.823511e-15 2.557981e-12
#31534_at   8.939350 -0.5113203 1.077008e-14 1.359723e-11

We show results in a volcano plot

ggplot(df, aes(diff, -log10(padj))) +
    geom_point() +
    labs(x = "Difference in mean probe intensity", y = "Adjusted p-value")

Mauritus - I cannot thank you enough! I'll write more once I've played around with the code. — Oars, Aug 29 '18 at 13:21

score 0 · Accepted Answer · answered Sep 02 '18 at 19:20

Thanks to Maurits I was able to use his code to answer my question. I also developed the following snipped that accomplished the task (I actually like Maurits's solution better but here is another way to complete the task:

> exprs<-exprs(ALL)
> pval<-numeric()
> p.dat<-pData(ALL)$sex
> r.sims<-nrow(exprs)
> for(gene in 1:r.sims) { 
+ gexprs<-exprs[gene,]
+ g.data<-data.frame(gexprs,p.dat)
+ ttest<-t.test(gexprs[p.dat=="M"],gexprs[p.dat=="F"])
+ pval[gene]<-ttest$p.value
+ }

score 0 · Answer 3 · answered Mar 16 '19 at 12:42

0

If you are allowed to use an external package then:

library(matrixTests)
row_t_welch(exprs[, pData(ALL)$sex == "M"], exprs[, pData(ALL)$sex == "F"])

This is assuming the genes are written in rows.

answered Mar 16 '19 at 12:42

Karolis Koncevičius

9,417
9
56
89

R: obtaining p-values from t.test for each gene

3 Answers3