Re-analysing pancancer dataset with dispersion priors

Published

April 27, 2026

Load dependencies

Load metadata

Show code

# WARNING! THIS INCLUDES UNPIBLISHED DATA
meta_path <- "./data/ash_pancancer/metadata_full.tsv"
meta <- read.csv(meta_path, sep = '\t') 
meta = cbind(run_acc = meta$run_accession, meta)
# From paper:
# RvsP = CR or PR VS. PD or cPD - excluded patients with a BOR of stable disease (SD) (n=29)

# Outcome
meta$RvsP = "R"
meta$RvsP[which(meta$BOR == "PD" | meta$BOR == "cPD")] = "NR" 
meta$RvsP = factor(meta$RvsP, levels = c("NR", "R"))

Load sylph outputs

Show code

sy <- read_sylph("./data/ash_pancancer/combined_q_99.tsv.gz") # q99)

Detected Sylph query output file.

Show code

# annoying renames to match meta V sylph file
colnames(sy) <- gsub("_1", "", colnames(sy))
colData(sy)$Sample_file <- gsub("_1", "", basename(colData(sy)$Sample_file))

# Reorder meta 
meta = meta[match(colnames(sy), meta$run_accession), ]
# get rid of SD
rmidx = which(meta$BOR == "SD")
if(length(rmidx) > 0){
  meta = meta[-rmidx, ]
  sy = sy[, -rmidx]
}

sy <- filter_by_presence(sy, min_nonzero = 8) # filter at 10%

Retained 19139 rows after filtering

Show code

dim(sy)

[1] 19139    77

Show code

# Checks before merging metadata
all(colnames(sy) %in% meta$run_accession)

[1] TRUE

Show code

all(meta$run_accession %in% colnames(sy))

[1] TRUE

Show code

sy = modify_metadata(sy, meta)

Empirical Bayes priors without dispersion term

Show code

design <- as.formula("~ RvsP + histology_cohort.x + age + sex + BMI") 

save_path <- "output_rds/ASH_prior_analysis_zib_q_99_ebp.rds"

if(file.exists(save_path)){
  ZB_fit <- readRDS(save_path)
} else {
  ebp = compute_eb_priors(sy, design, nthreads = parallel::detectCores(), low_cutoff = 0,
                          high_cutoff = Inf)
  # Run with weak prior for prediction
  ZB_fit <- glmZiBFit(sy, design, nthreads = parallel::detectCores(), MAP_prior = ebp)
  saveRDS(ZB_fit, save_path)
}

There are 167 top hits, some of them dodgy

Show code

th = top_hits(ZB_fit)

Found 167 tophits for RvsPR at alpha = 0.05 using holm

Show code

# Looking at the first 10
plot_ani_dist(sy, 'RvsP', th$Contig_name[1:10])

Warning: Removed 689 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Clearly, at least B. animalis and F. nucleatum should NOT be picked…

Post-hoc testing

Show code

th_ph = strainspy::comp_ani_diff_and_posthoc_test(sy, ZB_fit, th, progress = FALSE)

table(th_ph$Comment, useNA = 'always') # now there are only 16 hits


                 <3 nonzero <3 nonzero; too few nonzero 
                         17                          83 
            too few nonzero                        <NA> 
                         26                          41

Show code

plot_ani_dist(sy, 'RvsP', th_ph$Contig_name[is.na(th_ph$Comment)])

Warning: Removed 2554 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Empirical Bayes prior with dispersion terms

Weak

Show code

save_path <- "output_rds/ASH_prior_analysis_zib_q_99_ebp_disp5.rds"
ebp = ZB_fit@priors # priors are the same

# manually add the dispersion parameter
ebp@priors_df = rbind(ebp@priors_df, data.frame(prior = 'normal(0,5)', class = 'fixef_disp', coef = 1))

ebp@priors_df

            prior      class                  coef
1  normal(0,0.27)      fixef                 RvsPR
2  normal(0,0.35)      fixef histology_cohort.xNEN
3  normal(0,0.32)      fixef histology_cohort.xUGB
4  normal(0,0.18)      fixef                   age
5  normal(0,0.34)      fixef                  sexM
6  normal(0,0.12)      fixef                   BMI
7  normal(0,1.55)   fixef_zi                 RvsPR
8  normal(0,3.21)   fixef_zi histology_cohort.xNEN
9  normal(0,2.66)   fixef_zi histology_cohort.xUGB
10 normal(0,0.61)   fixef_zi                   age
11 normal(0,2.12)   fixef_zi                  sexM
12 normal(0,3.21)   fixef_zi                   BMI
13    normal(0,5) fixef_disp                     1

Show code

if(file.exists(save_path)){
  ZB_fit_d5 <- readRDS(save_path)
} else {
  ZB_fit_d5 <- glmZiBFit(sy, design, nthreads = parallel::detectCores(), MAP_prior = ebp)
  saveRDS(ZB_fit_d5, save_path)
}

th_d5 = top_hits(ZB_fit_d5)

Found 150 tophits for RvsPR at alpha = 0.05 using holm

We still have 150 top hits, prior is likely too weak. Sensitivity analysis:

Show code

plot_ani_dist(sy, 'RvsP', th_d5$Contig_name[1:20])

Warning: Removed 1376 rows containing non-finite outside the scale range
(`stat_boxplot()`).

B. animalis is also still there…

Strong

Show code

save_path <- "output_rds/ASH_prior_analysis_zib_q_99_ebp_disp1.rds"
ebp = ZB_fit@priors # priors are the same

# manually add the dispersion parameter
ebp@priors_df = rbind(ebp@priors_df, data.frame(prior = 'normal(0,1)', class = 'fixef_disp', coef = 1))

ebp@priors_df

            prior      class                  coef
1  normal(0,0.27)      fixef                 RvsPR
2  normal(0,0.35)      fixef histology_cohort.xNEN
3  normal(0,0.32)      fixef histology_cohort.xUGB
4  normal(0,0.18)      fixef                   age
5  normal(0,0.34)      fixef                  sexM
6  normal(0,0.12)      fixef                   BMI
7  normal(0,1.55)   fixef_zi                 RvsPR
8  normal(0,3.21)   fixef_zi histology_cohort.xNEN
9  normal(0,2.66)   fixef_zi histology_cohort.xUGB
10 normal(0,0.61)   fixef_zi                   age
11 normal(0,2.12)   fixef_zi                  sexM
12 normal(0,3.21)   fixef_zi                   BMI
13    normal(0,1) fixef_disp                     1

Show code

if(file.exists(save_path)){
  ZB_fit_d1 <- readRDS(save_path)
} else {
  ZB_fit_d1 <- glmZiBFit(sy, design, nthreads = parallel::detectCores(), MAP_prior = ebp)
  saveRDS(ZB_fit_d1, save_path)
}

th_d1 = top_hits(ZB_fit_d1)

Warning in top_hits(ZB_fit_d1): Multiple testing correction using `holm`: No
significant associations detected for coef = 2 at alpha = 0.050000

Now there are no hits. This has killed all signal. Let’s try with an empirical approach.

Ebp prior

Show code

save_path <- "output_rds/ASH_prior_analysis_zib_q_99_ebp_disp_ebp.rds"
ebp =  compute_eb_priors(sy, design, nthreads = parallel::detectCores(), low_cutoff = 0,
                          high_cutoff = Inf, est_disperion_prior = T, progress = FALSE)

plot_prior_bootstrap(ebp, 'dispersion')

Show code

ebp@priors_df

               prior      class                  coef
1     normal(0,0.27)      fixef                 RvsPR
2     normal(0,0.35)      fixef histology_cohort.xNEN
3     normal(0,0.32)      fixef histology_cohort.xUGB
4     normal(0,0.17)      fixef                   age
5     normal(0,0.34)      fixef                  sexM
6     normal(0,0.12)      fixef                   BMI
7     normal(0,1.55)   fixef_zi                 RvsPR
8     normal(0,3.22)   fixef_zi histology_cohort.xNEN
9     normal(0,2.66)   fixef_zi histology_cohort.xUGB
10     normal(0,0.6)   fixef_zi                   age
11    normal(0,2.12)   fixef_zi                  sexM
12    normal(0,3.26)   fixef_zi                   BMI
13 normal(6.99,1.25) fixef_disp                     1

Show code

if(file.exists(save_path)){
  ZB_fit_d_ebp <- readRDS(save_path)
} else {
  ZB_fit_d_ebp <- glmZiBFit(sy, design, nthreads = parallel::detectCores(), MAP_prior = ebp)
  saveRDS(ZB_fit_d_ebp, save_path)
}

th_debp = top_hits(ZB_fit_d_ebp)

Found 97 tophits for RvsPR at alpha = 0.05 using holm

Down to 97 hits…

Show code

plot_ani_dist(sy, 'RvsP', th_debp$Contig_name[1:20])

Warning: Removed 1367 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Not enough to kill B. animalis though…

Post-hoc testing

Show code

th_debp_ph = strainspy::comp_ani_diff_and_posthoc_test(sy, ZB_fit, th_debp, progress = FALSE)

table(th_debp_ph$Comment, useNA = 'always') # now there are only 15 hits, one less from before


                 <3 nonzero <3 nonzero; too few nonzero 
                         10                          48 
            too few nonzero                        <NA> 
                         12                          27

Show code

plot_ani_dist(sy, 'RvsP', th_debp_ph$Contig_name[is.na(th_debp_ph$Comment)])

Warning: Removed 1658 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Common/different hits

Show code

disp_hits = th_debp_ph$Contig_name[is.na(th_debp_ph$Comment)]
ndisp_hits = th_ph$Contig_name[is.na(th_ph$Comment)]

# Common ones
plot_ani_dist(sy, 'RvsP', intersect(disp_hits,ndisp_hits))

Warning: Removed 1542 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Show code

# Only in model without dispersion prior:
plot_ani_dist(sy, 'RvsP', setdiff(ndisp_hits, disp_hits))

Warning: Removed 1012 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Show code

# Only in model WITH dispersion prior:
plot_ani_dist(sy, 'RvsP', setdiff(disp_hits, ndisp_hits))

Warning: Removed 116 rows containing non-finite outside the scale range
(`stat_boxplot()`).