TCGA_survival_analysis

library(TCGAbiolinks)
library(biomaRt)

prepare_datasets <- function(project_name, filter = list(sample_type = c("Primary Tumor"))) {
  library(data.table)

  # from https://github.com/BioinformaticsFMRP/TCGAbiolinks/issues/493#issuecomment-1083311661
  # Defines the query to the GDC
  query <- GDCquery(project = project_name,
                    data.category = "Transcriptome Profiling",
                    data.type = "Gene Expression Quantification",
                    experimental.strategy = "RNA-Seq",
                    workflow.type = "STAR - Counts")
  
  # Get metadata matrix
  metadata <- query[[1]][[1]]
  
  # Download data using api
  GDCdownload(query, method = "api")
  
  # Get main directory where data is stored
  main_dir <- file.path("GDCdata", project_name)
  # Get file list of downloaded files
  file_list <- file.path("GDCdata", project_name,list.files(main_dir,recursive = TRUE)) 
  
  # Read first downloaded to get gene names
  test_tab <- read.table(file = file_list[1], sep = '\t', header = TRUE)
  # Delete header lines that don't contain usefull information
  test_tab <- test_tab[-c(1:4),]
  # STAR counts and tpm datasets
  tpm_data_frame <- data.frame(test_tab[,2])
  count_data_frame <- data.frame(test_tab[,2])
  
  # Append cycle to get the complete matrix
  for (i in c(1:length(file_list))) {
    # Read table
    test_tab <- fread(file = file_list[i],sep = "\t",header = T,data.table = F)
    # Delete not useful lines
    test_tab <- test_tab[-c(1:4),]
    # Column bind of tpm and counts data
    tpm_data_frame <- cbind(tpm_data_frame, test_tab[,9])
    sample_id = strsplit(file_list[i],split = "/")[[1]][6]
    sample_name = metadata[metadata$id == sample_id ,"cases.submitter_id"]
    colnames(tpm_data_frame)[i+1] = sample_name # col 1 is gene
    # count_data_frame <- cbind(count_data_frame, test_tab[,4])
    # Print progres from 0 to 1
    print(i/length(file_list))
  }
  # update col names
  colnames(tpm_data_frame)[1] = "Gene"
  genes = tpm_data_frame$Gene
  tpm_data_frame$Gene <- NULL
  
  
  
  # remove all 0 genes
  not_all_zero <-apply(tpm_data_frame, 1, function(row) any(row != 0))
  tpm_data_frame = tpm_data_frame[not_all_zero,]
  genes = genes[not_all_zero]
  
  
  # remove all rows that are not gene symbols
  ensembl = useMart(biomart="ensembl", dataset="hsapiens_gene_ensembl")
  results <- getBM(attributes=c("ensembl_gene_id","hgnc_symbol","transcript_biotype"), mart=ensembl)
  not_symbol = !genes %in% results$hgnc_symbol
  genes = genes[!not_symbol]
  tpm_data_frame = tpm_data_frame[!not_symbol,]
  rownames(tpm_data_frame) = make.unique(genes) 
  
  for (i in 1:length(filter)) {
    col = names(filter)[i]
    values = filter[[i]]
    primary_samples = metadata[metadata[[col]] %in% values,"cases.submitter_id"]
    tpm_data_frame = tpm_data_frame[,primary_samples]
  }
  return(tpm_data_frame)
}

set_clinical_data <- function(clin_data, genes, tpm_data_frame) {
  if (length(genes) == 1) {
    type = genes
  } else {
    type = "score"
  }
  
  gene_data = tpm_data_frame[rownames(tpm_data_frame) %in% genes, ]
  genes_average = as.data.frame(rowMeans(as.data.frame(t(gene_data))))
  colnames(genes_average) = "mean_score"
  gene_status = genes_average %>% mutate (
    gene_status = case_when(
      mean_score > quantile(mean_score)[4] ~ paste("High", type),
      mean_score < quantile(mean_score)[2] ~  paste("Low", type),
      TRUE ~ "same"
    )
  )
  rownames(gene_status) = gsub(x = rownames(gene_status),
                               pattern = "\\.",
                               replacement = "-")
  gene_status = gene_status[gene_status$gene_status != "same",]
  return(gene_status)
}

1 CECS

project_name = "TCGA-CESC"

# debugonce(prepare_datasets)
cesc_tpm = prepare_datasets(project_name = "TCGA-CESC")

clinical_data <- GDCquery_clinic(project_name, "clinical")
CESC_tissues = c("Basaloid squamous cell carcinoma",
"Papillary squamous cell carcinoma",
"Squamous cell carcinoma, keratinizing, NOS",
"Squamous cell carcinoma, large cell, nonkeratinizing, NOS",
"Squamous cell carcinoma, NOS")
clinical_data = clinical_data[clinical_data$primary_diagnosis %in% CESC_tissues,]

1.1 MYB

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = cesc_tpm)
gene_status

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10,
    width=10,filename = NULL
)
p[[1]]/p[[2]]

1.2 HPV+ HMSC signature

genes = "ANLN
TUBB6
AXL
IFT122
GFM1
SEPT2
PBRM1
INSIG1
ANXA2
CAPG"
genes = strsplit(x =genes,split = "\n")[[1]]

gene_status = set_clinical_data(clin_data = clinical_data,genes = genes,tpm_data_frame = cesc_tpm)
gene_status

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10,
    width=10,filename = NULL
)

p[[1]]/p[[2]]

2 OPSCC

project_name = "TCGA-HNSC"

HNSC_tpm = prepare_datasets(project_name = project_name,filter = list(sample_type = c("Primary Tumor")))

clinical_data <- GDCquery_clinic(project_name, "clinical")
OPSCC_tissues = c("Base of tongue, NOS",
"Oropharynx, NOS",
"Posterior wall of oropharynx",
"Tonsil, NOS")
clinical_data = clinical_data[clinical_data$tissue_or_organ_of_origin %in% OPSCC_tissues,]

2.1 MYB

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = HNSC_tpm)
gene_status

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10,
    width=10,filename = NULL
)

p[[1]]/p[[2]]

2.2 HPV+ HMSC signature

genes = "ANLN
TUBB6
AXL
IFT122
GFM1
SEPT2
PBRM1
INSIG1
ANXA2
CAPG"
genes = strsplit(x =genes,split = "\n")[[1]]

gene_status = set_clinical_data(clin_data = clinical_data,genes = genes,tpm_data_frame = HNSC_tpm)
gene_status

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10,
    width=10,filename = NULL
)
p[[1]]/p[[2]]