Beta-Binomial model for signals

library(signals)

Background

This vignette demonstrates how to use a beta-binomial model instead of a binomial model for the HMM emission model.

Inference

data("haplotypes")
data("CNbins")
haplotypes <- format_haplotypes_dlp(haplotypes, CNbins)
hscn <- callHaplotypeSpecificCN(CNbins, haplotypes, likelihood = "binomial")

print(hscn)
#> Haplotype specific copy number object 
#>  
#> Number of cells: 250
#> Bin size: 0.5 Mb 
#> Inferred LOH error rate: 0.018 
#> Emission model for HMM: binomial 
#> Average distance from median to expected BAF = 0.0027 
#> Average ploidy = 4 
#> Average number of segments = 79

hscn_bb <- callHaplotypeSpecificCN(CNbins, haplotypes, likelihood = "betabinomial")
#> Iteration 1: loglikelihood = -313171.9553
#> Iteration 2: loglikelihood = -309421.3946
#> Iteration 3: loglikelihood = -308910.1132
#> Iteration 4: loglikelihood = -308886.1312
#> Iteration 5: loglikelihood = -308885.7462
#> Iteration 6: loglikelihood = -308885.7422
#> Iteration 7: loglikelihood = -308885.7421

print(hscn_bb)
#> Haplotype specific copy number object 
#>  
#> Number of cells: 250
#> Bin size: 0.5 Mb 
#> Inferred LOH error rate: 0.018 
#> Emission model for HMM: betabinomial 
#>   Inferred over dispersion: 0.01
#>   Tarones Z score: 185.119
#> Average distance from median to expected BAF = 0.0027 
#> Average ploidy = 4 
#> Average number of segments = 79

We see that the model inferred a modest degree of overdispersion in the data, but that this was siginificantly different from 0 dispersion (Tarones Z-score is high).

Below I’ll plot the heatmaps for the two output, which in this case are very similar.

Binomial model

plotHeatmap(hscn, plotcol = "state_BAF", plottree = FALSE, spacer_cols = 15)

Beta-Binomial model

plotHeatmap(hscn_bb, plotcol = "state_BAF", plottree = FALSE, spacer_cols = 15)

QC plot

To interrogate the difference between the binomial and beta-binomial models we can plot the BAF and overlay the two models.

plotBBfit(hscn_bb)

We can check how similar the results are by comparing the two dataframes. We find for this data that the results are very similar.

print(dim(hscn$data))
#> [1] 1090750      20
print(dim(hscn_bb$data))
#> [1] 1090750      20
all.equal(orderdf(hscn$data), orderdf(hscn_bb$data))
#>  [1] "Component \"alleleA\": 'is.NA' value mismatch: 14554 in current 14537 in target"    
#>  [2] "Component \"alleleB\": 'is.NA' value mismatch: 14554 in current 14537 in target"    
#>  [3] "Component \"totalcounts\": 'is.NA' value mismatch: 14554 in current 14537 in target"
#>  [4] "Component \"BAF\": 'is.NA' value mismatch: 14554 in current 14537 in target"        
#>  [5] "Component \"state_min\": Mean relative difference: 0.6283041"                       
#>  [6] "Component \"A\": Mean relative difference: 0.3738565"                               
#>  [7] "Component \"B\": Mean relative difference: 0.5462495"                               
#>  [8] "Component \"state_AS_phased\": 7902 string mismatches"                              
#>  [9] "Component \"state_AS\": 2423 string mismatches"                                     
#> [10] "Component \"LOH\": 1978 string mismatches"                                          
#> [11] "Component \"phase\": 5350 string mismatches"                                        
#> [12] "Component \"state_phase\": 6335 string mismatches"                                  
#> [13] "Component \"state_BAF\": Mean relative difference: 0.7569649"

Another thing we can check is the total number of segments identified by the HMM. Using the beta-binomial model should reduce the influence of any noisy regions and thus we might expect fewer segments. We do observe fewer segments in the beta-binomial model.

segs <- create_segments(hscn$data, field = "state_AS_phased")

segs_bb <- create_segments(hscn_bb$data, field = "state_AS_phased")

print(dim(segs))
#> [1] 36082     5
print(dim(segs_bb))
#> [1] 32642     5

Summary

With this dataset, the difference between the nbinomial model and beta-binomial model is minimal but produces a statistically superior fit, in lower quality datasets with for example lower depth of coverage the beta-binomial model may may provide more noticeable differences.