Отредактировано: Как написать цикл для теста PostHoc Tukey и множественного сравнения

Question

Это мой образец матрицы больших данных, и каждый столбец назван с несколькими данными и разделен подчеркиванием.

structure(list(Gene = c("AGI4120.1_UBQ", "AGI570.1_Acin"), WT_Tissue_0T_1 = c(0.886461437, 1.093164915), WT_Tissue_0T_2 = c(1.075140682, 1.229862834), WT_Tissue_0T_3 = c(0.632903012, 1.094003128), WT_Tissue_1T_1 = c(0.883151274, 1.26322126), WT_Tissue_1T_2 = c(1.005627276, 0.962729188), WT_Tissue_1T_3 = c(0.87123469, 0.968078993), WT_Tissue_3T_1 = c(0.723601456, 0.633890322), WT_Tissue_3T_2 = c(0.392585237, 0.534819363), WT_Tissue_3T_3 = c(0.640185369, 1.021934772), WT_Tissue_5T_1 = c(0.720291294, 0.589244505), WT_Tissue_5T_2 = c(0.362131744, 0.475251717), WT_Tissue_5T_3 = c(0.549486925, 0.618177919), mut1_Tissue_0T_1 = c(1.464415756, 1.130533457), mut1_Tissue_0T_2 = c(1.01489573, 1.114915728), mut1_Tissue_0T_3 = c(1.171797418, 1.399956009), mut1_Tissue_1T_1 = c(0.927507448, 1.231911575), mut1_Tissue_1T_2 = c(1.089705396, 1.256782289 ), mut1_Tissue_1T_3 = c(0.993048659, 0.999044465), mut1_Tissue_3T_1 = c(1.000993049, 1.103486794), mut1_Tissue_3T_2 = c(1.062562066, 0.883617224 ), mut1_Tissue_3T_3 = c(1.037404833, 0.851875438), mut1_Tissue_5T_1 = c(0.730883813, 0.437440083), mut1_Tissue_5T_2 = c(0.480635551, 0.298762126 ), mut1_Tissue_5T_3 = c(0.85468388, 0.614923997)), row.names = c(NA, -2L), class = c("tbl_df", "tbl", "data.frame"), spec = structure(list( cols = list(Gene = structure(list(), class = c("collector_character", "collector")), WT_Tissue_0T_1 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_0T_2 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_0T_3 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_1T_1 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_1T_2 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_1T_3 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_3T_1 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_3T_2 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_3T_3 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_5T_1 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_5T_2 = structure(list(), class = c("collector_double", "collector")), WT_Tissue_5T_3 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_0T_1 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_0T_2 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_0T_3 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_1T_1 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_1T_2 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_1T_3 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_3T_1 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_3T_2 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_3T_3 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_5T_1 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_5T_2 = structure(list(), class = c("collector_double", "collector")), mut1_Tissue_5T_3 = structure(list(), class = c("collector_double", "collector"))), default = structure(list(), class = c("collector_guess", "collector"))), class = "col_spec"))

Я хочу следовать тесту Тьюки и построить гистограммы для каждого гена(Ответ против времени; заполнен обоими генотипами) с несколькими буквами сравнения.

Синтаксис

df1 <- df %>% 
         gather(var, response, WT_Tissue_0T_1:mut1_Tissue_5T_3) %>% 
         separate(var, c("Genotype", "Tissue", "Time"), sep = "_") %>% 
         arrange(desc(Gene))
df2 <- df1 %>% 
         group_by(`Gene`,Genotype,Tissue,Time) %>% 
         mutate(Response=mean(response),n=n(),se=sd(response)/sqrt(n))

Двусторонняя ANOVA

fit1 <- aov(Response ~ Genotype*Time, df2) 

В дальнейшем я хочучтобы перейти к тесту Тьюки (множественное сравнение), например, для гена "AGI4120.1_UBQ", построить график зависимости отклика от времени и посмотреть, как ведет себя каждый генотип (WT & mut1) в каждой временной точке (0T, 1T, 3T и 5T)?если ответ существенно отличается или нет и обозначается буквами на графике.

Как показано ниже, синтаксис lsmeans объединяет все гены в один и дает вывод, как я могу сделать так, чтобы он зацикливался для каждого гена отдельно (т. Е. "AGI4120.1_UBQ", "AGI570.1_Acin") и получал буквы для отображениястатистически различающиеся группы (также называемые «компактным отображением букв»)

 lsmeans(fit1, pairwise ~ Genotype | Time) 

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

df2$genotype <- factor(df2$genotype, levels = c("WT","mut1")) colours <- c('#336600','#ffcc00')library(ggplot2)ggplot(df2,aes( x=Time, y=Response, fill=Genotype))+geom_bar(stat='identity', position='dodge')+scale_fill_manual(values=colours)+geom_errorbar(aes(ymin=average_measure-se, ymax=average_measure+se)+facet_wrap(~`Gene`)+labs(x='Time', y='Response')

ОжидаетсяГрафик для обозначения компактного буквенного дисплея

Буду признателен за помощь, если это возможно.

Answer 1

Существует ряд проблем с вашим кодом.Я бы даже сказал, что это не совсем подходящая статья для StackOverflow, так как проблемы здесь разнообразны и не могут быть расширены за пределы этой специфической совокупности ошибок и синтаксических проблем.Но я отправлю несколько предложений в качестве ответа, чтобы вы начали - надеюсь, они помогут.

1.lsmeans():
Эта функция ожидает подобранную модель (например, из lm()) или объект ref.grid.Но вы передаете ему фрейм данных, без каких-либо свойств регрессии, необходимых для вычисления наименьших квадратов.(Подумайте: как lsmeans() узнает, какой должна быть зависимая переменная, когда вы запрашиваете парное сравнение с Genes в качестве независимой переменной?) Проверьте Использование lsmeans документации для получения более подробной информации.

Исходя из ваших данных, я бы сказал, что вы, вероятно, хотите запустить многоуровневую регрессию, используя что-то вроде пакета lme4, с вложенными Gene и Genotype и, возможно, Timeгруппировка уровней.

Но для демонстрации я сделаю это просто с lm().Передача подогнанного регрессионного объекта в lsmeans() работает как задумано:

fit <- lm(Response ~ Gene + Genotype + Time, data=df2)

lsmeans(fit, pairwise ~ Gene)[[2]] 

# Output
contrast                        estimate        SE df t.ratio p.value
AGI4120.1_UBQ - AGI570.1_Acin -0.0515123 0.0299492 42   -1.72  0.0928

Results are averaged over the levels of: Genotype, Time 

2.ggplot():
Вы не определили colours или average_measure в предоставленном вами коде;вызов этих необъявленных переменных вызовет сбой.

Конструктивно я бы предложил вам использовать df1 и разрешить ggplot выполнить группировку, а не group_by в df2.Затем используйте stat="summary" и fun.y="mean", чтобы выполнить вычисления summarise(), которые вы делали в df2.Это также позволяет вам использовать функцию mean_se для ваших панелей ошибок.Например:

ggplot(df1,aes( x=Time, y=response, fill=Genotype))+ 
  geom_bar(stat='summary', fun.y='mean', position=position_dodge(0.9))+
    stat_summary(fun.data = mean_se, geom = "errorbar", 
                 color="gray60", width=.1, position=position_dodge(0.9)) +
  scale_fill_manual(values=c("steelblue","orange"))+
  facet_wrap(~`Gene`)+ 
  labs(x='Time', y='Response')

Наконец, обратите внимание, что использование separate() в df1 выдаст предупреждение, хотя и не является ошибкой, поскольку есть дополнительное значение, созданное путем разделения на sep="_".Вы можете избежать этого (если это вызывает путаницу), добавив уровень для захвата окончательного значения (которое выглядит как индекс времени):

separate(var, c("Genotype", "Tissue", "Time", "Index"), sep = "_")