Хороший способ попарного t-теста для многих строк? (Р) - PullRequest
Новая
       62

Хороший способ попарного t-теста для многих строк? (Р)

0 голосов
/ 11 июля 2020

Я пишу код, который из таблицы необработанных данных (результат анализа QIIME микробиома) генерирует t-тест для каждой группы всех строк. В каждой строке есть бактерии и значения, соответствующие каждому образцу. Таблица может быть огромной, например, 80 столбцов на 400 строк.

Group_phylum_data: 

label_Group Bacteria_Firmicutes Archaea_Other Archaea_Euryarchaeota Bacteria_Other
HC       6.771703e-05             0           0.000000000   9.480385e-04
HC       3.362588e-05             0           0.016835356   5.604313e-05
HC       0.000000e+00             0           0.000000000   2.209945e-04
EPI       0.000000e+00             0           0.001121252   2.466755e-04
EPI       0.000000e+00             0           0.000000000   3.335038e-04

Итак, теперь это только первые строки с двумя группами (H C и EPI). Я хочу провести t-тест для каждой бактерии в столбцах среди групп. Я нашел этот pairwise_t_test из пакета rstatix, и он делает именно то, что я хочу, возвращая также скорректированное значение p. Поскольку групп может быть больше двух, я выбрал этот pairwise_t_test, потому что он может обрабатывать их и выполнять статистику для каждой комбинации. 

pwc1 <- Group_phylum_data %>%
  pairwise_t_test(Bacteria_Firmicutes ~ label_Group, p.adjust.method = "bonferroni")
pwc1

Проблема в том, что я не могу найти способ сделать это al oop для ввода названия каждой бактерии и получения полной таблицы с бактериями в каждой строке и статистикой в ​​соответствующих столбцах, что-то вроде 

.y.                    group1 group2    n1    n2     p p.signif p.adj p.adj.signif
  <chr>                  <chr>  <chr>  <int> <int> <dbl> <chr>    <dbl> <chr>       
1 Bacteria_Firmicutes    EPI    HC        46    28 0.82  ns       0.82  ns          
2 Archaea_Other EPI    HC        46    28 0.453 ns       0.453 ns

, которое я получил, вручную выполнив анализ, вставив названия бактерий. Я попытался сохранить имена в массиве и заменить одно имя (в примере «Bacteria_Firmicutes») чем-то вроде names [i], но это не сработало. Может быть, это предел этого скрипта, который работает только с указанным c именем ... а может я что-то сделал не так? Или есть другой и, возможно, лучший способ получить результат, который мне нужен для этого длинного набора данных? Спасибо!

1 Ответ

0 голосов
/ 11 июля 2020

Вы можете попробовать это (Archaea_Other равно нулю, поэтому вывод не производится). Надеюсь, это помогло. 

library(reshape2)
library(rstatix)
#Melt
Melted <- reshape2::melt(data,id.vars = 'label_Group')
#Stat test
pwc1 <- Melted %>% group_by(variable) %>% 
  pairwise_t_test(value ~ label_Group, p.adjust.method = "bonferroni")

# A tibble: 3 x 10
  variable              .y.   group1 group2    n1    n2     p p.signif p.adj p.adj.signif
* <fct>                 <chr> <chr>  <chr>  <int> <int> <dbl> <chr>    <dbl> <chr>       
1 Bacteria_Firmicutes   value EPI    HC         2     3 0.273 ns       0.273 ns          
2 Archaea_Euryarchaeota value EPI    HC         2     3 0.536 ns       0.536 ns          
3 Bacteria_Other        value EPI    HC         2     3 0.761 ns       0.761 ns

Data 

data <- structure(list(label_Group = c("HC", "HC", "HC", "EPI", "EPI"
), Bacteria_Firmicutes = c(6.771703e-05, 3.362588e-05, 0, 0, 
0), Archaea_Other = c(0L, 0L, 0L, 0L, 0L), Archaea_Euryarchaeota = c(0, 
0.016835356, 0, 0.001121252, 0), Bacteria_Other = c(0.0009480385, 
5.604313e-05, 0.0002209945, 0.0002466755, 0.0003335038)), class = "data.frame", row.names = c(NA, 
-5L))
...