T-критерий на основе перестановок для задачи множественного сравнения
Я имею дело с обширными протеомными данными между двумя независимыми группами (т.е. E1 и E2) с повторами. Мои данные структурированы следующим образом:
> shark_long_AN_E1_E2
# A tibble: 11,154 × 9
MatchGroup_AUTO genename `Unique peptides` name value Extraction Tissue replicate imputation
<fct> <chr> <fct> <chr> <dbl> <fct> <fct> <fct> <chr>
1 XM_038815272.1.p1 UNANNOTATED 8 E1_AN_1_Imp 1.10 E1 AN 1 E1_AN_1_Imp
2 XM_038815272.1.p1 UNANNOTATED 8 E1_AN_2_Imp 1.96 E1 AN 2 E1_AN_2_Imp
3 XM_038815272.1.p1 UNANNOTATED 8 E1_AN_3_Imp 0.994 E1 AN 3 E1_AN_3_Imp
4 XM_038815272.1.p1 UNANNOTATED 8 E2_AN_1_Imp1 7.05 E2 AN 1 Imp1
5 XM_038815272.1.p1 UNANNOTATED 8 E2_AN_2_Imp1 7.45 E2 AN 2 Imp1
6 XM_038815272.1.p1 UNANNOTATED 8 E2_AN_3_Imp1 8.07 E2 AN 3 Imp1
Мне нужно провести t-тест между E1 и E2 (извлечение) для каждого MatchGroup_AUTO отдельно (несколько гипотез). Поскольку у меня есть тысячи уровней в MatchGroup_AUTO, я хотел бы исправить FDR, выполнив перестановку. Однако большая часть информации, которую я нашел, не учитывает несколько гипотез, и я изо всех сил пытаюсь найти решение.
На данный момент я сделал следующее:
Я рассчитал разницу среднего значения между E1 и E2 для каждой MatchGroup_AUTO на основе неперестановочных данных.
mean_diff <- shark_filtered %>%
group_by(MatchGroup_AUTO) %>%
summarise(mean_diff = mean(value[Extraction == 'E1']) - mean(value[Extraction == 'E2']))
Затем я переставляю значения:
set.seed(1979)
n <- nrow(shark_filtered)
P <- 100
variable <- shark_filtered$value
PermSamples <- matrix(0, nrow = n, ncol = P)
for (i in 1:P){
PermSamples[, i] <- sample(variable, size = n, replace = FALSE)
}
После этого шага я теряюсь, потому что каждое переставленное значение мне нужно связать с уровнем в MatchGroup_AUTO, чтобы определить средние различия для каждого уровня в MatchGroup_AUTO.
Я был бы признателен за любые советы по этому поводу.
Спасибо!
Ответы 1
Сначала вычислите t.tests по всем группам совпадений; вы можете использовать vapply.
> match_groups <- unique(d$MatchGroup_AUTO)
>
> t_obs <- t(vapply(match_groups |> setNames(match_groups), \(x)
+ unlist(t.test(value ~ Extraction, d[d$MatchGroup_AUTO == x, ])[
+ c('estimate', 'statistic', 'p.value')
+ ]),
+ FUN.VALUE=numeric(4))) |>
+ as.data.frame() |>
+ transform(dif=`estimate.mean in group E1` - `estimate.mean in group E2`) |>
+ subset(select=c(1, 2, 5, 3, 4))
Далее вы делаете по сути то же самое, с той разницей, что вы меняете метку группы, выполняя sample(Extraction) (используется значение по умолчанию replace=FALSE) в формуле и replicateP раз. Обратите внимание, что ваша вселенная перестановок ограничена, и мы можем заранее проверить это choose().
> P <- 299
> k <- length(unique(d$Extraction))
> stopifnot(choose(median(table(d$MatchGroup_AUTO)), ## for safety
+ median(table(d$MatchGroup_AUTO))/k) > P)
> t_star <- vapply(match_groups |> setNames(match_groups), \(x)
+ replicate(P, unname(t.test(value ~ sample(Extraction),
+ d[d$MatchGroup_AUTO == x, ])$statistic)),
+ FUN.VALUE=numeric(P))
Наконец, чтобы получить так называемое значение P Монте-Карло, согласно
p_mc = (Σ (I(|t*_i| ≥ |t_obs|)) + 1) / (P + 1)
мы делаем:
> p_mc <- (colSums(t(t(abs(t_star)) >= abs(t_obs[, 'statistic.t']))) + 1)/(P + 1)
дает
Множественные t-тесты с наивными значениями p, значениями MC p и в качестве бонуса скорректированными с помощью BH:
> cbind(t_obs, p_mc, p_bh=p.adjust(t_obs[, 'p.value'], 'BH'))
estimate.mean in group E1 estimate.mean in group E2 dif statistic.t p.value p_mc p_bh
1 4.597076 4.597901 -0.0008253037 -0.0009231967 0.99926953 1.00000000 0.9992695
2 4.643047 6.321108 -1.6780613879 -1.8938434438 0.06823153 0.04666667 0.4747697
3 4.540582 5.031221 -0.4906389452 -0.5561249747 0.58232250 0.56000000 0.8873486
4 5.917108 4.755134 1.1619739948 1.2391911517 0.22495015 0.23666667 0.5141718
5 4.713262 5.386373 -0.6731112674 -0.7835043642 0.43986433 0.42000000 0.8279799
6 5.379462 4.391234 0.9882280469 1.0270542115 0.31320792 0.35666667 0.6681769
7 5.932996 6.479547 -0.5465517359 -0.6270228226 0.53545149 0.55000000 0.8816326
8 4.811961 4.727288 0.0846727284 0.0817098129 0.93542079 0.93000000 0.9992695
9 6.662828 4.894009 1.7688190528 1.8998771278 0.06765757 0.08666667 0.4747697
10 6.321617 4.503693 1.8179238454 1.8503059797 0.07418277 0.10000000 0.4747697
11 6.630730 5.102590 1.5281404531 1.4996137429 0.14428645 0.18000000 0.5102062
12 6.051498 6.492977 -0.4414792845 -0.4464648353 0.65846774 0.67666667 0.9030372
13 5.657552 5.671949 -0.0143969188 -0.0162466510 0.98715067 0.99666667 0.9992695
14 5.059838 5.427452 -0.3676143223 -0.4192296000 0.67805173 0.65333333 0.9030372
15 4.659772 5.417235 -0.7574631747 -0.8456147296 0.40448328 0.37666667 0.8089666
16 5.140637 5.445461 -0.3048237813 -0.3145401554 0.75535823 0.75333333 0.9296717
17 5.296246 6.581564 -1.2853180949 -1.3553979275 0.18540987 0.17000000 0.5102062
18 4.614938 5.165485 -0.5505467871 -0.6031703703 0.55102037 0.57000000 0.8816326
19 5.881448 4.002544 1.8789041787 2.1722002594 0.03810542 0.05666667 0.4747697
20 5.289458 4.232691 1.0567672827 1.3097270311 0.20036480 0.22666667 0.5102062
21 5.924970 5.339297 0.5856728367 0.6646850732 0.51142915 0.56666667 0.8816326
22 4.425766 4.419620 0.0061465230 0.0066061789 0.99477367 0.99666667 0.9992695
23 5.665854 4.161895 1.5039585807 1.6702859383 0.10526746 0.14333333 0.5102062
24 4.517780 5.883866 -1.3660867324 -1.3156712574 0.19834572 0.20000000 0.5102062
25 4.698634 5.051978 -0.3533435981 -0.3815466113 0.70549778 0.70333333 0.9030372
26 5.080060 6.424642 -1.3445819647 -1.2889702875 0.20727129 0.24000000 0.5102062
27 6.135660 5.016718 1.1189425180 1.3489422833 0.18751982 0.18666667 0.5102062
28 5.642470 5.226990 0.4154799856 0.3985038805 0.69311626 0.61666667 0.9030372
29 5.797826 5.887192 -0.0893660131 -0.1007856340 0.92039513 0.89666667 0.9992695
30 6.521022 4.888472 1.6325504372 1.6386794699 0.11224926 0.12000000 0.5102062
31 5.132761 5.020190 0.1125713450 0.1157839337 0.90861092 0.91666667 0.9992695
32 4.954973 7.008125 -2.0531525050 -2.3425604209 0.02676272 0.04666667 0.4747697
Примечания:
Для наборов геномных данных даже
vapplyможет работать вечно. Рассмотрите возможность использования многопоточности; в Unix-подобных версиях:parallel::mclapply, в Windowsparallel::parLapplyвместо этого.Однако я не уверен, что вы действительно можете справиться с проверкой нескольких гипотез, используя такого рода перестановки, но это может быть именно тот метод, который вы ищете. Вы можете задать вопросы по статистике на сайте Cross Validated.
Данные:
set.seed(42)
mg <- 32; ex <- 2; r <- 16
d <- expand.grid(MatchGroup_AUTO=1:mg, Extraction=paste0('E', 1:ex), rep=1:r) |>
transform(value=runif (mg*ex*r, .9, 10))
Сколько у вас реплик в реальности? Если у вас есть только 3 на каждое извлечение и группу совпадений, как это выглядит, у вас может быть только вселенная
choose(3*2, 3*2/2) == 20перестановок?