popPCR: EM Mixture Model fitting of dPCR droplet fluorescence

In popPCR: Classify Digital PCR Droplets by Fitting Fluorescence Populations

Description

Estimates target concentration by counting positive droplets with Poisson correction. Positive, negative, and rain populations are fitted using EM. Droplets are then classified using Maximum A Posteriori rule

Usage

popPCR(
  x,
  dist,
  volSamp = 20,
  volDrp = 0.85,
  maxComponents = Inf,
  negProbThres = 1e-07,
  useOnlyNegProbThres = FALSE
)

Arguments

x	numeric, vector of droplet fluorescence amplitude
dist	character, distribution of the mixture models ("normal", "skewed-normal", "t", "skewed-t")
volSamp	numeric, sample volume in microliter
volDrp	numeric, droplet (or partition) volume in nanoliter
maxComponents	numeric, maximum number of components (e.g. populations)
negProbThres	numeric, if only one population was detected, then its assumed as a negative population. Droplets will be classified as positive if its probability given the population < negProbThres.
useOnlyNegProbThres	logical, if TRUE, then droplets will be classified as positive if its probability given the leftmost population < negProbThres. Default is FALSE, i.e. classification is done by Maximum A Posteriori rule.

Value

Returns a result.popPCR S4 class object with attributes

• classification - character, vector of droplet classification

• dist - character, user-specified parameter for the mixture model

• dropletCount - list, droplet classification count

• em - list, returned value of EMMIXskew's EmSkew()

• estConc - list, estimated target concentration as lambda and sample concentration (with 95% CI)

• G - numeric, number of components fitted

• memberProb - list, component membership probability of all droplets

Examples

library(popPCR)

# Plot histograms of available data
hist(x_onePop, breaks = 100)
hist(x_twoPop, breaks = 100)
hist(x_multiPop, breaks = 100)

# ---- Mixture model fitting ---- #
# Example 1. One population sample
result <- popPCR(x_onePop, dist = "t")
printSummaryConc(result)
#    Output:
#        Populations detected : 1
#        Total droplets : 8000
#        Positive : 1 (0.01%)
#        Negative : 7999 (99.99%)
#
#        Target copies in sample          : 2.9414 ( 95% CI: [ -2.8237 , 8.7064 ] )
#        Mean target copies per partition : 1e-04 ( 95% CI: [ -1e-04 , 4e-04 ] )
printSummaryFit(result)
#    Output:
#        Results of fitting a 1-component t mixture model
#
#        Negative Population
#        Mix prop. : 1
#        Mu        : 1024.1614
#        Sigma     : 35253.1747
#        Dof       : 2.005

# (Option) increase negProbThres to classify negative droplets more strictly
result <- popPCR(x_onePop, dist = "t", negProbThres = 1e-4)
printSummaryConc(result)
#    Output:
#        Populations detected : 1
#        Total droplets : 8000
#        Positive : 691 (8.64%)
#        Negative : 7309 (91.36%)
#
#        Target copies in sample          : 2125.5312 ( 95% CI: [ 1966.9936 , 2284.0688 ] )
#        Mean target copies per partition : 0.0903 ( 95% CI: [ 0.0836 , 0.0971 ] )

# Example 2. Two population sample
result <- popPCR(x_twoPop, dist = "t")
printSummaryConc(result)
#    Output:
#        Populations detected : 2
#        Total droplets : 10254
#        Positive : 8693 (84.78%)
#        Negative : 1561 (15.22%)
#
#        Target copies in sample          : 44290.3819 ( 95% CI: [ 43215.6408 , 45365.1231 ] )
#        Mean target copies per partition : 1.8823 ( 95% CI: [ 1.8367 , 1.928 ] )
printSummaryFit(result)
#    Output:
#        Results of fitting a 2-component t mixture model
#
#        Negative Population
#        Mix prop. : 0.1522
#        Mu        : 2136.7435
#        Sigma     : 4126.8357
#        Dof       : 12.3562
#
#        Positive Population
#        Mix prop. : 0.8478
#        Mu        : 7580.1275
#        Sigma     : 42621.1894
#        Dof       : 2.415


# Example 3. Multiple population sample
result <- popPCR(x_multiPop, dist = "t", maxComponents = 4)
printSummaryConc(result)
#    Output:
#        Populations detected : 4
#        Total droplets : 1814
#        Positive : 44 (2.43%)
#        Negative : 1252 (69.02%)
#        Rain (1) : 258 (14.22%)
#        Rain (2) : 260 (14.33%)
#
#        Target copies in sample          : 8724.5195 ( 95% CI: [ 7999.0578 , 9449.9812 ] )
#        Mean target copies per partition : 0.3708 ( 95% CI: [ 0.34 , 0.4016 ] )

# In the output above, we see 2 rain populations! Let's examine its plot.
plot(stats::density(x_multiPop))
# We can see that Rain (1) is very close to the Negative population.
# Let's include droplets in Rain (1) in the negative droplet count.
nNegative <- result@dropletCount$neg + result@dropletCount$rain1
nTotal <- result@dropletCount$total
# Re-estimate concentration as follows
newEstimates <- calculateConc(nNegative, nTotal, volSamp = 20, volDrp = 0.85)
newEstimates
#    Output:
#       $lambda
#          lambda     lower     upper
#       0.1834247 0.1627763 0.2040731
#
#       $conc
#           conc    lower    upper
#       4315.875 3830.031 4801.719

popPCR documentation built on March 10, 2021, 1:08 a.m.