R: получение значений p из t.test для каждого гена

Question

Я использую набор Bioconductor (набор данных ALL) и пытаюсь провести t.test для каждого гена. Цель состоит в том, чтобы увидеть различия в экспрессии генов между полами. Я могу получить базовый t.test со следующим:

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

Но когда я пытаюсь применить функцию apply для извлечения p-значения для каждого гена, команда никогда не заканчивается, просто продолжает бесконечный цикл (я думаю).

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

вот пример самцов, 12625 строк (т. Е. Идентификаторы зондов).

> 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

Answer 1

Благодаря Маурицу я смог использовать его код для ответа на мой вопрос. Я также разработал следующий фрагмент кода, который выполнил задачу (мне действительно нравится решение Морица, но есть еще один способ выполнить задачу:

> 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
+ }

Answer 2

Если вам разрешено использовать внешний пакет, тогда:

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

Предполагается, что гены записаны в строках.

Answer 3

Вот кое-что, с чего можно начать. (С риском повторения ;-), пожалуйста, обратите внимание, что это больше статистическое / вычислительное упражнение, чем то, что вы действительно должны делать; как объясняется в моем комментарии, существуют сложные методы для характеристики дифференциальной экспрессии генов. T-тест (или ANOVA) является очень грубым методом для сравнения.

1. Загружаем ВСЕ библиотеку и данные.

   library(ALL)
   data(ALL)
   

2. Чтобы охарактеризовать различия в средней интенсивности проб между мужчинами и женщинами, мы выполняем двухсторонний t-критерий и сохраняем результаты в list.

   lst <- apply(exprs(ALL), 1, function(x)
       t.test(x[which(pData(ALL)$sex == "M")], x[which(pData(ALL)$sex == "F")]))
   

3. Мы извлекаем t-статистику для каждого зонда, разницу в средней интенсивности зонда и значение p и сохраняем результаты в 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))))
   

4. Мы корректируем p-значения для проверки нескольких гипотез, используя метод FDR Бенджамини и Хохберга .

   df <- transform(df, padj = p.adjust(pval, method = "BH"))
   

5. Мы проверяем первые 10 строк df, отсортированных от наименьшего к наибольшему скорректированного значения p.

   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
   

6. Мы показываем результаты на участке вулкана

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