Statistical errata and post-hoc tests

Setup

library(ithi.utils)
load_base_libs()

library(methods)
library(grid)
library(gridExtra)
library(gridBase)

library(ithi.meta)
library(ithi.figures)
library(ithi.utils)
library(ithi.seq)
library(ithi.clones)
library(ithi.supp)
library(ithi.xcr)

ihc_table_path <- snakemake@input$ihc_table
xcr_table_path <- snakemake@input$xcr_table
clone_tree_file <- snakemake@input$clone_tree_file
clone_branch_length_file <- snakemake@input$clone_branch_length_file
clone_prevalence_file <- snakemake@input$clone_prevalence_file
tcr_diversity_file <- snakemake@input$tcr_diversity
bcr_diversity_file <- snakemake@input$bcr_diversity
molecular_subtype_file <- snakemake@input$molsubtypes
clonal_measures_file <- snakemake@input$ith_stats_file
mutsig_dir <- snakemake@input$mutsig_dir
ith_icgc_bc_file <- snakemake@input$ith_icgc_bc
nanostring_annotations_path <- snakemake@input$nanostring_annotations
icgc_subtype_file <- snakemake@input$icgc_subtypes
icgc_specimen_file <- snakemake@input$icgc_specimen
master_variant_file <- snakemake@input$snv_table
master_breakpoint_file <- snakemake@input$breakpoint_table

db_path <- snakemake@params$db
tils_for_cluster <- snakemake@params$tils_for_cluster
stat_types <- snakemake@params$ith_stat_types

nclusts <- 3

annotation_colours <- ithi.figures::get_annotation_colours()

ihc_table <- fread(ihc_table_path)
xcr_table <- read_clonotypes(xcr_table_path, duplicates = FALSE, db_path = db_path)

Read 19.7% of 304822 rows
Read 62.3% of 304822 rows
Read 88.6% of 304822 rows
Read 304822 rows and 18 (of 18) columns from 0.070 GB file in 00:00:05

molsubtypes <- fread(molecular_subtype_file)

tree_branch_data <- read_clone_tree_data(clone_tree_file, clone_branch_length_file, 
    clone_prevalence_file, db_path)

xcr_diversity <- ithi.supp::get_xcr_diversity(tcr_diversity_file, bcr_diversity_file, 
    db_path, xcr_table)

clonal_measures <- ithi.clones::read_ith_stats(clonal_measures_file, db_path, 
    duplicates = FALSE)

ith_icgc_bc <- fread(ith_icgc_bc_file)
icgc_specimen_data <- fread(icgc_specimen_file)

nanostring_labels <- fread(nanostring_annotations_path)

til_clusters <- ithi.figures:::get_til_clusters(ihc_table, molsubtypes, tils_for_cluster = tils_for_cluster, 
    nclusts = 3)

master_variant_table <- read_variant_file(master_variant_file, db_path)
master_breakpoint_table <- read_variant_file(master_breakpoint_file, db_path)

sig_results <- produce_signature_results(mutsig_dir, master_variant_table, master_breakpoint_table, 
    db_path)
sig_results_snv <- sig_results$snv
sig_results_sv <- sig_results$sv

Analysis

Reviewer 3 has several comments about statistical tests made in the manuscript. These are all valid + easy to address – some are mistakes/typos on my part, others are things we could/should have done, etc.

Fig S4B

Fixed to only show patients with samples in both groups.

cor_tilsubtype_clonal_res <- supp_cor_tilsubtype_clonal(ihc_table, tree_branch_data, 
    tils_for_cluster, nclusts, molsubtypes, db_path)

grid.newpage()
grid.draw(cor_tilsubtype_clonal_res$plots$tilclust_clonalsim)

There, no more ‘orphaned’ points.

Permutation test

Q: What analysis was performed in the sentence “(p > 0.3, permutation test, Figure S4B)” on page 6, which permutation test? Figure S4B is not a permutation test.

This was actually a nested ranks test. It is incorrectly noted as a Wilcoxon signed-rank within the figure legend – this needs to be fixed.

I can see why this would be confused for a MW test, though. What isn’t clearly noted here is that I plot the MEANS for each patient. In fact, all pairwise comparisons in each patient are being compared, not just the mean by cluster equality.

The reason why a Wilcoxon signed-rank is because, to make that type of comparison, patient needs to be accounted for as a random effect. The nested ranks test is one possible extension of a MW/Wilcoxon-type test that uses bootstrapping to compute p-values. See the R package nestedRanksTest.

Tumour purity vs. ITH

Q: In figure S4A, the correlation between “proportion sub-clonal” and “cellularity” is significant (p=0.0169), yet in the manuscript the authors say: “and none of the clonal measures were confounded by tumor purity (all p > 0.2, Figure S4A)”. Page 6.

This is a wording issue. Technically we’re right, but we should make it more obvious by listing out the CLONAL measures, i.e. not including proportion subclonal CN. This is a valid point though – we want to avoid using ambiguous wording wherever possible – not worth saving a few characters to make things more confusing.

Figure 1C

Q: Figure 1C does not say which statistical test was used, and if a Kruskal-Wallis was used, a post hoc analysis should be employed.

Point taken. We’ll do some post-hoc tests to compare pairwise significance between groups.

Since this is a Kruskal-Wallis, the most appropriate post-hoc test to use is the Dunn test.

May also be a good idea to annotate on the plot itself.

nonalpha_cluster_colours <- stringr::str_extract(annotation_colours$til_cluster_colours, 
    "#[0-9A-Z]{6}")
ith_boxplots <- ithi.figures:::plot_ith_boxplots(clonal_measures, stat_types, 
    til_clusters, nonalpha_cluster_colours, force_font = FALSE, scale_factor = 7/18, 
    orientation = "wide", post_hoc = TRUE)

grid.newpage()
grid.draw(ith_boxplots)

TODO: Add a legend for the asterisks. Two asterisks corresponds to 0.05, one to 0.1, and 3 to 0.001.

As per the results above, all measures except clone divergence show a significant difference between all other classes and ES-TIL.

Figure 6B

Q: Statistical analysis of 6B should include post hoc test.

This is going to be removed from the text (because the entire mutation signatures section is going to go poof), but let’s do it anyways for interest (and because I need this for a slide).

label_file <- Sys.glob(file.path(mutsig_dir, "output", "*_labels.tsv"))
labels <- data.table::fread(label_file)
labels <- subset(labels, select = c(library, project, histotype))
labels$library <- stringr::str_replace(labels$library, "^patient_", "")
labels <- labels %>% plyr::rename(c(library = "patient_id"))

signature_labels <- stringr::str_extract(c(colnames(sig_results_snv), colnames(sig_results_sv)), 
    "SN?V-[0-9]+")
signature_labels <- signature_labels[!is.na(signature_labels)]

sigheat_patient_res <- ithi.supp:::get_formal_clusters(sig_results, labels, 
    signature_labels, annotation_colours)

col_clust <- sigheat_patient_res$hc
sample_order <- col_clust$labels[col_clust$order]
clusters <- cutree(col_clust, 3)
general_labels <- as.data.frame(clusters) %>% tibble::rownames_to_column(var = "patient_id") %>% 
    plyr::rename(c(clusters = "cluster"))
general_labels$cluster <- as.character(general_labels$cluster)

message("Plotting boxplots ...")
nonalpha_cluster_colours <- stringr::str_extract(annotation_colours$mutsig_clusters, 
    "#[0-9A-Z]{6}")

# Technically icgc_subtypes hasn't been assigned, but we don't use it
# anymore so w/e
immune_boxplot_fig <- ithi.figures:::plot_immune_boxplot(ith_icgc_bc, level = "patient", 
    nanostring_labels, general_labels, molsubtypes, icgc_specimen_data, annotation_colours, 
    sample_order, db_path, icgc_subtypes, nonalpha_cluster_colours, force_font = FALSE, 
    post_hoc = TRUE)

grid.newpage()
grid.draw(immune_boxplot_fig)

TODO: Modify text to acknowledge that there is significant difference between H-HRD and H-FBI-2 in SOME immune pathways. Or alternatively just be more specific in saying that there isn’t significance in the other categories (instead of saying that they’re similar overall).