Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper I am trying to replicate the results of this paper roughly guided by this pipeline.
Basically, we are trying to detect a differential expression between mesenchymal and epithelial cells under cisplatin treatment vs control.
However, I am facing a problem when trying to calculate the p-value for the Mann-Whitney U-test. This is the part of the paper I am referring to:

Robust Multichip Average (RMA) normalization was performed at the
transcript level on the results from the Affymetrix Human Gene ST1.0
arrays using Affymetrix Power Tool 1.14.4.1 for all 46 sham- or
cisplatin-treated ovarian cancer cell lines. The normalized data were
subsequently standardized using ComBat71 to remove the batch effect.
In this experiment, the cisplatin treatment assay was performed in
triplicate on 20 cell lines, while single assays (without replicate)
were performed on the remaining 26 cell lines. Taking advantage of the
triplicate data, potentially fragile probes with strong variations (an
s.d. of &gt;0.2) within the triplicates were removed, decreasing the
probe number from 33 297 to 21 329. To perform a fair comparison, the
triplicate data were then log-averaged into one value so that one
result for each cell line could be used in the following analyses. To
detect differential responses to cisplatin between epithelial- and
mesenchymal-like cell lines, the transcriptomic responses to cisplatin
were computed by subtracting the gene expression value of control
(cisplatin untreated) cells from that of cisplatin-treated cells.
Mann–Whitney U-test (P&lt;0.01 as a cut-off value) was subsequently used
to detect the differential transcriptomic responses between the
expression changes by cisplatin treatment in epithelial-like cell
lines with those in mesenchymal-like cell lines (Supplementary Table
3).

My p-values for a given gene are astronomically different to that of the paper...
Here is my code:
expr_f is the expression values matrix from the Expression set (log-transformed).
SDRF_f is the pData from the same expression set.
I am feeding into the wilcox.test function a vector of expression from epithelial cells (epi_values) and mesencymal cells (mes_values) for the same probe (8143663).
 genelist=rownames(expr_f)
epi_samples = rownames(subset(SDRF_f, SDRF_f[,&quot;Classification&quot;]==&quot;Epithelial-like&quot;))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,&quot;Classification&quot;]==&quot;Mesenchymal-like&quot;))

epi_values = expr_f[&quot;8143663&quot;, epi_samples]
mes_values = expr_f[&quot;8143663&quot;, mes_samples]


wilcox.test(mes_values, epi_values, exact=TRUE)

In the paper, their p-value for this probe is 0.000372259258165627 while mine is 0.8785. The difference makes no sense at all to me.
I don't know if I am not quite getting the concept of the u-test, if I am using the function wrong, or if it's something else. Any input is appreciated
I think you may need to check your pipeline in this regard. I ran the analysis and replicated as closely as I could their M&amp;M (I didn't filter the &quot;fragile probes&quot; since you anyway care only about a specific probe). In any case, I get something relatively close to their result (P=0.004059), which is still one order of magnitude off, but without more detailed M&amp;M, this is probably as close as I can get to their analysis. Here's what I did (if you prefer to read the notebook html, that is here)
library(Biobase)
library(GEOquery)
library(sva) # ComBat
library(oligo)
library(googlesheets4)
library(stringr)
library(dplyr)

#' Get the GEO dataset mostly for the phenoData annotation
gset &lt;- getGEO(&quot;GSE47856&quot;, GSEMatrix =TRUE, getGPL=FALSE)
adf &lt;- Biobase::phenoData(gset[[1]])
sampleNames(adf)

#' Read RAW cel files from the GEO data supplementary rar archive
#' https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE47856&amp;format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles &lt;- sapply(sampleNames(adf), function(f) {
  pat &lt;- paste0(f, &quot;.*.CEL.gz$&quot;)
  dir(&quot;~/Downloads/cel&quot;, pattern = pat, full.names = TRUE)
})
raw_data &lt;- oligo::read.celfiles(celfiles,
                                 phenoData = adf, 
                                 experimentData = experimentData(gset[[1]]))

#' Get the raw expression and RMA on transcript clusters
rma_expr &lt;- oligo::rma(raw_data, target = &quot;core&quot;)
rma_expr_assay &lt;- assayData(rma_expr)$exprs

#' Use ComBat for batch effect removal
batch_data &lt;- phenoData(rma_expr)[[&quot;batch:ch1&quot;]]
cb_expr &lt;- sva::ComBat(rma_expr_assay, batch_data)

#' Get supplementary table 1 from the paper (which I copied on this Google Sheet)
supp_data &lt;- read_sheet(&quot;https://docs.google.com/spreadsheets/d/142__1O35zvlttejvz62XQX33uLRJhqIDvb604smVSRo&quot;)

#' Match the 'cell line' phenotype from the phenoData to the sample info of
#' the supplementary table and extract the cell classification. Also remove
#' all white space and convert everything to lower case for easier matching.
#' Some cell lines also have a/b attached, so I'm dropping that to match them
#' up. Mayb not be correct to do that here and treat all a/b the same (?)
cell_line &lt;- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[[&quot;cell line:ch1&quot;]]),
    &quot; +&quot;),
  &quot;^tay$&quot;, &quot;taya&quot;),
  &quot;dov13[ab]&quot;, &quot;dov13&quot;
)

sample_info &lt;- str_remove(str_to_lower(supp_data$Sample), " +")
classification &lt;- supp_data$Classification[match(cell_line, sample_info)]

#' Get the treament from the phenoData object
treatment &lt;- adf[[&quot;treatment:ch1&quot;]]

#' We'll need the geometric mean later (R has no built-in)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x &gt; 0]), na.rm=na.rm) / length(x))
}

#' Set up a tibble (==data frame) with all that info for your probe
tibble(Expression = cb_expr[&quot;8143663&quot;, ], 
       Line = cell_line, 
       Treatment = treatment, 
       Classification = classification) %&gt;%
  # Summarise replicated data with the geometric mean
  group_by(Line, Treatment, Classification) %&gt;%
  summarise(Emean = gm_mean(Expression)) %&gt;% 
  # Make sure data is in order
  arrange(Line, Treatment, Classification)%&gt;%
  # Calculate the difference between the cell types
  group_by(Line, Classification) %&gt;%
  summarise(EmeanDiff = diff(Emean)) -&gt; xdata

#' Test the difference
xdata &lt;- split(xdata$EmeanDiff, xdata$Classification)
wilcox.test(x = xdata$`Epithelial-like`, y = xdata$`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
```