text: user7382: 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
  
  user9693: 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()
  ``` 
  
  
original_text: 7382: <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>
   9693: <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>
  