Here we are demonstrating the feasibility of analyzing genomic data using Stan. The first use case is to analyze somatic mutations for association with survival, after adjusting for key clinical variables with known prognostic status.
knitr::opts_chunk$set(fig.width=8, fig.height=6, fig.path='Figs/', echo=TRUE, warning=FALSE, message=FALSE) # library(SuMu) devtools::load_all('.') library(dplyr) library(survminer) library(survival) library(ggplot2) library(scales) library(tidyr) library(rstanarm) options(mc.cores = 4)
First, download the clinical data. Here we are using the TCGA skin cutaneous melanoma (SKCM) cohort.
clin_df <- SuMu::get_tcga_clinical(cohort = "SKCM") ## format some clinical data variables clin_df2 <- clin_df %>% dplyr::mutate(stage_part1 = gsub(pathologic_stage, pattern = '(Stage [0I]+).*', replacement = '\\1'), diagnosis_year_group = cut(year_of_initial_pathologic_diagnosis, breaks = c(1975, 1990, 1995, 2000, 2005, 2010, 2015, 2020), include.lowest = TRUE), os_10y = ifelse(OS_IND == 1 & OS <= 10*365.25, 1, 0), sample = sampleID )
For this analysis we will consider the survival time in r print(unique(clin_df['OS_UNIT']))
since initial pathologic diagnosis.
fit <- survfit(Surv(OS, OS_IND) ~ 1, data = clin_df2) survminer::ggsurvplot(fit) + ggtitle('Survival since diagnosis in full cohort')
Plotting by stage, although the time of 'stage' determination may be confounded if not collected at time of initial diagnosis.
fit <- survfit(Surv(OS, OS_IND) ~ pathologic_stage, data = clin_df2) survminer::ggsurvplot(fit, legend = "right")
There also seem to be differences by tumor type.
fit <- survfit(Surv(OS, OS_IND) ~ sample_type, data = clin_df2) survminer::ggsurvplot(fit, legend = "right")
(Aside: I wonder how similar tumor type is to sample type? For example, we could have a metastatic patient where the sample was obtained from the primary tumor. We will want to adjust our genetic data analysis for the sample type but may want to estimate prognosis according to the tumor type?)
A variable like year_of_initial_pathologic_diagnosis
is guaranteed to be unconfounded since we can safely assume it was collected at the time of diagnosis.
fit <- survfit(Surv(OS, OS_IND) ~ diagnosis_year_group, data = clin_df2) survminer::ggsurvplot(fit, legend = 'right')
This makes it pretty clear that we have a strong "survival" bias to our data. This would suggest that, among people whose diagnosis was made in the 90s, only those who survived long enough to be enrolled were included in the study.
Let's look at a histogram of years of initial diagnosis.
ggplot(clin_df2, aes(x = year_of_initial_pathologic_diagnosis, fill = diagnosis_year_group)) + geom_histogram() + theme_minimal()
Let's look at the time since initial diagnosis (presumably, the time from enrollment to diagnosis).
Finally, we can visualize a more comprehesive set of clinical variables.
fit <- survival::coxph(Surv(OS, OS_IND) ~ age_at_initial_pathologic_diagnosis + sample_type + breslow_depth_value + initial_weight + strata(year_of_initial_pathologic_diagnosis), data = clin_df2) print(fit)
We can download the somatic mutations to supplement the phenotypes.
mut_df <- SuMu::get_tcga_somatic_mutations(cohort = "SKCM") %>% dplyr::mutate(gene_aa = paste0(gene, ".", Amino_Acid_Change), gene_effect = paste0(gene, ".", effect) )
Check the most frequent mutations.
mut_df_missense = mut_df %>% dplyr::filter(effect == "Missense_Mutation") mut_df_missense$gene_aa = paste0(mut_df_missense$gene, ":", mut_df_missense$Amino_Acid_Change) mut_df_missense %>% select(gene_aa) %>% table %>% sort %>% rev %>% as.data.frame %>% head(10)
Filter to top genes
top_genes <- mut_df %>% dplyr::group_by(gene) %>% dplyr::mutate(gene_count = n()) %>% dplyr::ungroup() %>% dplyr::distinct(gene, .keep_all = TRUE) %>% dplyr::top_n(gene_count, n = 10) %>% dplyr::select(gene) mut_df_topgenes <- mut_df %>% dplyr::semi_join(top_genes)
clin_df2_nonmiss <- clin_df2 %>% dplyr::mutate( revised_breslow_depth = ifelse(is.na(breslow_depth_value), 0, breslow_depth_value)) %>% tidyr::drop_na(os_10y, age_at_initial_pathologic_diagnosis, initial_weight, revised_breslow_depth, sample_type, diagnosis_year_group) mutation_matrix <- SuMu:::prep_biomarker_data( biomarker_data = mut_df_topgenes, data = clin_df2_nonmiss, biomarker_formula = 1 ~ gene_effect, .fun = sum, id = 'sample' ) glm_df2 <- mutation_matrix %>% dplyr::left_join(clin_df2_nonmiss %>% dplyr::select(sample, os_10y), by = 'sample')
stan-glm
model to these genetic data# construct input formula gene_names2 <- mutation_matrix %>% head(1) %>% dplyr::select(-sample) %>% names() gene_subformula <- stringr::str_c('`', stringr::str_c(gene_names2, collapse = '` + `'), '`') my_formula2 <- stringr::str_c('os_10y ~ ', gene_subformula) # call to `stan_glm` glmfit2 <- rstanarm::stan_glm( data = glm_df2, formula = my_formula2, sparse = TRUE, family = binomial(), chains = 4, prior = rstanarm::hs_plus() )
rescale <- function(x) { (x - mean(x, na.rm=T))/(2*sd(x, na.rm=T)) } clin_df3 <- clin_df2 %>% dplyr::mutate( rescale_age_at_initial_pathologic_diagnosis = rescale(age_at_initial_pathologic_diagnosis), rescale_initial_weight = rescale(initial_weight), rescale_breslow_depth_value = rescale(breslow_depth_value) ) glmfit_clin <- rstanarm::stan_glmer( os_10y ~ rescale_age_at_initial_pathologic_diagnosis + sample_type + rescale_breslow_depth_value + rescale_initial_weight + ( rescale_age_at_initial_pathologic_diagnosis + sample_type + rescale_breslow_depth_value + rescale_initial_weight | diagnosis_year_group ), data = clin_df3, init_r = 1, family = binomial() ) print(glmfit_clin)
glm_df3 <- clin_df3 %>% dplyr::inner_join(mutation_matrix, by = 'sample') # construct input formula clinical_formula <- os_10y ~ rescale_age_at_initial_pathologic_diagnosis + sample_type + rescale_breslow_depth_value + rescale_initial_weight + `__BIOMARKERS__` + ( rescale_age_at_initial_pathologic_diagnosis + sample_type + rescale_breslow_depth_value + rescale_initial_weight + `__BIOMARKERS__` | diagnosis_year_group ) gene_subformula <- stringr::str_c('`', stringr::str_c(gene_names2, collapse = '` + `'), '`') my_formula3 <- stringr::str_c( as.character(clinical_formula)[2], as.character(clinical_formula)[3], sep = as.character(clinical_formula)[1]) my_formula3 <- as.formula(gsub(my_formula3, pattern = '`__BIOMARKERS__`', replacement = gene_subformula)) update(clinical_formula, stringr::str_c('~ . ', gene_subformula, stringr::str_c('(', gene_subformula, '| diagnosis_year_group)'), sep = '+') ) # call to `stan_glm` glmfit_clingen <- rstanarm::stan_glmer( data = glm_df3, formula = my_formula3, sparse = TRUE, family = binomial(), chains = 4, prior = rstanarm::hs_plus() )
summary_table=feature_table(glmfit_clingen) feature_graph(glmfit_clingen) view_feature(mutation_matrix,clin_df2,gsub("`","",rownames(summary_table)[1])) view_feature(mutation_matrix,clin_df2,gsub("`","",rownames(summary_table)[nrow(summary_table)]))
auc(clin_df3 %>% dplyr::select(os_10y, rescale_age_at_initial_pathologic_diagnosis, sample_type, rescale_breslow_depth_value, rescale_initial_weight, diagnosis_year_group) %>% tidyr::drop_na(), "os_10y", glmfit_clin, h_gram = T, roc_plot = T)
Our function syntax will look like the following.
fit <- fit_glm( data = clin_df, formula = os_10y ~ rescale_.. + `__BIO__`, biomarker_data = mut_df, biomarker_formula = 1 ~ (1|gene_aa) + (1|gene) + (1|effect), id = 'sample' )
It currently supports the following
fit <- fit_glm( data = clin_df, formula = os_10y ~ age + weight + tumor_type, biomarker_data = mut_df, biomarker_formula = DNA_VAF ~ gene_aa, id = 'sample' )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.