DirichletProcess Model Clustering on TCR Distances

[GWMcElfresh]

Hi folks,

This addresses #55 - specifically, offering a clearer way to parameterize the tcr distance/dianaHeight cutoff one may use for clustering.

There is one high level entrypoint DirichletClusterAnalysis() which runs a nonparameteric bayesian clustering method on the distribution of TCR distances:

One first needs to run the typical processing steps, including choosing a dianaHeight threshold.

seuratObj <- CalculateTcrDistances(
  inputData = seuratObj,
  chains = c("TRA", "TRB"),
  minimumCloneSize = 2,
  calculateChainPairs = TRUE
)

seuratObj <- RunTcrClustering(
  seuratObj_TCR = seuratObj,
  dianaHeight = 20, 
  clusterSizeThreshold = 1
)

Then, you can run this function to estimate clusters and their means/variances/proportion within the cluster

dp <- DirichletClusterAnalysis(
    seuratObj   = seuratObj,
    assayName   = "TRA_fl",
    splitField  = "Population",
    maxSamples  = 1000,
    nIterations = 500,
    verbose = FALSE, 
    nBins = 20, samplesPerBin = 150
)

Optionally, you can also choose to sample along deciles to ensure you sample rare modes. In the above data, a very clear mode is present in the small X values, but that is less likely with unsorted data. You can use the parameters: nBins (10 for deciles, 100 for percentiles) and the number of points per bin samplesPerBin. If nBins * samplesPerBin > maxSamples, then it's downsampled randomly to equal maxSamples from new, more uniform sampled distribution.

dp <- DirichletClusterAnalysis(
    seuratObj   = seuratObj,
    assayName   = "TRA_fl",
    splitField  = "Population",
    maxSamples  = 1000,
    nIterations = 500,
    verbose = FALSE, 
    nBins = 20, samplesPerBin = 150
)

The returned tcrDirichletResult object is a list which enables two helper functions PlotClusterMeans() and PlotMixingProportions() to display clusters, their TCR distances, their proportional abundance, and a metadata variable to compare.

e.g.:

library(patchwork)
(PlotClusterMeans(dp) + Seurat::NoLegend() ) + 
PlotMixingProportions(dp) + 
plot_layout(guides = 'collect') 

This allows me to reparameterize the dianaHeight to a cutoff of ~40 rather than 20, such that I capture the clones that were captured by the tetramers.

I think it's worthwhile to support poisson models as well, but willing to punt for now.