0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Cross-posted on biostars: biostars.org/p/467695 $\endgroup$ – Ram RS Oct 16 at 15:24
  • $\begingroup$ Wow that is some difference particularly for the Mann-Whitney. You are using Wilcox of course. Regardless there is one possibility which I will mention below if I get chance. $\endgroup$ – Michael Oct 16 at 19:18
1
$\begingroup$

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()
``` 
| improve this answer | |
$\endgroup$
  • $\begingroup$ Hi @CSP it is certainly an interesting discussion and some good R coding. Could you look through this very detailed reply and comment please? $\endgroup$ – Michael Oct 18 at 14:19
  • $\begingroup$ You were right, I had a minor mistake in my pipeline! Thank you for this answer, very detailed and useful! $\endgroup$ – CSP Oct 19 at 8:41
  • 1
    $\begingroup$ Great that it worked. If this is answered, please accept the answer $\endgroup$ – Bastian Schiffthaler Oct 19 at 12:24
  • $\begingroup$ Agreed @CSP please accept the above answer $\endgroup$ – Michael Oct 19 at 12:33
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ In R, a two-sample wilcox.test (i.e. providing both x and y) is the same as a Mann-Whitney U test. From the docs: >Performs one- and two-sample Wilcoxon tests on vectors of data; the latter is also known as ‘Mann-Whitney’ test. $\endgroup$ – Bastian Schiffthaler Oct 18 at 7:46
  • 1
    $\begingroup$ Also they cannot do a paired test, since they have uneven groups (Epithelial vs Mesenchymal) $\endgroup$ – Bastian Schiffthaler Oct 18 at 7:47
  • $\begingroup$ Okay fair point @BastianSchiffthaler , I'm wrong. I was aware of the ability of Wilcox to approximate Mann-Whitney ... if the samples are naturally unpaired ... again simply I'm wrong. Mann-Whitney is naturally unpaired so I was never really sure. The only thing I'd point out is that generally any bug can result in p<<0.01 but when p>>0.1 it difficult for a bug to achieve this. $\endgroup$ – Michael Oct 18 at 14:13
  • $\begingroup$ The samples are unpaired, it was a mistake in my prior pipeline, but thank you for you very useful input! $\endgroup$ – CSP Oct 19 at 8:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.