Identifying differentially expressed latent variables in medulloblastoma

J. Taroni 2018

Here, we’re testing for latent variable differential expression between medulloblastoma subgroups (Group 3, Group 4, and SHH).

As in 37-medulloblastoma_recount2_model, where we prepped the data we are analyzing here, this is using Northcott, et al. and Robinson, et al. data.

Set up

# pipe is required for LVTestWrapper
`%>%` <- dplyr::`%>%`

We have several custom functions that we’ve written and previously used in our ANCA-associated vasculitis analyses. They are general enough that we can use them again here.

source(file.path("util", "test_LV_differences.R"))

Directories

# directories specific to this notebook
plot.dir <- file.path("plots", "38")
dir.create(plot.dir, recursive = TRUE, showWarnings = FALSE)
results.dir <- file.path("results", "38")
dir.create(results.dir, recursive = TRUE, showWarnings = FALSE)

Read in files

We need two things to do this analysis: a B matrix (latent variable values) and a data.frame that contains the sample labels for the groups that we’d like to test for differences between.

GSE37382 (Northcott, et al.)

northcott.b.file <- file.path("results", "37", "GSE37382_recount2_B.RDS")
northcott.sample.file <- file.path("data", "sample_info", 
                                   "GSE37382_cleaned_metadata.tsv")

Read in B (MultiPLIER)

northcott.b.matrix <- readRDS(northcott.b.file)

Read in the data.frame with the subgroup information and check that the sample names match between these two files.

northcott.sample.df <- readr::read_tsv(northcott.sample.file)

Parsed with column specification:
cols(
  source_name = col_character(),
  sample_identifier = col_character(),
  tissue_type = col_character(),
  age = col_character(),
  histology = col_character(),
  sex = col_character(),
  subgroup = col_character()
)

all(northcott.sample.df$source_name %in% colnames(northcott.b.matrix))

[1] TRUE

LVTestWrapper requires that the sample names/identifiers are in a column named Sample, so we’ll change source_name to Sample.

colnames(northcott.sample.df)[1] <- "Sample"

GSE37418 (Robinson, et al.)

We’ll repeat that process for the Robinson, et al. data

robinson.b.file <- file.path("results", "37", "GSE37418_recount2_B.RDS")
robinson.sample.file <- file.path("data", "sample_info", 
                                  "metadata_GSE37418.tsv")

robinson.b.matrix <- readRDS(robinson.b.file)
robinson.sample.df <- readr::read_tsv(robinson.sample.file)

Parsed with column specification:
cols(
  .default = col_character(),
  channel_count = col_integer(),
  contact_phone = col_double(),
  `contact_zip/postal_code` = col_integer(),
  data_row_count = col_integer(),
  taxid_ch1 = col_integer()
)
See spec(...) for full column specifications.

The subgroup information in this data set is a bit different, let’s look at the counts

robinson.sample.df %>% 
  dplyr::group_by(subgroup) %>%
  dplyr::tally()

To do differential expression analysis, we’ll remove the SHH OUTLIER and U samples

sample.index <- which(robinson.sample.df$subgroup %in% c("U",
                                                         "SHH OUTLIER"))
samples.to.remove <- robinson.sample.df$refinebio_accession_code[sample.index]
samples.to.remove

[1] "GSM918628" "GSM918618" "GSM918644"

remove.column.index <- which(colnames(robinson.b.matrix) %in% samples.to.remove)

Do the filtering

robinson.b.matrix <- robinson.b.matrix[, -remove.column.index]
robinson.sample.df <- robinson.sample.df %>%
  dplyr::filter(!(refinebio_accession_code %in% samples.to.remove)) %>%
  dplyr::select(c("refinebio_accession_code", "subgroup", "Sex", "age",
                  "title"))

all(robinson.sample.df$refinebio_accession_code %in% 
      colnames(robinson.b.matrix))

[1] TRUE

colnames(robinson.sample.df)[1] <- "Sample"

Test for differential expression

LVTestWrapper gives us 3 things: 1) differential expression results from limma, 2) boxplot + jitter plots of the latent variable expression and 3) the B matrix in “long” format – this is what is used for plotting.

We’ll use "BH" correction for multiple hypotheses testing (the default).

northcott.results <- 
  LVTestWrapper(b.matrix = northcott.b.matrix,
                sample.info.df = northcott.sample.df,
                phenotype.col = "subgroup",
                # the boxplot output, the "long" format B data.frame (useful for 
                # plotting), and the limma differential expression results
                # will be output in files that begin with this string
                file.lead = "GSE37382_subgroup_recount2_model",
                plot.dir = plot.dir,
                results.dir = results.dir)

Column `Sample` joining factor and character vector, coercing into character vector

robinson.results <- 
  LVTestWrapper(b.matrix = robinson.b.matrix,
                sample.info.df = robinson.sample.df,
                phenotype.col = "subgroup",
                file.lead = "GSE37418_subgroup_recount2_model",
                plot.dir = plot.dir,
                results.dir = results.dir)

Column `Sample` joining factor and character vector, coercing into character vector

Any overlap at all?

northcott.delvs <- 
  northcott.results$limma$LV[which(northcott.results$limma$adj.P.Val < 0.05)]

robinson.delvs <- 
  robinson.results$limma$LV[which(robinson.results$limma$adj.P.Val < 0.05)]

VennDiagram::venn.diagram(list(Robinson = robinson.delvs,
                               Northcott = northcott.delvs),
                          file.path(plot.dir, "Medulloblastoma_DELV_Venn.tiff"))

[1] 1