03. Exploratory Data Analysis

Overview

See the Glossary for term definitions used throughout this project.

  • Patient demographics, biospecimen characteristics, and RNA-seq quality metrics
  • DuckDB queries against underlying parquet files
  • Cytogenetic alteration frequencies and co-occurrence patterns

All summaries are pre-computed as targets (vig_*). Placeholder messages appear in place of plots when pre-computed outputs are not yet available.

Clinical Demographics

  • Patients: 995
  • Variables: 88
  • Data completeness: 95.8%

Age at Diagnosis

Generating code 
{
    if (is.null(eda_clinical_summary)) 
        return(NULL)
    clin <- eda_clinical_summary
    if (is.null(clin$age_years) || length(clin$age_years) == 
        0) 
        return(NULL)
    age_df <- data.frame(age = clin$age_years)
    med_age <- sprintf("%.1f", median(clin$age_years, na.rm = TRUE))
    ggplot2::ggplot(age_df, ggplot2::aes(x = age)) + ggplot2::geom_histogram(bins = 30, 
        fill = "steelblue", colour = "white") + ggplot2::geom_vline(xintercept = median(clin$age_years, 
        na.rm = TRUE), linetype = "dashed", colour = "red") + 
        ggplot2::labs(x = "Age (years)", y = "Number of Patients", 
            title = "Age at Diagnosis", subtitle = paste0("n = ", 
                length(clin$age_years), " patients, median = ", 
                med_age, " years"), caption = paste0("Age at diagnosis for ", 
                length(clin$age_years), " CoMMpass patients. ", 
                "x-axis: age in years (converted from GDC days via / 365.25); ", 
                "y-axis: number of patients. ", "Median = ", 
                med_age, " years (dashed red line). ", "Source: GDC clinical. ", 
                "See glossary for age unit conventions; ", "survival-analysis Cox models include age as a covariate.")) + 
        theme_commpass_dark() + ggplot2::theme(plot.caption = ggplot2::element_text(size = 7, 
        hjust = 0, lineheight = 1.2))
}

Gender Distribution

Generating code 
{
    if (is.null(eda_clinical_summary)) 
        return(NULL)
    clin <- eda_clinical_summary
    if (is.null(clin$gender_table)) 
        return(NULL)
    gt <- clin$gender_table
    gt$Percentage <- sprintf("%.1f%%", 100 * gt$Freq/sum(gt$Freq))
    caption <- paste0("Gender distribution for ", clin$n_patients, 
        " CoMMpass patients. ", "Var1 = GDC gender category (male/female). ", 
        "Freq = number of patients. ", "Percentage = proportion of total cohort. ", 
        "Data: GDC clinical ", "(https://portal.gdc.cancer.gov/projects/MMRF-COMMPASS).")
    DT::datatable(gt, rownames = FALSE, filter = "top", options = list(pageLength = 10, 
        scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
        caption))
}

Race Distribution

Generating code 
{
    if (is.null(eda_clinical_summary)) 
        return(NULL)
    clin <- eda_clinical_summary
    if (is.null(clin$race_table)) 
        return(NULL)
    rt <- clin$race_table
    rt$Percentage <- sprintf("%.1f%%", 100 * rt$Freq/sum(rt$Freq))
    caption <- paste0("Race distribution for ", clin$n_patients, 
        " CoMMpass patients. ", "Var1 = self-reported race category (GDC definitions). ", 
        "Freq = number of patients. ", "Percentage = proportion of total cohort. ", 
        "Categories from GDC clinical data.")
    DT::datatable(rt, rownames = FALSE, filter = "top", options = list(pageLength = 10, 
        scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
        caption))
}

Vital Status

Generating code 
{
    if (is.null(eda_clinical_summary)) 
        return(NULL)
    clin <- eda_clinical_summary
    if (is.null(clin$vital_table)) 
        return(NULL)
    vt <- clin$vital_table
    vt$Percentage <- sprintf("%.1f%%", 100 * vt$Freq/sum(vt$Freq))
    caption <- paste0("Vital status for ", clin$n_patients, " CoMMpass patients. ", 
        "Status = Alive (censored at last follow-up) or Dead (event observed). ", 
        "Freq = number of patients. ", "Percentage = proportion of total cohort. ", 
        "Data: GDC clinical.")
    DT::datatable(vt, rownames = FALSE, filter = "top", options = list(pageLength = 10, 
        scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
        caption))
}

ISS Staging

The International Staging System (ISS) classifies multiple myeloma into three stages based on serum beta-2 microglobulin and albumin levels. Higher stages indicate more advanced disease and worse prognosis.

  • Patients with ISS data: 995
  • Missing ISS: 0
  • Source column: iss_stage

ISS Stage Distribution

Generating code 
{
    if (is.null(eda_iss_summary)) 
        return(NULL)
    iss <- eda_iss_summary
    if (!isTRUE(iss$available)) 
        return(NULL)
    caption <- paste0("ISS stage distribution for ", iss$n_with_iss, 
        " patients with staging data (", iss$n_missing, " missing). ", 
        "ISS = International Staging System ", "(https://doi.org/10.1200/JCO.2005.04.242). ", 
        "Stage I = beta-2 microglobulin < 3.5 mg/L and albumin >= 3.5 g/dL. ", 
        "Stage II = not I or III. ", "Stage III = beta-2 microglobulin >= 5.5 mg/L. ", 
        "Freq = number of patients in each stage.")
    DT::datatable(iss$iss_table, rownames = FALSE, filter = "top", 
        options = list(pageLength = 10, scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}
Generating code 
{
    if (is.null(eda_iss_summary)) 
        return(NULL)
    iss <- eda_iss_summary
    if (!isTRUE(iss$available) || nrow(iss$iss_table) == 0) 
        return(NULL)
    iss_plot <- iss$iss_table[!is.na(iss$iss_table$ISS_Stage), 
        ]
    if (nrow(iss_plot) == 0) 
        return(NULL)
    ggplot2::ggplot(iss_plot, ggplot2::aes(x = ISS_Stage, y = Freq)) + 
        ggplot2::geom_col(fill = "steelblue") + ggplot2::labs(x = "ISS Stage", 
        y = "Number of Patients", title = "ISS Stage Distribution", 
        subtitle = paste0("n = ", sum(iss_plot$Freq), " patients with staging data"), 
        caption = paste0("ISS stage distribution for ", sum(iss_plot$Freq), 
            " patients with staging data. ", "Stage I: B2M < 3.5 mg/L + albumin >= 3.5 g/dL. ", 
            "Stage II: not I or III. Stage III: B2M >= 5.5 mg/L. ", 
            "x-axis: ISS stage; y-axis: number of patients. ", 
            "Source: ISS (Greipp et al. 2005, doi:10.1200/JCO.2005.04.242). ", 
            "See survival-analysis KM curves by ISS for prognostic impact.")) + 
        theme_commpass_dark() + ggplot2::theme(plot.caption = ggplot2::element_text(size = 7, 
        hjust = 0, lineheight = 1.2))
}

Age by ISS Stage

Generating code 
{
    if (is.null(eda_iss_summary)) 
        return(NULL)
    iss <- eda_iss_summary
    if (!isTRUE(iss$available) || is.null(iss$age_by_iss)) 
        return(NULL)
    age_tbl <- do.call(rbind, lapply(names(iss$age_by_iss), function(stage) {
        x <- iss$age_by_iss[[stage]]
        data.frame(ISS_Stage = stage, N = x$n, Mean_Age = sprintf("%.1f", 
            x$mean), Median_Age = sprintf("%.1f", x$median), 
            SD = sprintf("%.1f", x$sd))
    }))
    caption <- paste0("Age at diagnosis (years) by ISS stage. ", 
        "ISS_Stage = International Staging System stage (I/II/III). ", 
        "N = number of patients. ", "Mean_Age, Median_Age = central tendency in years. ", 
        "SD = standard deviation. ", "GDC stores age in days; converted via age / 365.25.")
    DT::datatable(age_tbl, rownames = FALSE, filter = "top", 
        options = list(pageLength = 10, scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}

Survival Endpoints

Days to Death

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    8.0   209.5   475.0   525.5   723.5  1753.0
  • Censoring rate: 80.8% (patients still alive at last follow-up)

Days to Last Follow-up

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   -3.0   402.5   778.0   776.1  1085.5  1984.0

Biospecimen Overview

  • Total samples: 2,119
  • Variables: 31

Sample Types

Samples per Patient

RNA-seq Quality

  • Genes: 60,660
  • Samples: 100
  • Total counts: 6,045,546,169
  • Median library size: 58,030,298
  • Sparsity (% zero entries): 49.1%

Library Size Distribution

Genes Detected per Sample

Sample Clustering (PCA)

Principal component analysis on the top 500 most variable genes (VST-transformed) reveals sample-level structure. Points are colored by clinical covariates to check for batch effects or biological groupings.

Cytogenetic Landscape

Cytogenetic markers from FISH testing are key to myeloma risk stratification (IMWG 2014 criteria). High-risk alterations include t(4;14), t(14;16), del(17p), and gain(1q).

Alteration Frequencies

Generating code 
{
    if (is.null(cyto_frequency_summary) || nrow(cyto_frequency_summary) == 
        0) 
        return(NULL)
    caption <- paste0("Cytogenetic alteration frequencies from FISH testing. ", 
        "Marker = FISH probe target. ", "Positive = patients with alteration detected. ", 
        "Tested = patients with FISH data for that marker. ", 
        "% = prevalence among tested patients. ", "High-risk markers per IMWG 2014: t(4;14), t(14;16), del(17p), gain(1q). ", 
        "t() = translocation, del() = deletion, gain() = extra copy. ", 
        "See data dictionary for marker definitions.")
    DT::datatable(cyto_frequency_summary, rownames = FALSE, filter = "top", 
        options = list(pageLength = 15, scrollX = TRUE), colnames = c("Marker", 
            "Positive", "Tested", "%"), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}

Oncoprint

Co-occurrence Analysis

Generating code 
{
    if (is.null(cyto_cooccurrence) || nrow(cyto_cooccurrence) == 
        0) 
        return(NULL)
    marker_labels <- c(t_4_14 = "t(4;14)", t_11_14 = "t(11;14)", 
        t_14_16 = "t(14;16)", t_14_20 = "t(14;20)", del_17p = "del(17p)", 
        del_1p = "del(1p)", gain_1q = "gain(1q)")
    display <- cyto_cooccurrence
    display$marker1 <- ifelse(display$marker1 %in% names(marker_labels), 
        marker_labels[display$marker1], display$marker1)
    display$marker2 <- ifelse(display$marker2 %in% names(marker_labels), 
        marker_labels[display$marker2], display$marker2)
    display_df <- display[, c("marker1", "marker2", "odds_ratio", 
        "pvalue", "padj", "n_both", "tendency")]
    display_df$odds_ratio <- round(display_df$odds_ratio, 2)
    display_df$pvalue <- round(display_df$pvalue, 4)
    display_df$padj <- round(display_df$padj, 4)
    caption <- paste0("Pairwise co-occurrence analysis of cytogenetic markers ", 
        "(Fisher's exact test, BH-adjusted). ", "Marker 1/2 = FISH probe targets tested pairwise. ", 
        "Odds Ratio > 1 = co-occurrence; < 1 = mutual exclusivity. ", 
        "p-value = Fisher's exact test. ", "Adjusted p = Benjamini-Hochberg corrected. ", 
        "Both+ = patients positive for both markers. ", "Tendency = co-occurrence or mutual exclusivity. ", 
        "* = padj < 0.05.")
    DT::datatable(display_df, rownames = FALSE, filter = "top", 
        options = list(pageLength = 15, scrollX = TRUE), colnames = c("Marker 1", 
            "Marker 2", "Odds Ratio", "p-value", "Adjusted p", 
            "Both+", "Tendency"), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}

Cytogenetic markers are used for survival stratification and inform differential expression analysis.

DuckDB Query Examples

The package provides two functions for querying the parquet files:

  • query_commpass_parquet() – one-shot query that returns a data frame
  • get_commpass_tbl() – returns a lazy dplyr::tbl() for chained operations

Example: Simple Query

Show code 
safe_tar_read("code_eda_simple_query")

One-shot query (returns data frame)

clinical <- query_commpass_parquet(“clinical”) head(clinical)

With dplyr/tbl syntax

con <- DBI::dbConnect(duckdb::duckdb()) clinical_tbl <- get_commpass_tbl(“clinical”, con = con) females <- clinical_tbl |> dplyr::filter(gender == “female”) |> dplyr::collect() DBI::dbDisconnect(con, shutdown = TRUE)

Generating code 
{
    if (is.null(eda_duckdb_demo)) 
        return(NULL)
    demo <- eda_duckdb_demo
    caption <- paste0("First 10 patients from clinical data via dplyr/DuckDB query. ", 
        "age_at_diagnosis is in days (GDC convention; divide by 365.25 for years). ", 
        "This demonstrates zero-copy DuckDB access to parquet files.")
    DT::datatable(demo$demo_result, rownames = FALSE, filter = "top", 
        options = list(pageLength = 10, scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}

Example: Aggregation with DuckDB

Show code 
safe_tar_read("code_eda_aggregation")

Lazy tbl for dplyr-style queries

con <- DBI::dbConnect(duckdb::duckdb()) clinical_tbl <- get_commpass_tbl(“clinical”, con = con)

result <- clinical_tbl |> dplyr::filter(!is.na(gender)) |> dplyr::group_by(gender) |> dplyr::summarise( n = dplyr::n(), mean_age_years = mean(age_at_diagnosis / 365.25, na.rm = TRUE) ) |> dplyr::collect()

DBI::dbDisconnect(con, shutdown = TRUE)

Generating code 
{
    if (is.null(eda_duckdb_demo)) 
        return(NULL)
    demo <- eda_duckdb_demo
    caption <- paste0("Patient count and mean age (years) by gender via ", 
        "dplyr/DuckDB aggregation. ", "n = number of patients. ", 
        "mean_age = average age at diagnosis in years.")
    DT::datatable(demo$agg_result, rownames = FALSE, filter = "top", 
        options = list(pageLength = 10, scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            caption))
}

Cross-dataset Integration

Patient ID overlap between clinical and biospecimen datasets:

  • Clinical patients: 995
  • Biospecimen patients: 2,119
  • Shared (in both): 0
  • Clinical only: 995
  • Biospecimen only: 2,119

Clinical-Informed Normalisation

Standard z-scoring treats all variation equally, but for clinical biomarkers this compresses the clinically important normal range and expands the uninformative extreme range. The SCOPE method (Hussain et al. 2024) uses clinical reference ranges to normalise biomarkers, preserving clinically meaningful variation.

Reference Ranges

Clinical vs Z-Score Comparison

Data not available: Target vig_normalisation_comparison returned NULL. See Data Sources.

Missingness as Signal

Missing data in clinical datasets is rarely random. In multiple myeloma, FISH testing is expensive and may be selectively ordered for patients with suspected high-risk disease. If missingness correlates with survival, it suggests the data is missing not at random (MNAR) — the missingness itself is clinically informative.

Variable Missingness

Does Missingness Predict Survival?

Cox regression testing whether patients with missing values for each variable have different survival outcomes. A significant HR suggests the missingness mechanism is informative (MNAR).

Data Sources

Results in this vignette are derived from the MMRF CoMMpass study (MMRF-COMMPASS, ~1,143 patients), downloaded via TCGAbiolinks. The pipeline runs with a configurable sample_limit (default 200; CI uses 20).

For full citations, data access tiers, and the distinction between pipeline data and synthetic test data, see the Data Sources vignette.

Recent Changes

Recent project commits with lines added, files changed, and change categories.

Reproducibility

Git Commit Info (click to expand)
Generating code 
{
    if (!requireNamespace("gert", quietly = TRUE)) 
        return(NULL)
    tryCatch({
        info <- gert::git_info()
        log1 <- gert::git_log(max = 1)
        git_df <- data.frame(Item = c("Commit Hash", "Author", 
            "Time", "Branch"), Value = c(info$commit, log1$author, 
            as.character(log1$time), info$shorthand), stringsAsFactors = FALSE)
        DT::datatable(git_df, rownames = FALSE, options = list(pageLength = 10, 
            dom = "t", scrollX = TRUE), caption = htmltools::tags$caption(style = "caption-side: top; text-align: left;", 
            "Git repository state at pipeline build time."))
    }, error = function(e) NULL)
}
Session Info (click to expand)
Show code 
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] targets_1.12.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] base64url_1.4      gtable_0.3.6       jsonlite_2.0.0     dplyr_1.2.1       
#>  [5] compiler_4.6.1     tidyselect_1.2.1   jquerylib_0.1.4    callr_3.8.0       
#>  [9] scales_1.4.0       yaml_2.3.12        fastmap_1.2.0      ggplot2_4.0.3     
#> [13] R6_2.6.1           labeling_0.4.3     generics_0.1.4     igraph_2.3.3      
#> [17] knitr_1.51         htmlwidgets_1.6.4  backports_1.5.1    tibble_3.3.1      
#> [21] maketools_1.3.2    bslib_0.11.0       pillar_1.11.1      RColorBrewer_1.1-3
#> [25] rlang_1.3.0        DT_0.34.0          cachem_1.1.0       xfun_0.60         
#> [29] sass_0.4.10        sys_3.4.3          S7_0.2.2           otel_0.2.0        
#> [33] cli_3.6.6          withr_3.0.3        magrittr_2.0.5     crosstalk_1.2.2   
#> [37] ps_1.9.3           digest_0.6.39      grid_4.6.1         processx_3.9.0    
#> [41] secretbase_1.3.0   lifecycle_1.0.5    prettyunits_1.2.0  vctrs_0.7.3       
#> [45] evaluate_1.0.5     glue_1.8.1         data.table_1.18.4  farver_2.1.2      
#> [49] codetools_0.2-20   buildtools_1.0.0   rmarkdown_2.31     tools_4.6.1       
#> [53] pkgconfig_2.0.3    htmltools_0.5.9