Stable Release

flowPloidy is part of the BioConductor repository. You can install the stable version of the package directly from there. First, install bioconductor in your R session:

source("https://bioconductor.org/biocLite.R")
biocLite()

Then install flowPloidy:

biocLite("flowPloidy")

The examples in this vignette, and the flowPloidy help files, use data in the flowPloidyData package, which is also on BioConductor. (You don’t need flowPloidyData to use flowPloidy on your own data, only for the examples).

biocLite("flowPloidyData")

Development Release

The Bioconductor system emphasizes stable packages. As a consequence, updates are only released twice a year. If you need the latest version of flowPloidy, you may prefer to use the developmental version. This is available directly from the GitHub repository, but requires a few more steps to install.

First, you need to install bioconductor, if you haven’t already:

source("https://bioconductor.org/biocLite.R")
biocLite()

Next, you’ll need to install two dependencies from bioconductor:

biocLite("flowCore")
biocLite("flowPloidyData")

You also need the devtools package, if you don’t have it already:

install.packages("devtools", dependencies = TRUE)

Now you can use devtools to import the latest version of flowPloidy direct from the repository:

library("devtools")
install_github("plantarum/flowPloidy", dependencies = TRUE,
               build_vignettes = TRUE)

This should install flowPloidy along with all its dependencies. However, sometimes R gets confused and will complain about missing packages, even though we’ve asked it to automatically install everything that flowPloidy needs (i.e., dependencies = TRUE). If you see messages in the terminal indicating a package name that isn’t found, try installing that package directly and then re-run the install_github line above. You may need to repeat this process several times to get all the dependencies installed.

Histogram presentation

Calling browseFlowHist() opens a window in your internet browser that will display your histograms. The first one looks like this:

On the right we see the histogram. The graphical elements are:

• shaded area: the raw data
• dotted gray line: the initial model estimate
• red line: the fitted model
• green line: the fitted debris component
• blue line: the fitted components for the G1 and G2 peaks of sample A
• orange line: the fitted components for the G1 and G2 peaks of sample B (the G2 peak is not visible on this plot)
• magenta line: the fitted s-phase component for sample A (too small to see in this plot)
• turquoise line: the fitted s-phase component for sample B (too small to see in this plot)
• purple line: the fitted aggregate component (a very small peak around 240)

You also see coloured circles. These are the initial peak estimates, blue for sample A, and orange for sample B. The solid circles are the G1 estimates, and there is a hollow circle to indicate the G2 peak (only visible for sample A in this plot).

The upper right corner of the plot contains some additional information regarding the model fit:

• RCS: the residual Chi-Square value. This is a rough goodness-of-fit value. Bagwell (1993) suggests that values between 0.7 and 4 indicate an acceptable model fit. See ?RCS for more details in the online manual.
• A and B: the parameter estimates for samples A and B. The three numbers are the peak position, the number of nuclei counted, and the coefficient of variation.
• Linearity: The linearity parameter. This is the ratio between G1 and G2 peaks. By default, we fit it as a model parameter bounded between 1.9 and 2.1¹.

For this example, the default settings have worked perfectly. The model fit is acceptable (RCS 2.208), and the CVs for the peaks are well below the usual 5% threshold for acceptable.

One thing we don’t see is an indication of which peak is the standard, and which is our unknown sample. That is because flowPloidy doesn’t know, and has no way to infer this automatically. You will need to indicate this yourself, based on your understanding of the expected peak sizes and ratios, and perhaps running the standard on its own without an unknown. You can tell flowPloidy which peak is the standard by selecting the value in the “Standard Peak” drop-down selector, on the left side of the plot. By default this is “X”, indicating the standard has not been identified. If you set it to “A” or “B”, the value will be stored in the results. If you don’t set the value, you will still be able to recover the peak ratio from your results, but will have to calculate the pg values yourself. We’ll return to this below.

We’ll ignore the other options in the navigation panel for now (and the gating interface below), but will discuss them further below. For now, click the “Next” button to move on to the next histogram.

Correcting a Failed Model Fit

The second histogram didn’t work at all:

Since the model fitting failed, the plot shows only the raw data and the initial estimates. The fitting parameters are also absent.

Looking at the initial estimates, it’s easy to see what went wrong - a false peak in the noise near position 60 was incorrectly identified as a G1 peak. We can correct this by re-positioning the peak estimates. The “Peak” drop-down list at the left allows you to select which peak to move, and left-clicking on the plot will set the new position for that peak.

Moving Peak B to the second big peak, and Peak A to the first peak, the analysis works as expected:

In most cases, the exact placement of the initial peak isn’t critical. As long as you put it somewhere near the top of the appropriate peak, flowPloidy will be able to find the true value. You can confirm this by clicking around within the peaks: you end up with the same result (within a small rounding error) for quite a range of initial peak positions.

Note that the A peak must have a lower mean (i.e., be further left on the histogram) than the B peak. If you try to put the A peak above the B peak, flowPloidy will swap them around for you, to make the lower peak A and the upper peak B. This is necessary for some of the internal model-fitting code to work properly.

Changing Model Components: Debris Model

Moving on to the next histogram, we see that there were also problems with peak detection here:

We can fix this the same way:

Quick and easy. However, our RCS value is now off the charts! Two things are contributing to this. First, RCS is sensitive to a number of factors, not all of which are related to model fit (see ?RCS). Second, the most appropriate way to model histogram debris varies for different samples. As botanists, we will analyze samples from many different species, tissue types, and using different preparation protocols. In some cases, the Single Cut debris model will be most effective; in others, the Multi-Cut model works better. When your analysis produces high RCS values, as in this case, it may be a clue that the other debris model is more appropriate. We can check this by selecting “MC” from the “Debris Model” drop-down list at the left:

With the Multi-Cut debris model, we have a much more sensible RCS value. You can also see the improved fit visually - note the top of the fitted model is much closer to the raw data with the “MC” debris model, especially in the region between 80 and 140.

For the nitty-gritty details of the debris models, see ?DebrisModels. In practice, use whichever one gives you a lower RCS value. Also keep in mind that there is usually only a very small difference between the two options when it comes to peak means and CVs, even if the difference in RCS can be very large.

Local Minima Traps

Despite its power, non-linear regressions can get stuck in ‘local minima’. Unlike linear regression, there is no single unique solution to a non-linear regression. Fitting a model requires testing out different parameters before determining which combination is best. Sometimes, the algorithm gets stuck on incorrect parameters. Luckily, this is usually easily detected by humans.

An example is file 11, “734.LMD” in our sample data:

Once again, we need to correct our peak estimates. If we move the B peak over we get the following fit:

Other than the RCS value being a little high, the stats for this fit look ok. However, notice that the G2 peak for sample A doesn’t match well with the raw data. That’s strange, as the G1 peak looks fine, as does the sample B peak. This is an example of the model getting trapped in a local minimum.

We can fix this by changing the initial estimates. Notice the initial estimate for the A peak is actually stuck on a ‘noise’ peak just below the real peak. If we shift that point closer to the true peak, we get a more sensible fit:

The new fit is much better, and our RCS value is now down below 4 where we like to see it. Notice, however, that the parameters of most interest (peak position, cell count, and CV) don’t change appreciably; if we hadn’t noticed this problem, our results would have been nearly identical.

The Gating Interface

At the bottom of the browser window you’ll see this:

At the left side is a simple control panel. We can select the X and Y variables using the dropdown menus. All of the data channels in your raw data are available. The zoom slider allows you to zoom in on the lower portion of the plot.

The Y axis dropdown allows you to toggle between two plotting relationships. Y/X plots your chosen X channel on the X axis, and the ratio of the Y and X channels on the Y axis. This is often useful in ‘straightening out’ the plot, making it easier to see the boundaries between cell populations. If you prefer, you can set the Y axis drop down to Y instead, which will give you the usual Y ~ X plot.

Using the Y/X ratio is also useful in this context, because the interface we’re using (built with the R Shiny package) only allows for rectangular selections. So ‘straightening out’ our plot simplifies fitting our target points into a rectangular space.

The final element on the control panel is the Set Gate button. Clicking that will apply the current gate to our data (or remove an existing gate if none is selected).

Interpreting the Gate Plot

On the lower right we have three plot windows. The left-most is the actual gate window, which displays the raw data. The middle window is a preview of the raw data after the gate is applied. The right-most window shows the gate ‘residuals’, the complement of the middle window - it shows a histogram of all the observations that are excluded by the gate.

The plot we’re most interested in is the one on the left. For most of our gating, we set the X-axis to the fluorescence data channel; in this case that’s FL2.A. We plot side scatter (actually, the ratio side scatter/fluorescence) on the Y-axis.

The amount of fluorescence emitted by each particle tells us how much dye they have incorporated; we expect clusters corresponding to the G1 nuclei of our sample and the standard. There may also be clusters for the G2 nuclei. This is exactly what we saw in the previous examples where we didn’t need gating: the unaltered histograms were clear enough to fit our model and extract the parameters we need to determine genome size.

In this case, we have a lot of debris particles that are emitting the same fluorescence intensities as our intact nuclei, so we need another parameter to separate them. That’s where side scatter comes in. Side scatter is a measure of ‘optical complexity’ or ‘granularity’. Simply put, it tells us something about the shape of the particles. And our nuclei will have a characteristic shape. Typically, this shape has lower ‘optical complexity’ than the irregular debris particles.

In order to see this, we’ll need to zoom in on the gate plot. Dragging the zoom slider about halfway across² we see this:

Notice the dense cloud of points on the upper left. These particles have low fluorescence and high granularity. They could be bits of damaged nuclei, or irregularly-shaped cell components that are weakly fluorescent in the FCM laser. This is the debris we need to exclude.

When we do, we’ll uncover our intact nuclei, which are visible along the bottom of the plot in three discrete clusters. The two near the center of the X-axis are the G1 peaks for the sample and standard. The smaller cluster further right is the G2 peak for the standard. Zooming in further we can see these more clearly:

Here, notice that the nuclei form tight, dark clusters, with particles trailing away upwards. Those upward ‘side-scatter tails’ are probably intact nuclei with bits of debris or fluorescent compounds stuck to them. This makes their shapes more irregular, increasing their side scatter values. Unfortunately, there is no clear dividing line between ‘clean nuclei’ and ‘goopy nuclei’.

In practice, we place our gates high enough to retain all the particles in the densest section of the clusters, and extend them upwards to include the bottom portion of the tails. There is an element of subjectivity here: we are balancing biasing our results with goopy nuclei by including too much of the tail, and reducing the amount of precious data in the analysis by being too aggressive with our gate. Including more of the goopy nuclei increases our CV, and excluding them lowers our cell count.

Applying a Gate

Use the mouse to left-click and drag a box around the cell clouds:

The upper boundary of the gate is important - it will have a direct influence on both our cell counts and CVs. However, the left and right boundaries have very little impact. The model-fitting process does a very nice job of determining where our peaks start and stop, so we don’t need to worry about drawing ‘correct’ boundaries on the left and right. Consequently, we set our gates to include the full width of the plot.

We can see the effect of the gate in the middle window, and the events that we’re excluding in the right window. At this point, we can see three peaks in the histogram preview. To apply this gate to the analysis, click the Set Gate button. The result should look something like this:

The histogram plot now includes the –GATED– flag. This is a caution to remind us that applying a gate may influence the model-fitting process in a number of ways, including how effective the debris model components (SC and MC) capture the remaining debris, and the absolute value of RCS. In general, these issues are outweighed by the improved fit for high-debris samples.

The final outcome here is we have a usable result. Sadly for those of us working with Vaccinium, we’re usually limited to very low cell counts. But given the challenges presented by this difficult tissue, we have succeeded in isolating the nuclei necessary to get a GC estimate.