README.md

In pontikos/flowBeads: flowBeads: Analysis of flow bead data

flowBeads

This is the public repository for the flowBeads Bioconductor package for working with calibration beads in flow cytometry, based on flowCore.

flowBeads has received some recent attention, in particular I had a few questions about relative normalisation for channels in which the bead manufacturer does not specify the bead MEF.

I think relative normalisation is possible, provided that the number of peaks is consistent across samples. However I believe it is easier to achieve this without using the Bioconductor package which can be quite rigid and cumbersome.

Instead, here is some R code of how to go about it, which also makes use of flowClust for clustering on forward and side scatter in order to filter doublets:

require(flowCore)
# We will use cluster for the pam function (implementation of k medoids).
require(cluster)
# flowClust is a very useful Bioconductor package to do clustering by fitting a mixture of normal distributions.
# It also ignores outliers from the clustering.
require(flowClust)
# This package provides the download.file function.
# If download file doesn't work (returns status code 127) then you can just download the file and save it in the directory
# where you run the script.
require(RCurl)
# This is another repository of mine containing some plotting functions for flow data.
download.file("https://raw.githubusercontent.com/pontikos/FCS/master/fcs.R", destfile = "fcs.R", method = "curl")
source('fcs.R')

This function will retrieve the MFIs of the 8 bead populations:

get.MFI <- function(X) {
  print(length(scatter.channels <- as.character(grep('FSC|SSC',colnames(X),value=TRUE))))
  print(length(fluo.channels <- as.character(grep('^SSC|FSC|Time',colnames(X),invert=T,value=T))))
  # use PCA to reduce the number of dimensions from 6 to 2
  pca <- princomp(X[,scatter.channels])
  pca.X <- pca$scores[,1:2]
  #smoothPlot(pca.X)
  # This might not always be necessary depending.
  # A good way of finding out is to plot the FSC-A against SSC-A.
  res <- flowClust(pca.X,K=4)
  #plot to check result
  #plotClustRes(pca.X,res=res,outliers=FALSE)
  # Pick population which contains the most beads.
  X <- X[which(res@label==which.max(res@w)),]
  X.trans <- apply(X[,fluo.channels],2,logicleTransform())
  res <- pam(X.trans[,fluo.channels],8)
  #check results
  plotClusters(X.trans[,fluo.channels[1:4]],outliers=TRUE, classification=res$clustering,chulls=FALSE) 
  plotClusters(X.trans[,fluo.channels[4:8]],outliers=TRUE,classification=res$clustering,chulls=FALSE)
  # these are now the MFIs
  MFIs <-do.call('rbind',by(X[,fluo.channels],res$clustering,colMeans))
  MFIs <- MFIs[order(MFIs[,1]),]
  return(MFIs)
}

beads1 <- flowCore::read.FCS(file.path("lucas","QC8PeaksBeads_ After Capture Beads_SAS_ARIAIII_CORDOBA_19112014.fcs"))
MFI.beads1 <- get.MFI(beads1@exprs)
beads2 <- flowCore::read.FCS(file.path("lucas","QC8PeaksBeads_32140219_SAS_ARIA_19MAR2015_19MAR2015.fcs"))
MFI.beads2 <- get.MFI(beads2@exprs)

Next it's simply a question of defining the linear transforms which maps the peaks between samples:

fluo.channels <- colnames(MFI.beads1)
trans <-do.call('rbind', lapply( fluo.channels, function(chan) coefficients(lm(log10(MFI.beads1[,chan])  ~ log10(MFI.beads2[,chan]))) ) )
rownames(trans) <- fluo.channels
colnames(trans) <- c('alpha','beta')
# if beta is not an integer, this normalisation is only defined for positive x
# since beta is usually very close to 1 it may be worth just setting to 1
trans[,'beta'] <- 1
normalisation <- lapply(fluo.channels , function(n) return(function(x) 10**trans[n,'alpha'] + x**trans[n,'beta']) )
names(normalisation) <- fluo.channels

Now normalisation contains the transform to compare samples analysed on the same day as beads2 with those analysed at the same day as beads1. So for example if x contains your data from day 2 then you can simply do this to normalise it to day 1:

  chan <- 'APC-AF750-A'
  x <- beads2@exprs[,chan]
  x.norm <- normalisation[[chan]](x)
  plot(density(logicleTransform()(x)))
  lines(density(logicleTransform()(x.norm)),col='red')
  lines(density(logicleTransform()(beads1@exprs[,chan])),col='green')

Ok so these curves don't align well. Why? This is the regression on log10 scale:

plot(log10(MFI.beads1[,chan])  ~ log10(MFI.beads2[,chan]))
abline(lm(log10(MFI.beads1[, chan]) ~ log10(MFI.beads2[, chan])))

This is the regression on linear scale:

plot((MFI.beads1[,chan])  ~ (MFI.beads2[,chan]))
abline(lm((MFI.beads1[, chan]) ~ (MFI.beads2[, chan])))

R my.log10 <- function(x) ifelse(x<=0,0,log10(x)) chan <- 'APC-AF750-A' x <- beads2@exprs[,chan] x.norm <- normalisation[[chan]](x) plot(density(my.log10(x))) lines(density(my.log10(x.norm)),col='red') lines(density(my.log10(beads1@exprs[,chan])),col='green')

(to be continued...)

pontikos/flowBeads documentation built on May 25, 2019, 11:23 a.m.