ggloor / ALDEx_bioc

ALDEx_bioc is the working directory for updating bioconductor
27 stars 13 forks source link

Rearranging dataframe changes output and Benjamini-Hochberg corrected p-value estimates differ between Bioconductor and developer version #70

Open wiederkf opened 8 months ago

wiederkf commented 8 months ago

Hi there,

I ran into two issues and would appreciate your feedback.

Issue 1

In my testing script below, rearranging the dataframe changes the output (Case 1 vs Case 2 and Case 3 vs Case 4). Why is this?

Issue 2

While the estimated p-values between the Bioconductor (Cases 1 and 2) and developer versions (Cases 3 and 4) are similar, I obtain vastly different Benjamini-Hochberg corrected p-values. Why is this? Which should I trust?

Versions

Testing script

require(tidyverse)

# Bioconductor version

remove.packages("ALDEx2")
BiocManager::install("ALDEx2")

require(ALDEx2)
data(selex)

## Case 1
set.seed(2)
selex_sub_1<-selex[1200:1600,]
conds_1<-c(rep("NS",7),rep("S",7))
x_1<-aldex.clr(selex_sub_1,conds_1,mc.samples=16)
x_tt_1<-aldex.ttest(x_1)
x_effect_1<-aldex.effect(x_1)
x_all_1<-data.frame(x_tt_1,x_effect_1)
x_all_1[1:5,]

## Case 2
set.seed(2)
selex_sub_2<-selex_sub_1 %>%
  dplyr::select(8:14,1:7)
conds_2<-c(rep("S",7),rep("NS",7))
x_2<-aldex.clr(selex_sub_2,conds_2,mc.samples=16)
x_tt_2<-aldex.ttest(x_2)
x_effect_2<-aldex.effect(x_2)
x_all_2<-data.frame(x_tt_2,x_effect_2)
x_all_2[1:5,]

# Developer version

detach("package:ALDEx2",unload=TRUE)
remove.packages("ALDEx2")
devtools::install_github("ggloor/ALDEx_bioc")

require(ALDEx2)
data(selex)

## Case 3
set.seed(2)
selex_sub_3<-selex[1200:1600,]
conds_3<-c(rep("NS",7),rep("S",7))
x_3<-aldex.clr(selex_sub_3,conds_3,mc.samples=16)
x_tt_3<-aldex.ttest(x_3)
x_effect_3<-aldex.effect(x_3)
x_all_3<-data.frame(x_tt_3,x_effect_3)
x_all_3[1:5,]

## Case 4
set.seed(2)
selex_sub_4<-selex_sub_3 %>%
  dplyr::select(8:14,1:7)
conds_4<-c(rep("S",7),rep("NS",7))
x_4<-aldex.clr(selex_sub_4,conds_4,mc.samples=16)
x_tt_4<-aldex.ttest(x_4)
x_effect_4<-aldex.effect(x_4)
x_all_4<-data.frame(x_tt_4,x_effect_4)
x_all_4[1:5,]

Output

Cases 1-4