Survival Analysis

Colorectal Cancer | Pyroptosis | Unique Genes per RCD Type | Log Fold Change

I. Preliminaries

Loading libraries

library("tidyverse")
library("tibble")
library("msigdbr")
library("ggplot2")
library("TCGAbiolinks")
library("RNAseqQC")
library("DESeq2")
library("ensembldb")
library("purrr")
library("magrittr")
library("vsn")
library("matrixStats")
library("dplyr")
library("grex")
library("survminer")
library("survival")

II. Downloading the TCGA gene expression data

Create a function for downloading TCGA gene expression data.

For more detailed documentation, refer to 2. Differential Gene Expression Analysis - TCGA.Rmd.

query_and_filter_samples <- function(project) {
  query_tumor <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = "Primary Tumor"
  )
  tumor <- getResults(query_tumor)

  query_normal <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = "Solid Tissue Normal"
  )
  normal <- getResults(query_normal)

  submitter_ids <- inner_join(tumor, normal, by = "cases.submitter_id") %>%
    dplyr::select(cases.submitter_id)
  tumor <- tumor %>%
    dplyr::filter(cases.submitter_id %in% submitter_ids$cases.submitter_id)
  normal <- normal %>%
    dplyr::filter(cases.submitter_id %in% submitter_ids$cases.submitter_id)

  samples <- rbind(tumor, normal)
  unique(samples$sample_type)

  query_project <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = c("Solid Tissue Normal", "Primary Tumor"),
    barcode = as.list(samples$sample.submitter_id)
  )

  # If this is your first time running this notebook (i.e., you have not yet downloaded the results of the query in the previous block),
  # uncomment the line below

  # GDCdownload(query_project)

  return(list(samples = samples, query_project = query_project))
}

Download the TCGA gene expression data for colorectal cancer (TCGA-COAD).

projects <- c("TCGA-COAD")

with_results_projects <- c()

samples <- list()
project_data <- list()

for (project in projects) {
  result <- tryCatch(
    {
      result <- query_and_filter_samples(project)
      samples[[project]] <- result$samples
      project_data[[project]] <- result$query_project

      with_results_projects <- c(with_results_projects, project)
    },
    error = function(e) {

    }
  )
}

Running the code block above should generate and populate a directory named GDCdata.

III. Data preprocessing

Construct the RNA-seq count matrix for each cancer type.

tcga_data <- list()
tcga_matrix <- list()

projects <- with_results_projects
for (project in projects) {
  tcga_data[[project]] <- GDCprepare(project_data[[project]], summarizedExperiment = TRUE)
}

for (project in projects) {
  count_matrix <- assay(tcga_data[[project]], "unstranded")

  # Remove duplicate entries
  count_matrix_df <- data.frame(count_matrix)
  count_matrix_df <- count_matrix_df[!duplicated(count_matrix_df), ]
  count_matrix <- data.matrix(count_matrix_df)
  rownames(count_matrix) <- cleanid(rownames(count_matrix))
  count_matrix <- count_matrix[!(duplicated(rownames(count_matrix)) | duplicated(rownames(count_matrix), fromLast = TRUE)), ]

  tcga_matrix[[project]] <- count_matrix
}

Format the samples table so that it can be fed as input to DESeq2.

for (project in projects) {
  rownames(samples[[project]]) <- samples[[project]]$cases
  samples[[project]] <- samples[[project]] %>%
    dplyr::select(case = "cases.submitter_id", type = "sample_type")
  samples[[project]]$type <- str_replace(samples[[project]]$type, "Solid Tissue Normal", "normal")
  samples[[project]]$type <- str_replace(samples[[project]]$type, "Primary Tumor", "tumor")
}

DESeq2 requires the row names of samples should be identical to the column names of count_matrix.

for (project in projects) {
  colnames(tcga_matrix[[project]]) <- gsub(x = colnames(tcga_matrix[[project]]), pattern = "\\.", replacement = "-")
  tcga_matrix[[project]] <- tcga_matrix[[project]][, rownames(samples[[project]])]

  # Sanity check
  print(all(colnames(tcga_matrix[[project]]) == rownames(samples[[project]])))
}

IV. Differential gene expression analysis

For more detailed documentation on obtaining the gene set, refer to 7. Differential Gene Expression Analysis - TCGA - Pan-cancer - Unique Genes.Rmd.

RCDdb <- "temp/unique_genes/necroptosis_ferroptosis_pyroptosis/"

Write utility functions for filtering the gene sets, performing differential gene expression analysis, plotting the results, and performing variance-stabilizing transformation.

filter_gene_set_and_perform_dgea <- function(genes) {
  tcga_rcd <- list()

  for (project in projects) {
    rownames(genes) <- genes$gene_id
    tcga_rcd[[project]] <- tcga_matrix[[project]][rownames(tcga_matrix[[project]]) %in% genes$gene_id, ]
    tcga_rcd[[project]] <- tcga_rcd[[project]][, rownames(samples[[project]])]
  }

  dds_rcd <- list()
  res_rcd <- list()

  for (project in projects) {
    print(project)
    print("=============")
    dds <- DESeqDataSetFromMatrix(
      countData = tcga_rcd[[project]],
      colData = samples[[project]],
      design = ~type
    )
    dds <- filter_genes(dds, min_count = 10)
    dds$type <- relevel(dds$type, ref = "normal")
    dds_rcd[[project]] <- DESeq(dds)
    res_rcd[[project]] <- results(dds_rcd[[project]])
  }

  deseq.bbl.data <- list()

  for (project in projects) {
    deseq.results <- res_rcd[[project]]
    deseq.bbl.data[[project]] <- data.frame(
      row.names = rownames(deseq.results),
      baseMean = deseq.results$baseMean,
      log2FoldChange = deseq.results$log2FoldChange,
      lfcSE = deseq.results$lfcSE,
      stat = deseq.results$stat,
      pvalue = deseq.results$pvalue,
      padj = deseq.results$padj,
      cancer_type = project,
      gene_symbol = genes[rownames(deseq.results), "gene"]
    )
  }

  deseq.bbl.data.combined <- bind_rows(deseq.bbl.data)
  deseq.bbl.data.combined <- dplyr::filter(deseq.bbl.data.combined, abs(log2FoldChange) >= 1.5 & padj < 0.05)

  return(deseq.bbl.data.combined)
}

plot_dgea <- function(deseq.bbl.data.combined) {
  sizes <- c("<10^-15" = 4, "10^-10" = 3, "10^-5" = 2, "0.05" = 1)

  deseq.bbl.data.combined <- deseq.bbl.data.combined %>%
    mutate(fdr_category = cut(padj,
      breaks = c(-Inf, 1e-15, 1e-10, 1e-5, 0.05),
      labels = c("<10^-15", "10^-10", "10^-5", "0.05"),
      right = FALSE
    ))

  top_genes <- deseq.bbl.data.combined %>%
    group_by(cancer_type) %>%
    mutate(rank = rank(-abs(log2FoldChange))) %>%
    dplyr::filter(rank <= 10) %>%
    ungroup()

  ggplot(top_genes, aes(y = cancer_type, x = gene_symbol, size = fdr_category, fill = log2FoldChange)) +
    geom_point(alpha = 0.5, shape = 21, color = "black") +
    scale_size_manual(values = sizes) +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", limits = c(min(deseq.bbl.data.combined$log2FoldChange), max(deseq.bbl.data.combined$log2FoldChange))) +
    theme_minimal() +
    theme(
      axis.text.x = element_text(size = 9, angle = 90, hjust = 1)
    ) +
    theme(legend.position = "bottom") +
    theme(legend.position = "bottom") +
    labs(size = "Adjusted p-value", fill = "log2 FC", y = "Cancer type", x = "Gene")
}

perform_vsd <- function(genes) {
  tcga_rcd <- list()

  for (project in projects) {
    rownames(genes) <- genes$gene_id
    tcga_rcd[[project]] <- tcga_matrix[[project]][rownames(tcga_matrix[[project]]) %in% genes$gene_id, ]
    tcga_rcd[[project]] <- tcga_rcd[[project]][, rownames(samples[[project]])]
  }

  vsd_rcd <- list()

  for (project in projects) {
    print(project)
    print("=============")
    dds <- DESeqDataSetFromMatrix(
      countData = tcga_rcd[[project]],
      colData = samples[[project]],
      design = ~type
    )
    dds <- filter_genes(dds, min_count = 10)

    # Perform variance stabilization
    dds <- estimateSizeFactors(dds)
    nsub <- sum(rowMeans(counts(dds, normalized = TRUE)) > 10)
    vsd <- vst(dds, nsub = nsub)
    vsd_rcd[[project]] <- assay(vsd)
  }

  return(vsd_rcd)
}

Pyroptosis

Fetch the gene set of interest.

genes <- read.csv(paste0(RCDdb, "Pyroptosis.csv"))
print(genes)

genes$gene_id <- cleanid(genes$gene_id)
genes <- distinct(genes, gene_id, .keep_all = TRUE)
genes <- subset(genes, gene_id != "")
genes

Filter the genes to include only those in the gene set of interest, and then perform differential gene expression analysis.

deseq.bbl.data.combined <- filter_gene_set_and_perform_dgea(genes)

[1] "TCGA-COAD"
[1] "============="

Warning: some variables in design formula are characters, converting to factorsestimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

deseq.bbl.data.combined

Plot the results.

plot_dgea(deseq.bbl.data.combined)

Perform variance-stabilizing transformation for further downstream analysis (i.e., for survival analysis).

vsd <- perform_vsd(genes)

[1] "TCGA-COAD"
[1] "============="

V. Downloading the clinical data

Download clinical data from TCGA, and perform some preprocessing: - The deceased column should be FALSE if the patient is alive and TRUE otherwise - The overall_survival column should reflect the follow-up time if the patient is alive and the days to death otherwise

download_clinical_data <- function(project) {
  clinical_data <- GDCquery_clinic(project)
  clinical_data$deceased <- ifelse(clinical_data$vital_status == "Alive", FALSE, TRUE)
  clinical_data$overall_survival <- ifelse(clinical_data$vital_status == "Alive",
    clinical_data$days_to_last_follow_up,
    clinical_data$days_to_death
  )

  return(clinical_data)
}

tcga_clinical <- list()
for (project in projects) {
  tcga_clinical[[project]] <- download_clinical_data(project)
}

VI. Performing survival analysis

Write utility functions for performing survival analysis.

construct_gene_df <- function(gene_of_interest, project) {
  normal_df <- tcga_matrix[[project]] %>%
    as.data.frame() %>%
    rownames_to_column(var = "gene_id") %>%
    gather(key = "case_id", value = "counts", -gene_id) %>%
    left_join(., genes, by = "gene_id") %>%
    dplyr::filter(gene == gene_of_interest) %>%
    dplyr::filter(case_id %in% rownames(samples[[project]] %>% dplyr::filter(type == "normal")))
  normal_df$case_id <- paste0(sapply(strsplit(as.character(normal_df$case_id), "-"), `[`, 1), '-',
                            sapply(strsplit(as.character(normal_df$case_id), "-"), `[`, 2), '-', 
                            sapply(strsplit(as.character(normal_df$case_id), "-"), `[`, 3))
  
  tumor_df <- tcga_matrix[[project]] %>%
      as.data.frame() %>%
      rownames_to_column(var = "gene_id") %>%
      gather(key = "case_id", value = "counts", -gene_id) %>%
      left_join(., genes, by = "gene_id") %>%
      dplyr::filter(gene == gene_of_interest) %>%
      dplyr::filter(case_id %in% rownames(samples[[project]] %>% dplyr::filter(type == "tumor")))
  tumor_df$case_id <- paste0(sapply(strsplit(as.character(tumor_df$case_id), "-"), `[`, 1), '-',
                            sapply(strsplit(as.character(tumor_df$case_id), "-"), `[`, 2), '-', 
                            sapply(strsplit(as.character(tumor_df$case_id), "-"), `[`, 3))
  
  gene_df <- inner_join(normal_df, tumor_df, by = c("gene_id", "case_id", "deathtype", "gene", "description", "gene_biotype", "pmid", "comment"))
  gene_df$log_fold = log2(gene_df$counts.y / gene_df$counts.x)
  
  gene_df$strata <- ifelse(abs(gene_df$log_fold) >= 1.5, "HIGH", "LOW")
  gene_df <- merge(gene_df, tcga_clinical[[project]], by.x = "case_id", by.y = "submitter_id")
  
  return(gene_df)
}

compute_surival_fit <- function(gene_df) {
  return (survfit(Surv(overall_survival, deceased) ~ strata, data = gene_df))
}

compute_cox <- function(gene_df) {
  return (coxph(Surv(overall_survival, deceased) ~ strata, data=gene_df))
}

plot_survival <- function(fit) {
  return(ggsurvplot(fit,
    data = gene_df,
    pval = T,
    risk.table = T,
    risk.table.height = 0.3
  ))
}

compute_survival_diff <- function(gene_df) {
  return(survdiff(Surv(overall_survival, deceased) ~ strata, data = gene_df))
}

Perform survival analysis by testing for the difference in the Kaplan-Meier curves using the G-rho family of Harrington and Fleming tests: https://rdrr.io/cran/survival/man/survdiff.html

Our genes of interest are GSDMD (the primary executor of pyroptosis) and the differentially expressed genes..

significant_projects <- c()
significant_genes <- c()

ctr <- 1
for (project in projects) {
  for (gene in c("GSDMD", genes$gene)) {
    cat(project, gene, "\n\n")
    error <- tryCatch (
      {
        gene_df <- construct_gene_df(gene, project)
      },
      error = function(e) {
        cat("\n\n============================\n\n")
        e
      }
    )
    
    if(inherits(error, "error")) next

    if (nrow(gene_df) > 0) {
      fit <- compute_surival_fit(gene_df)
      tryCatch (
        {
          survival <- compute_survival_diff(gene_df)
          cox <- compute_cox(gene_df)
          print(ctr)
          ctr <- ctr + 1
          print(survival)
          cat("\n")
          print(cox)
          print(plot_survival(fit))
          if (pchisq(survival$chisq, length(survival$n)-1, lower.tail = FALSE) < 0.05) {
            significant_projects <- c(significant_projects, project)
            significant_genes <- c(significant_genes, gene)
          }
        },
        error = function(e) {
        }
      )
      
    }
    
    cat("\n\n============================\n\n")
  }
}

TCGA-COAD GSDMD 

[1] 1
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 12        3     3.16   0.00854    0.0125
strata=LOW  34        9     8.84   0.00306    0.0125

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW 0.07733   1.08040  0.69102 0.112 0.911

Likelihood ratio test=0.01  on 1 df, p=0.9104
n= 46, number of events= 12 


============================

TCGA-COAD CHMP7 

[1] 2
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 16        3     4.59     0.549     0.971
strata=LOW  30        9     7.41     0.340     0.971

 Chisq= 1  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.6636    1.9417   0.6849 0.969 0.333

Likelihood ratio test=1.02  on 1 df, p=0.3128
n= 46, number of events= 12 


============================

TCGA-COAD GSDMC 

[1] 3
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 32        8     8.06  0.000383   0.00145
strata=LOW  14        4     3.94  0.000782   0.00145

 Chisq= 0  on 1 degrees of freedom, p= 1 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z    p
strataLOW 0.02602   1.02636  0.68250 0.038 0.97

Likelihood ratio test=0  on 1 df, p=0.9696
n= 46, number of events= 12 


============================

TCGA-COAD ELANE 

[1] 4
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 44       11   11.428     0.016     0.343
strata=LOW   2        1    0.572     0.321     0.343

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.6125    1.8451   1.0623 0.577 0.564

Likelihood ratio test=0.28  on 1 df, p=0.5939
n= 46, number of events= 12 


============================

TCGA-COAD IRF1 

[1] 5
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 18        6     5.16     0.136     0.276
strata=LOW  28        6     6.84     0.103     0.276

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3196    0.7265   0.6109 -0.523 0.601

Likelihood ratio test=0.27  on 1 df, p=0.6032
n= 46, number of events= 12 


============================

TCGA-COAD CYCS 

[1] 6
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 19        6     4.49     0.509     0.937
strata=LOW  27        6     7.51     0.304     0.937

 Chisq= 0.9  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z    p
strataLOW -0.6206    0.5376   0.6510 -0.953 0.34

Likelihood ratio test=0.93  on 1 df, p=0.3342
n= 46, number of events= 12 


============================

TCGA-COAD GSDMA 

[1] 7
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 26        7        8     0.125     0.432
strata=LOW  20        5        4     0.250     0.432

 Chisq= 0.4  on 1 degrees of freedom, p= 0.5 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.4221    1.5252   0.6467 0.653 0.514

Likelihood ratio test=0.42  on 1 df, p=0.5147
n= 46, number of events= 12 


============================

TCGA-COAD CASP4 

[1] 8
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 11        4     3.99  3.47e-05  6.53e-05
strata=LOW  35        8     8.01  1.73e-05  6.53e-05

 Chisq= 0  on 1 degrees of freedom, p= 1 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

               coef exp(coef)  se(coef)      z     p
strataLOW -0.005542  0.994474  0.685863 -0.008 0.994

Likelihood ratio test=0  on 1 df, p=0.9936
n= 46, number of events= 12 


============================

TCGA-COAD BAK1 

[1] 9
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 20        4     4.92     0.171     0.345
strata=LOW  26        8     7.08     0.119     0.345

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.3811    1.4640   0.6526 0.584 0.559

Likelihood ratio test=0.35  on 1 df, p=0.5556
n= 46, number of events= 12 


============================

TCGA-COAD NOD1 

[1] 10
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH  9        3     2.01    0.4895     0.617
strata=LOW  37        9     9.99    0.0984     0.617

 Chisq= 0.6  on 1 degrees of freedom, p= 0.4 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.5387    0.5835   0.6938 -0.776 0.438

Likelihood ratio test=0.56  on 1 df, p=0.4552
n= 46, number of events= 12 


============================

TCGA-COAD NLRP7 

[1] 11
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 37       11    10.66    0.0107     0.101
strata=LOW   9        1     1.34    0.0852     0.101

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3377    0.7134   1.0657 -0.317 0.751

Likelihood ratio test=0.11  on 1 df, p=0.7409
n= 46, number of events= 12 


============================

TCGA-COAD CASP3 

[1] 12
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 15        3     3.44    0.0563    0.0878
strata=LOW  31        9     8.56    0.0226    0.0878

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2062    1.2291   0.6972 0.296 0.767

Likelihood ratio test=0.09  on 1 df, p=0.7646
n= 46, number of events= 12 


============================

TCGA-COAD GSDMB 

[1] 13
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 25        6     6.34    0.0187    0.0448
strata=LOW  21        6     5.66    0.0209    0.0448

 Chisq= 0  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.1292    1.1379   0.6106 0.212 0.832

Likelihood ratio test=0.04  on 1 df, p=0.8328
n= 46, number of events= 12 


============================

TCGA-COAD GZMB 

[1] 14
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 30        6     7.16     0.188     0.497
strata=LOW  16        6     4.84     0.277     0.497

 Chisq= 0.5  on 1 degrees of freedom, p= 0.5 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)   z     p
strataLOW 0.4243    1.5286   0.6064 0.7 0.484

Likelihood ratio test=0.49  on 1 df, p=0.4829
n= 46, number of events= 12 


============================

TCGA-COAD GSDME 

[1] 15
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 17        4     4.73    0.1141     0.201
strata=LOW  29        8     7.27    0.0744     0.201

 Chisq= 0.2  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2803    1.3235   0.6275 0.447 0.655

Likelihood ratio test=0.2  on 1 df, p=0.6515
n= 46, number of events= 12 


============================

TCGA-COAD CHMP3 

[1] 16
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 13        3     2.66    0.0438    0.0602
strata=LOW  33        9     9.34    0.0125    0.0602

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.1702    0.8435   0.6941 -0.245 0.806

Likelihood ratio test=0.06  on 1 df, p=0.8085
n= 46, number of events= 12 


============================

TCGA-COAD DPP9 

[1] 17
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH  9        1     1.73    0.3108     0.377
strata=LOW  37       11    10.27    0.0525     0.377

 Chisq= 0.4  on 1 degrees of freedom, p= 0.5 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.6365    1.8899   1.0547 0.604 0.546

Likelihood ratio test=0.43  on 1 df, p=0.5114
n= 46, number of events= 12 


============================

TCGA-COAD NOD2 

[1] 18
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 24        5     6.26     0.254     0.583
strata=LOW  22        7     5.74     0.277     0.583

 Chisq= 0.6  on 1 degrees of freedom, p= 0.4 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.4755    1.6089   0.6288 0.756 0.449

Likelihood ratio test=0.59  on 1 df, p=0.4417
n= 46, number of events= 12 


============================

TCGA-COAD NLRC4 

[1] 19
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 17        5     5.24   0.01107    0.0223
strata=LOW  29        7     6.76   0.00858    0.0223

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW 0.09364   1.09817  0.62729 0.149 0.881

Likelihood ratio test=0.02  on 1 df, p=0.8809
n= 46, number of events= 12 


============================

TCGA-COAD GSDMD 

[1] 20
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 12        3     3.16   0.00854    0.0125
strata=LOW  34        9     8.84   0.00306    0.0125

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW 0.07733   1.08040  0.69102 0.112 0.911

Likelihood ratio test=0.01  on 1 df, p=0.9104
n= 46, number of events= 12 


============================

TCGA-COAD TIRAP 

[1] 21
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 10        4     2.44     0.999      1.33
strata=LOW  36        8     9.56     0.255      1.33

 Chisq= 1.3  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.7315    0.4812   0.6484 -1.128 0.259

Likelihood ratio test=1.19  on 1 df, p=0.2751
n= 46, number of events= 12 


============================

TCGA-COAD SCAF11 

[1] 22
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 15        3     3.48    0.0674     0.104
strata=LOW  31        9     8.52    0.0276     0.104

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2228    1.2496   0.6920 0.322 0.747

Likelihood ratio test=0.11  on 1 df, p=0.744
n= 46, number of events= 12 


============================

TCGA-COAD NLRP6 

[1] 23
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 28        9     7.21     0.446       1.2
strata=LOW  18        3     4.79     0.671       1.2

 Chisq= 1.2  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.7302    0.4818   0.6805 -1.073 0.283

Likelihood ratio test=1.26  on 1 df, p=0.2611
n= 46, number of events= 12 


============================

TCGA-COAD AIM2 

[1] 24
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 25        7     5.14     0.672      1.33
strata=LOW  21        5     6.86     0.504      1.33

 Chisq= 1.3  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW -0.7302    0.4818   0.6462 -1.13 0.258

Likelihood ratio test=1.34  on 1 df, p=0.2474
n= 46, number of events= 12 


============================

TCGA-COAD CASP6 

[1] 25
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 17        3     3.55    0.0846     0.134
strata=LOW  29        9     8.45    0.0355     0.134

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2549    1.2903   0.6972 0.366 0.715

Likelihood ratio test=0.14  on 1 df, p=0.7105
n= 46, number of events= 12 


============================

TCGA-COAD NLRP2 

[1] 26
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 38       12     9.86     0.463      2.69
strata=LOW   8        0     2.14     2.138      2.69

 Chisq= 2.7  on 1 degrees of freedom, p= 0.1 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

                coef  exp(coef)   se(coef)      z     p
strataLOW -1.937e+01  3.862e-09  9.797e+03 -0.002 0.998

Likelihood ratio test=4.79  on 1 df, p=0.02864
n= 46, number of events= 12 


============================

TCGA-COAD IRF2 

[1] 27
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 17        3     3.82    0.1759     0.293
strata=LOW  29        9     8.18    0.0821     0.293

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z    p
strataLOW 0.3758    1.4562   0.6980 0.538 0.59

Likelihood ratio test=0.3  on 1 df, p=0.5824
n= 46, number of events= 12 


============================

TCGA-COAD PJVK 

[1] 28
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH  8        3     2.27    0.2377     0.301
strata=LOW  38        9     9.73    0.0553     0.301

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3706    0.6903   0.6792 -0.546 0.585

Likelihood ratio test=0.28  on 1 df, p=0.5963
n= 46, number of events= 12 


============================

TCGA-COAD CASP5 

[1] 29
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 36        9     9.39    0.0158     0.083
strata=LOW  10        3     2.61    0.0567     0.083

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2059    1.2287   0.7155 0.288 0.773

Likelihood ratio test=0.08  on 1 df, p=0.7766
n= 46, number of events= 12 


============================

TCGA-COAD NLRP1 

[1] 30
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 22        8     5.12      1.62      3.14
strata=LOW  24        4     6.88      1.21      3.14

 Chisq= 3.1  on 1 degrees of freedom, p= 0.08 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z      p
strataLOW -1.1687    0.3108   0.6930 -1.687 0.0917

Likelihood ratio test=3.23  on 1 df, p=0.07214
n= 46, number of events= 12 


============================

TCGA-COAD CASP9 

[1] 31
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 15        3     3.19    0.0108    0.0161
strata=LOW  31        9     8.81    0.0039    0.0161

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW 0.08831   1.09232  0.69600 0.127 0.899

Likelihood ratio test=0.02  on 1 df, p=0.8985
n= 46, number of events= 12 


============================

TCGA-COAD PLCG1 

[1] 32
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 16        3     3.58    0.0934      0.14
strata=LOW  30        9     8.42    0.0397      0.14

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.2538    1.2889   0.6803 0.373 0.709

Likelihood ratio test=0.14  on 1 df, p=0.7041
n= 46, number of events= 12 


============================

TCGA-COAD IL18 

[1] 33
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 28        6     6.43    0.0286    0.0702
strata=LOW  18        6     5.57    0.0331    0.0702

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.1617    1.1755   0.6109 0.265 0.791

Likelihood ratio test=0.07  on 1 df, p=0.7919
n= 46, number of events= 12 


============================

TCGA-COAD DPP8 

[1] 34
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

             N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 16        4     3.95   0.00057  0.000906
strata=LOW  30        8     8.05   0.00028  0.000906

 Chisq= 0  on 1 degrees of freedom, p= 1 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

              coef exp(coef) se(coef)     z     p
strataLOW -0.01903   0.98115  0.63232 -0.03 0.976

Likelihood ratio test=0  on 1 df, p=0.976
n= 46, number of events= 12 


============================

Display the results only for genes where a significant difference in survival has been reported.

significant_genes

NULL

num_significant_genes <- length(significant_genes)

if (num_significant_genes > 0) {
  for (i in 1 : num_significant_genes) {
    project <- significant_projects[[i]]
    gene <- significant_genes[[i]]
    
    cat(project, gene, "\n\n")
    gene_df <- construct_gene_df(gene, project)
    
    survival <- compute_survival_diff(gene_df)
    cox <- compute_cox(gene_df)
    print(survival)
    cat("\n")
    print(cox)
    print(plot_survival(fit))
    
    cat("\n\n============================\n\n")
  } 
}

---

1. De La Salle University, Manila, Philippines, gonzales.markedward@gmail.com↩︎

2. De La Salle University, Manila, Philippines, anish.shrestha@dlsu.edu.ph↩︎