test_diff returns p-values in wrong order

[const-ae]

Hi Arne,

first thanks for the helpful package.

I was recently playing around with the test_diff function, comparing different imputation methods, when I noticed it sometimes returns the p-values in the wrong order. Below I have created a small reproducible example, to show the problem:

set.seed(1)
suppressPackageStartupMessages({
  library(SummarizedExperiment)
  library(DEP)
  library(tidyverse)
})

syn_data <- matrix(c(
  rep(42, times=6),
  c(1:3, 11:13),
  c(5:7, 8:6)), ncol=6, byrow = TRUE)


# This is important, because if the row names are sorted
# alphabetically there is no error!
rownames(syn_data) <- paste0("name", sample(1:nrow(syn_data)))
colnames(syn_data) <- c(paste0("cond_a-", 1:3), paste0("cond_b-", 1:3))
syn_data
#>       cond_a-1 cond_a-2 cond_a-3 cond_b-1 cond_b-2 cond_b-3
#> name2       42       42       42       42       42       42
#> name3        1        2        3       11       12       13
#> name1        5        6        7        8        7        6

syn_exp_df <- data.frame(sample=colnames(syn_data)) %>%
  mutate(label=sample) %>%
  separate(sample, into=c("condition", "replicate"), remove=FALSE, sep="-")

syn_row_df <- data.frame(name=rownames(syn_data), ID=seq_len(nrow(syn_data)))

syn_se <- SummarizedExperiment(syn_data, colData=syn_exp_df, rowData=syn_row_df)


syn_res <- test_diff(syn_se, type="manual", test="cond_a_vs_cond_b")
#> Tested contrasts: cond_a_vs_cond_b
#> Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There
#> may be too few input test statistics for reliable FDR calculations!
#> Warning: Censored sample for null model estimation has only size 2 !
rowData(syn_res)
#> DataFrame with 3 rows and 7 columns
#>           name        ID cond_a_vs_cond_b_CI.L cond_a_vs_cond_b_CI.R
#>       <factor> <integer>             <numeric>             <numeric>
#> name2    name1         3     -3.12424921339827      1.12424921339827
#> name3    name2         1  -0.00823518949045511   0.00823518949045511
#> name1    name3         2     -12.1242492133983     -7.87575078660174
#>       cond_a_vs_cond_b_diff cond_a_vs_cond_b_p.adj cond_a_vs_cond_b_p.val
#>                   <numeric>              <numeric>              <numeric>
#> name2                    -1      0.239304093291693      0.268028898798365
#> name3                     0      0.666666666666667                      1
#> name1                   -10   4.44089209850063e-16   0.000140843608180863
assay(syn_res)
#>       cond_a-1 cond_a-2 cond_a-3 cond_b-1 cond_b-2 cond_b-3
#> name2       42       42       42       42       42       42
#> name3        1        2        3       11       12       13
#> name1        5        6        7        8        7        6

^{Created on 2018-11-15 by the reprex package (v0.2.1)}

As you can see the rownames of rowData(syn_res) do not match the name column of that dataframe and the p-value of 1, which I would expect in the first row where each element if 42, appears in the second row.

From what I understand the problem is that in the test_diff() function you use the merge method, but don't set sort=FALSE.

Best, Constantin

[arnesmits]

Thanks a lot Constantin!