<p>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 <a href="https://bschiffthaler.s3-eu-west-1.amazonaws.com/Public/sebio_14571.html" rel="nofollow noreferrer">here</a>)</p>
<pre><code>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<span class="math-container">$Sample), " +")
classification &lt;- supp_data$</span>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<span class="math-container">$EmeanDiff, xdata$</span>Classification)
wilcox.test(x = xdata<span class="math-container">$`Epithelial-like`, y = xdata$</span>`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
<span class="math-container">```</span> 
</code></pre>

CC BY-SA 4.0
2020-10-18T14:19:31.227Z
Hi @CSP it is certainly an interesting discussion and some good R coding. Could you look through this very detailed reply and comment please?
M__
CC BY-SA 4.0
2020-10-19T08:41:37.000Z
You were right, I had a minor mistake in my pipeline! Thank you for this answer, very detailed and useful!
CSP
CC BY-SA 4.0
2020-10-19T12:24:59.987Z
Great that it worked. If this is answered, please accept the answer
Bastian Schiffthaler
CC BY-SA 4.0
2020-10-19T12:33:27.353Z
Agreed @CSP please accept the above answer
M__
CC BY-SA 4.0
2020-10-18T07:43:36.243Z
2020-10-18T07:43:36.243Z
Bastian Schiffthaler
<p>The most likely explanation here is paired versus unpaired.</p>
<p>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.</p>

CC BY-SA 4.0
2020-10-18T07:46:10.460Z
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.
Bastian Schiffthaler
CC BY-SA 4.0
2020-10-18T07:47:01.020Z
Also they cannot do a paired test, since they have uneven groups (Epithelial vs Mesenchymal)
Bastian Schiffthaler
CC BY-SA 4.0
2020-10-18T14:13:05.330Z
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.
M__
CC BY-SA 4.0
2020-10-19T08:43:21.617Z
The samples are unpaired, it was a mistake in my prior pipeline, but thank you for you very useful input!
CSP
CC BY-SA 4.0
2020-10-17T19:46:38.447Z
2020-10-17T19:46:38.447Z
M__
<p>I am trying to replicate the results of <a href="https://www.nature.com/articles/onc2014136" rel="nofollow noreferrer">this</a> paper roughly guided by <a href="https://www.bioconductor.org/packages/devel/workflows/vignettes/maEndToEnd/inst/doc/MA-Workflow.html#" rel="nofollow noreferrer">this</a> 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:</p>
<blockquote>
<p>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).</p>
</blockquote>
<p>My p-values for a given gene are astronomically different to that of the paper...
Here is my code:</p>
<p>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).</p>
<pre><code> 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)
</code></pre>
<p>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</p>

CC BY-SA 4.0
2020-10-16T15:24:05.307Z
Cross-posted on biostars: https://www.biostars.org/p/467695/
Ram RS
CC BY-SA 4.0
2020-10-16T19:18:10.963Z
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.
M__
<r><microarray>
Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper
# Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper

(Asked by: username_0 on 2020-10-16)

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

**Tags:** r, microarray
> **Comments:**
> **username_1**: Cross-posted on biostars: https://www.biostars.org/p/467695/  
> **username_2**: 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.  

---

## 2 answers

### Answer by username_3 on 2020-10-18. Score: 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()
```

> **Comments:**
> **username_2**: Hi @CSP it is certainly an interesting discussion and some good R coding. Could you look through this very detailed reply and comment please?  
> **username_0**: You were right, I had a minor mistake in my pipeline! Thank you for this answer, very detailed and useful!  
> **username_3**: Great that it worked. If this is answered, please accept the answer  
> **username_2**: Agreed @CSP please accept the above answer  

---

### Answer by username_2 on 2020-10-17. Score: 0

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.

> **Comments:**
> **username_3**: 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.  
> **username_3**: Also they cannot do a paired test, since they have uneven groups (Epithelial vs Mesenchymal)  
> **username_2**: 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.  
> **username_0**: The samples are unpaired, it was a mistake in my prior pipeline, but thank you for you very useful input!  

---

