Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper

Question

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 >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<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[,"Classification"]=="Epithelial-like"))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Mesenchymal-like"))

epi_values = expr_f["8143663", epi_samples]
mes_values = expr_f["8143663", 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

Answer 1

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&M (I didn't filter the "fragile probes" 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&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 <- getGEO("GSE47856", GSEMatrix =TRUE, getGPL=FALSE)
adf <- 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&format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles <- sapply(sampleNames(adf), function(f) {
  pat <- paste0(f, ".*.CEL.gz$")
  dir("~/Downloads/cel", pattern = pat, full.names = TRUE)
})
raw_data <- oligo::read.celfiles(celfiles,
                                 phenoData = adf, 
                                 experimentData = experimentData(gset[[1]]))

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

#' Use ComBat for batch effect removal
batch_data <- phenoData(rma_expr)[["batch:ch1"]]
cb_expr <- sva::ComBat(rma_expr_assay, batch_data)

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

#' 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 <- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[["cell line:ch1"]]),
    " +"),
  "^tay$", "taya"),
  "dov13[ab]", "dov13"
)

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

#' Get the treament from the phenoData object
treatment <- adf[["treatment:ch1"]]

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

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

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

#' Session info
sessionInfo()
``` 

Answer 2

The most likely explanation here is paired versus unpaired.

Mann-Whitney is normally associated with unpaired data, but the authors may have performed a paired test. You performed an unpaired test and thus there is a massive discrepancy in results. The authors should have stated the effect of pairing - it that is what they did - but they are under no strict obligations to do so. You however have performed a Wilcoxan rank sign test, albeit the results shouldn't be so different.