Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper
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()
```