cervical_TCGA

# print params
paged_table(as.data.frame(unlist(params)))

1 Data

library(TCGAbiolinks)
library(biomaRt)
library(kableExtra)
library(rmarkdown)
source("./Notebooks/TCGA/TCGA_functions.R")

# get hpv signature genes
scc_deg = readRDS(file = "./Data_out/scc_deg.rds")

avg_diff_cutoff = 0.6
fdr_cutoff = 0.1
signature = scc_deg %>% filter(avg_diff > avg_diff_cutoff &
                                  fdr < fdr_cutoff) %>% rownames()
signature

 [1] "AQP3"    "TXN"     "IGFBP2"  "RPS2"    "RPS5"    "MT-CO2"  "KRT19"   "KRT5"    "ATP1B3"  "MT-CYB"  "DEK"     "MT-CO3" 
[13] "S100A2"  "ITM2B"   "SFN"     "NUCKS1"  "IFI16"   "SERBP1"  "FABP5"   "RPS10"   "MT-ND4"  "HMGN2"   "TNFSF10" "MT-ATP6"
[25] "ID1"     "HMGB1"   "HES1"    "MT-ND3"  "PERP"    "TACSTD2" "JUN"     "CLDN4"   "MT-ND5"  "FOS"     "IFI27"   "MT1X"   
[37] "LY6D"    "RPS11"   "FOSB"    "DSP"     "MT-ND1"  "KRT15"  

project_name = "TCGA-CESC"

# cesc_tpm = prepare_datasets(project_name = "TCGA-CESC")
# saveRDS(cesc_tpm,file = "./Data_out/TCGA_survival_analysis/TCGA_CESC_TPM.RDS")
cesc_tpm = readRDS("./Data_out/TCGA_survival_analysis/TCGA_CESC_TPM.RDS")

if (params$log_tpm) {
  cesc_tpm = log2(cesc_tpm+1)
  
}

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,]

2 MYB

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = cesc_tpm,stratify = params$stratify)
paged_table(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,pval.method = T
)
p[[1]]/p[[2]]

Warning in grSoftVersion() :
  unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
  libXt.so.6: cannot open shared object file: No such file or directory

3 HPV+ signature

# load data 
library("readxl")
library(stringr)
library(ggpubr)
TCGA_HPV_data <- read_excel("./Input_data/TCGA/PMC3806554_ncomms3513-s2.xlsx",skip = 3)
TCGA_HPV_data  %<>% filter(`Tumor type` == "Primary solid")
TCGA_HPV_data$submitter_id = TCGA_HPV_data$`Sample barcode` %>% str_sub( start = 1, end = 12)

TCGA_HPV_data = TCGA_HPV_data %>%
    filter(Cancer == "CESC") %>%
    mutate (hpv_status = 
            if_else(
              condition = ppm == "N/A",
              true = "HPV-",
              false =  "HPV+"))
clinical_data = inner_join(x = clinical_data,y = TCGA_HPV_data,by = c("submitter_id"="submitter_id"))

p = TCGAanalyze_survival(
    data = clinical_data,
    clusterCol = "hpv_status",
    main = "TCGA Set\n GBM",
    filename = NULL,pval.method = T
)

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

4 MYB- HPV correlation

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

gene_status = gene_status[TCGA_HPV_data$submitter_id,] %>% cbind(HPV_PPM = as.numeric(TCGA_HPV_data$ppm))

Warning in eval(quote(list(...)), env) : NAs introduced by coercion

paged_table(gene_status)

sp <- ggscatter(gene_status, x = "MYB", y = "HPV_PPM",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   )
# Add correlation coefficient
sp + stat_cor(method = "pearson")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 12 rows containing non-finite values (stat_smooth).
Warning: Removed 12 rows containing non-finite values (stat_cor).
Warning: Removed 12 rows containing missing values (geom_point).

5 HPV signature-HPV correlation

gene_status = set_clinical_data(clin_data = clinical_data,genes = signature,tpm_data_frame = cesc_tpm,stratify = params$stratify,signature_name = "HPV+ signature")

gene_status = gene_status[TCGA_HPV_data$submitter_id,] %>% cbind(HPV_PPM = TCGA_HPV_data$ppm)
paged_table(gene_status)

sp <- ggscatter(gene_status, x = "HPV+ signature", y = "HPV_PPM",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   )
# Add correlation coefficient
sp + stat_cor(method = "pearson")

6 MYB- HPV T test

# create data

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = cesc_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","hpv_status")],by = c("patient"="submitter_id"))

paged_table(gene_status)

NA
NA
NA

# plot
library(ggpubr)
library(rstatix)
stat.test <- gene_status %>%
  t_test(MYB ~ hpv_status) %>%
  add_significance()
stat.test

bxp <- ggboxplot(gene_status, x = "hpv_status", y = "MYB", fill = "#00AFBB",add = "jitter")
stat.test <- stat.test %>% add_xy_position(x = "hpv_status")
bxp + 
  stat_pvalue_manual(stat.test,  label = "T-test, p = {p}") 

7 MYB and HPV survival plot

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = cesc_tpm,stratify = params$stratify)
gene_status$HPV_status   = if_else(
  condition = row.names(gene_status)  %in% TCGA_HPV_data$submitter_id,
  true = "HPV+",
  false = "HPV-"
)

gene_status$gene_status = paste(gene_status$gene_status,gene_status$HPV_status,sep = ", ")
gene_status$HPV_status = NULL
paged_table(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,pval.method = T
)
p[[1]]/p[[2]]

8 shared signature

shared_signature = readRDS(file = "./Data_out/temp/shared_signature.RDS")

gene_status = set_clinical_data(clin_data = clinical_data,genes = shared_signature,tpm_data_frame = cesc_tpm,stratify = params$stratify,signature_name = "HPV+ signature")
paged_table(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,pval.method = T
)

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

9 HMSC signature

# load
hmsc_deg = read.table(file = "./Data_out/HPV_analysis/hpv_deg_df.csv",row.names = 1)
  avg_diff_cutoff = 1
  fdr_cutoff = 0.1
  signature = hmsc_deg %>% filter(avg_diff>avg_diff_cutoff & fdr<fdr_cutoff) %>% rownames()
  signature

gene_status = set_clinical_data(clin_data = clinical_data,genes = signature,tpm_data_frame = cesc_tpm,stratify = params$stratify,signature_name = "HPV+ signature")
paged_table(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,pval.method = T
)

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