Metadata

Sample metadata

date_types <- cols(lib_date = "c", seq_date = "c")

metadata         <- list()
metadata$patient <- read_tsv_quiet(paths$patient_md)
metadata$mrna    <- read_tsv_quiet(paths$mrna_md, col_types = date_types)

# Append site information to mRNA metadata
metadata$mrna <- 
  metadata$patient %>% 
  select(patient, site, clinical_variant) %>% 
  left_join(metadata$mrna, ., by = "patient")

We have incomplete HIV status for this cohort. For this reason, we only consider the sporadic and endemic variants. We will exclude any HIV-positive cases (that we know of).

hiv_patients <- filter(metadata$patient, hiv_status == "Positive") %$% patient

metadata <- purrr::map(metadata, ~filter(.x, !patient %in% hiv_patients))

There happens to be a few patients with more than one RNA-seq dataset because both the FF and FFPE tumor tissue underwent library construction and sequencing. For now, we’ll simply be excluding the FFPE samples when an FF sample is available for the same patient.

metadata$mrna %<>% 
  group_by(patient) %>% 
  mutate(num_ff = sum(ff_or_ffpe == "FF")) %>% 
  filter(xor(num_ff > 0, ff_or_ffpe == "FFPE")) %>% 
  ungroup() %>% 
  select(-num_ff)

testthat::expect_false(any(duplicated(metadata$mrna$patient)))

Reference

Gene/transcript metadata

This gene-to-transcript mapping will allow tximport to summarize the transcript counts reported by Salmon at the gene level. We also use this data frame to convert gene IDs into gene symbols.

tx2gene_cols <- c("transcript", "gene_id", "gene")
tx2gene      <- 
  read_tsv_quiet(paths$tx2gene, col_names = tx2gene_cols) %>%
  mutate(gene_id = fix_gene_ids(gene_id))

# Useful functions
id2gene     <- select(tx2gene, gene_id, gene) %>% distinct()
get_gene_id <- function(name) 
                 filter(id2gene, gene %in% name) %$% unique(gene_id)
get_gene    <- function(id) 
                 slice(id2gene, match(id, gene_id)) %$% gene

Gene lists

genes        <- list()
genes$mbl    <- read_lines(paths$mbl_genes)
genes$morgan <- read_lines(paths$morgan_genes)

geneids        <- map(genes, ~get_gene_id(.x))
geneids$wright <- read_lines(paths$wright_genes)

genes <- map(geneids, ~get_gene(.x))

Results

Sex and EBV status

We inferred sex and EBV status from the sequencing data because the clinical annotations were incomplete.

status <- list()

status$sex <- 
  read_tsv_quiet(paths$sex_status) %>% 
  mutate(patient = get_patient_id(sample)) %>% 
  select(-sample, -chrY_chrX_read_ratio)

status$ebv <- 
  read_tsv_quiet(paths$ebv_status) %>%
  mutate(patient = get_patient_id(sample)) %>% 
  select(-sample, -ebv_read_ratio)

# Check if inferred sex status is consistent with clinical annotations
is_sex_consistent <- 
  metadata$patient %>% 
  inner_join(status$sex, by = "patient") %$% 
  ifelse(!is.na(annotated_sex), annotated_sex == sex, TRUE)
testthat::expect_true(all(is_sex_consistent))

# Fill in any blanks using the clinical annotations
status$sex <- 
  metadata$patient %>% 
  left_join(status$sex, by = "patient") %>% 
  mutate(sex = ifelse(is.na(sex), annotated_sex, sex)) %>% 
  select(patient, sex)

status$ebv <- 
  metadata$patient %>% 
  left_join(status$ebv, by = "patient") %>% 
  mutate(ebv_status = ifelse(clinical_variant == "Centroblast", 
                             "Negative", ebv_status)) %>% 
  select(patient, ebv_status)

metadata <- map(metadata, ~left_join(.x, status$sex, by = "patient"))
metadata <- map(metadata, ~left_join(.x, status$ebv, by = "patient"))

rm(status)

Annotations

Now that the sex and EBV status are incorporated into the metadata data frames, we can create an annotations data frame. This is useful for certain tools like DEseq2 and pheatmap.

annotations <- 
  metadata$mrna %>% 
  select(
    patient, clinical_variant, tissue_status, ebv_status, 
    sex, ff_or_ffpe, site, lib_date) %>% 
  mutate_at(c("lib_date"), function(x) substr(x, 1, 7)) %>% 
  mutate_each(funs(fill_na), -patient) %>% 
  mutate_each("as.factor", -patient) %>% 
  as.data.frame() %>% 
  column_to_rownames("patient")

Legend colours

Another thing we can do now that the metadata is loaded is define the legend colours.

colours <- get_legend_colours(annotations)

display_colours(colours, 5)

Salmon

We measured expression using Salmon (quasi-alignment). The transcriptome we used it Gencode v25 supplemented with EBV transcripts.

# Load Salmon data
salmon     <- list()
salmon$txi <- load_salmon(paths$salmon, tx2gene, metadata$mrna$biospecimen_id)

# Load raw data as DESeq2 dataset
salmon$raw        <- list()
salmon$raw$dds    <- DESeqDataSetFromTximport(salmon$txi, 
                       annotations[colnames(salmon$txi$counts),], ~ 1)
salmon$raw$dds    <- estimateSizeFactors(salmon$raw$dds)
salmon$raw$counts <- counts(salmon$raw$dds, normalized = TRUE)