Data Loading and Tidying

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

metadata         <- list()
metadata$patient <- 
  read_tsv_quiet(paths$patient_md) %>% 
  mutate(clinical_variant = fct_relevel(clinical_variant, "Sporadic"))

metadata$genome  <- read_tsv_quiet(paths$genome_md, col_types = date_types)
metadata$mrna    <- read_tsv_quiet(paths$mrna_md, col_types = date_types)

# Append site information to genome and mrna metadata
metadata <- map_at(metadata, c("genome", "mrna"), append_info, "site", "clinical_variant")

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

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

metadata %<>% map_if(~ "ff_or_ffpe" %in% names(.x), remove_ffpe)

testthat::expect_true(all(table(metadata$mrna$patient) == 1))
testthat::expect_true(all(table(metadata$genome$patient) == 2))

biomart <- useEnsembl("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl",
                               host = "oct2016.archive.ensembl.org")

symbols <-
  getBM(c("ensembl_gene_id", "hgnc_symbol"), mart = biomart) %>% 
  rename(gene_id = ensembl_gene_id, symbol = hgnc_symbol) %>% 
  distinct() %>% 
  mutate(symbol = ifelse(symbol == "", gene_id, symbol)) %>% 
  group_by(gene_id) %>% 
  slice(1)

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)) %>% 
  left_join(symbols, by = "gene_id") %>% 
  mutate(symbol = ifelse(grepl("ebv_|_PAR_Y", transcript), gene, symbol)) %>% 
  select(transcript, symbol, gene_id) %>% 
  as.data.frame()

# Useful functions and variables
ebv_genes   <- filter(tx2gene, grepl("^ebv_", transcript)) %$% symbol
get_gene_id <- function(name) filter(tx2gene, symbol %in% name) %$% unique(gene_id)
get_gene    <- function(id) slice(tx2gene, match(id, gene_id)) %$% symbol

genes         <- list()
genes$bl      <- read_lines(paths$bl_genes)
genes$mbl     <- read_lines(paths$mbl_genes) %>% intersect(tx2gene$symbol)
genes$morgan  <- read_lines(paths$morgan_genes) %>% intersect(tx2gene$symbol)
genes$wright  <- read_lines(paths$wright_genes) %>% get_gene()
genes$malaria <- read_lines(paths$malaria_genes) %>% get_gene()

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),
    ebv_type   = ifelse(clinical_variant == "Centroblast", "None", ebv_type)) %>% 
  select(patient, ebv_status, ebv_type)

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

rm(status)

metadata$genome <- 
  read_tsv_quiet(paths$tc) %>% 
  select(biospecimen_id = Tumor_Sample_Barcode, tumour_content = Tumour_Content) %>% 
  left_join(metadata$genome, ., by = "biospecimen_id")

maf <- read_maf(paths$maf)

# Update variant classification with simpler categories
maf@data$Variant_Classification <- get_categories(maf@data$Consequence)
maf@maf.silent$Variant_Classification <- get_categories(maf@maf.silent$Consequence)

# Restrict to cases with metadata
maf <- 
  metadata$genome %$% 
  gsub("-", ".", biospecimen_id) %>% 
  subsetMaf(maf, tsb = ., includeSyn = TRUE, mafObj = TRUE) %>% 
  process_maf()

# Estimate tumour content
maf@data <- 
  estimate_tc(maf@data) %>% 
  as.data.table() %>% 
  setkey(patient) %>% 
  merge(setkey(maf@data, patient), .)

# Add FF/FFPE status and tumour_content to MAF
maf@data <- 
  metadata$genome %>% 
  filter(tissue_status == "Tumor") %>% 
  select(patient, ff_or_ffpe, ebv_type, clinical_variant, sex) %>% 
  as.data.table() %>% 
  setkey(patient) %>% 
  merge(maf@data, ., all = FALSE)

maf_all <- suppressWarnings(fread(paths$maf_all, showProgress = FALSE))
  
maf_all[, patient := get_patient_id(Tumor_Sample_Barcode)] %>% invisible()
maf_all[, t_vaf := t_alt_count / (t_alt_count + t_ref_count)] %>% invisible()
setkey(maf_all, patient)

maf_all <- 
  estimate_tc(maf_all) %>% 
  as.data.table() %>% 
  setkey(patient) %>% 
  merge(maf_all, .)

maf_all <- 
  metadata$genome %>% 
  filter(biospecimen_id %in% maf_all$Tumor_Sample_Barcode) %>% 
  select(patient, ff_or_ffpe) %>%
  distinct() %>% 
  as.data.table() %>% 
  setkey(patient) %>% 
  merge(maf_all, .)

mmaf <- read_maf(paths$mmaf)

# Update variant classification with simpler categories
mmaf@data$Variant_Classification       <- get_categories(mmaf@data$Consequence)
mmaf@maf.silent$Variant_Classification <- get_categories(mmaf@maf.silent$Consequence)

mmaf <- subsetMaf(
  mmaf, query = "Variant_Classification %in% c('Truncation', 'Missense', 'Splicing')", 
  includeSyn = TRUE, mafObj = TRUE)

setkey(mmaf@data, Tumor_Sample_Barcode)

mmaf@data <-
  metadata$genome %>% 
  select(Tumor_Sample_Barcode = biospecimen_id, clinical_variant, ebv_type) %>% 
  mutate(Tumor_Sample_Barcode = gsub("-", ".", Tumor_Sample_Barcode)) %>% 
  as.data.table() %>% 
  setkey(Tumor_Sample_Barcode) %>% 
  merge(mmaf@data, ., all.x = TRUE) %>% 
  replace_na(list(clinical_variant = "Sporadic", ebv_type = "None"))

smgs <- 
  read_tsv_quiet(paths$smgs, col_names = c("gene", "num_methods")) %>% 
  filter(num_methods > 1) %$%
  gene

ffpe_samples <- metadata$genome %$% biospecimen_id[ff_or_ffpe == "FFPE"]

sequenza <- 
  Sys.glob(file.path(paths$sequenza, "*", "*_segments.txt")) %>% 
  set_names(get_sample_id(.)) %>% 
  map_df(read_tsv_quiet, .id = "biospecimen_id") %>% 
  rename(start = start.pos, end = end.pos) %>% 
  filter(
    !biospecimen_id %in% ffpe_samples,
    chromosome %in% paste0("chr", seq(1, 22)),
    !is.na(CNt)) %>% 
  mutate(patient = get_patient_id(biospecimen_id))

manta_raw <- 
  Sys.glob(paths$manta) %>% 
  set_names(get_sample_id(.)) %>% 
  map(VariantAnnotation::readVcf) %>% 
  map(filter_pass) %>% 
  map(calc_tvaf) %>% 
  map_df(vcf_to_df, .id = "biospecimen_id") %>% 
  as_tibble() %>% 
  inner_join(select(metadata$genome, biospecimen_id, patient, ff_or_ffpe), 
             by = "biospecimen_id") %>% 
  inner_join(estimate_tc(maf@data), by = "patient")

annotations <- 
  metadata$mrna %>% 
  mutate(cell_sorting = ifelse(is_cell_sorted(biospecimen_id), "Sorted", "Unsorted")) %>% 
  select(patient, clinical_variant, ebv_type, sex, ff_or_ffpe, 
         cell_sorting, site, tissue, lib_date) %>% 
  mutate_at(c("lib_date"), function(x) substr(x, 1, 7)) %>% 
  mutate_if(function(x) !is.factor(x), funs(fill_na)) %>% 
  mutate_each("as.factor", -patient) %>% 
  mutate(
    clinical_variant = fct_relevel(clinical_variant, "Sporadic", "Endemic"),
    ebv_type = fct_relevel(ebv_type, "Type 1", "Type 2")) %>% 
  as.data.frame() %>% 
  column_to_rownames("patient")

gannotations <- 
  metadata$genome %>% 
  filter(tissue_status == "Tumor") %>% 
  mutate(cell_sorting = ifelse(is_cell_sorted(biospecimen_id), "Sorted", "Unsorted")) %>% 
  select(patient, clinical_variant, ebv_type, sex, 
         ff_or_ffpe, cell_sorting, site, tissue) %>% 
  mutate_if(function(x) !is.factor(x), funs(fill_na)) %>% 
  mutate_each("as.factor", -patient) %>% 
  mutate(
    clinical_variant = fct_relevel(clinical_variant, "Sporadic", "Endemic"),
    ebv_type = fct_relevel(ebv_type, "Type 1", "Type 2")) %>% 
  as.data.frame() %>% 
  column_to_rownames("patient")

colours <- get_legend_colours(annotations)

display_colours(colours, 5)

# 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$col_data   <- annotations[colnames(salmon$txi$counts),]
salmon$raw$dds    <- DESeqDataSetFromTximport(salmon$txi, salmon$col_data, ~ 1)
salmon$raw$dds    <- estimateSizeFactors(salmon$raw$dds)
salmon$raw$counts <- counts(salmon$raw$dds, normalized = TRUE)