R/afRemover.R

In ellispatrick/AFremover: A method to identify and remove autofluorescent signals from multi-channel images post acquisition

#' Calculate textural features for each region of interest (ROI). Measures size, pixel correlation, standard deviation, and kurtosis of ROIs across two channels, output as a dataframe.
#' @param im1 A matrix containing the intensities of image 1.
#' @param im2 A matrix containing the intensities of image 2.
#' @param mask A matrix of the mask with labelled ROIs.
#' @param imBit bit depth of im1 and im2.
#' @import moments
#' @export
#' @return A dataframe of relevant measurements for clustering.

afMeasure <- function(im1, im2, mask, imBit = 16) {

  imDimension <- dim(im1) # Currently assumes all images the same size

  im1PixelVals <- split(im1, mask)
  im2PixelVals <- split(im2, mask)

  im1PixelVals$"0" <- NULL
  im2PixelVals$"0" <- NULL

  No <- names(im1PixelVals)
  Size <- NULL
  Corr <- NULL
  SD1 <- NULL
  SD2 <- NULL
  Kurt1 <- NULL
  Kurt2 <- NULL

  for (i in names(im1PixelVals)) {
    im1Pixels <- im1PixelVals[[i]]
    im2Pixels <- im2PixelVals[[i]]

    Size[i] <- length(im1Pixels)
    Corr[i] <- cor(im1Pixels, im2Pixels)
    SD1[i] <- sd(im1Pixels)
    SD2[i] <- sd(im2Pixels)
    Kurt1[i] <- moments::kurtosis(im1Pixels)
    Kurt2[i] <- moments::kurtosis(im2Pixels)
  }

  return(data.frame(No, Size, Corr, SD1, SD2, Kurt1, Kurt2))

}


#' Classifies ROIs as autofuorescent and non-autofluorescent based on the measurements obtained by afMeasure.
#' @param mask A matrix of 0's and 1's describing the mask of the ROIs
#' @param df A data frame generated by afMeasure.
#' @param minSize The minimum ROI size (pixel area) to classify - should be greater than 1.
#' @param maxSize The maximum ROI size (pixel area) to classify.
#' @param corr A correlation cut-off for autofluorescence; all objects above this cut-off within the identified autofluorescence cluster is classified as auofluorescence.
#' @param kAuto Logical to estimate the number of clusters between 3 and a user provided number of clusters defined by 'k'; if TRUE, will perform this estimate.
#' @param k If 'kAuto' = FALSE: number of clusters used when performing k-means; if 'kAuto' = TRUE : maximum number of clusters for estimating an optimal k.
#' @import moments
#' @import EBImage
#' @export
#' @return A matrix containing a mask of all autofluorescent objects.
#' @examples
#' set.seed(51773)
#' ## Read in images.
#' imageFile1 = system.file("extdata","ImageB.CD3.tif", package = "AFremover")
#' imageFile2 = system.file("extdata","ImageB.CD11c.tif", package = "AFremover")
#' im1 <- EBImage::readImage(imageFile1)
#' im2 <- EBImage::readImage(imageFile2)
#'
#' ## Transform the image.
#' im1 = im1/max(im1)
#' im2 = im2/max(im2)
#'
#' combined <- EBImage::rgbImage(green=sqrt(im1), red=sqrt(im2))
#' EBImage::display(combined, all = TRUE, method = 'raster')
#'
#'
#' ## Create masks using EBImage.
#'
#' # Find tissue area
#' tissue1 = im1 > 2*min(im1)
#' tissue2 = im2 > 2*min(im2)
#'
#' # Calculate thresholds
#' imThreshold1 <- mean(im1[tissue1]) + 2*sd(im1[tissue1])
#' imThreshold2 <- mean(im2[tissue2]) + 2*sd(im2[tissue2])
#'
#' # Calculate masks.
#' mask1 <- EBImage::bwlabel(im1 > imThreshold1)
#' mask2 <- EBImage::bwlabel(im2 > imThreshold2)
#'
#' # Calculate intersection mask
#' mask <- intMask(mask1,mask2)
#'
#'
#' ## Calculate textural features.
#' df <- afMeasure(im1, im2, mask)
#'
#' ## Alternatively
#' ## Correlation only
#' # afMask <- afIdentify(mask, df, minSize = 100, maxSize = Inf, corr = 0.6)
#'
#' ## Clustering with given k
#' # afMask <- afIdentify(mask, df, minSize = 100, maxSize = Inf, k = 6)
#'
#' ## Clustering with estimated k.
#' afMask <- afIdentify(mask, df, minSize = 100, maxSize = Inf, k = 20, kAuto = TRUE)
#'
#'
#' ## Remove autofluorescence from images
#' im1AFRemoved <- im1
#' im2AFRemoved <- im2
#' im1AFRemoved[afMask != 0] <- 0
#' im2AFRemoved[afMask != 0] <- 0
#'
#' combinedRemoved <- EBImage::rgbImage(green = sqrt(im1AFRemoved), red = sqrt(im2AFRemoved))
#' img_comb = EBImage::combine(combined, combinedRemoved)
#' EBImage::display(img_comb, all = TRUE, method = 'raster',nx = 2)
#'
#' ##Or
#' ##Exclude AF ROIs
#'
#' exclude1 = unique(mask1[afMask>0])
#' mask1Removed = mask1
#' mask1Removed[mask1Removed%in%exclude1] = 0
#' exclude2 = unique(mask2[afMask>0])
#' mask2Removed = mask2
#' mask2Removed[mask2Removed%in%exclude2] = 0

afIdentify <- function(mask, df, minSize = 100, maxSize = Inf, corr = -1, kAuto = FALSE, k = 1) {

  dfToCluster <- df[which(df$Size >= minSize & df$Size <= maxSize),]

  corrAF <- dfToCluster$Corr > corr

  if (k > 1) {
    dfToCluster$Corr <- atanh(dfToCluster$Corr)
    dfToCluster$SD1 <- log(dfToCluster$SD1)
    dfToCluster$SD2 <- log(dfToCluster$SD2)
    dfToCluster$Kurt1 <- log(dfToCluster$Kurt1)
    dfToCluster$Kurt2 <- log(dfToCluster$Kurt2)

    corrVals <- dfToCluster$Corr

    cols <- c('Corr', 'SD1', 'SD2', 'Kurt1', 'Kurt2')

    dfToCluster[, cols] <-
      t(t(dfToCluster[, cols]) / apply(dfToCluster[, cols], 2, sd, na.rm = TRUE))

    if (kAuto) {
      kVals <- 2:k
      testVals <- numeric(length(kVals))

      clusterResults <- data.frame(matrix(length(kVals), nrow = nrow(dfToCluster)))
      afClustNum <- numeric(length(kVals))

      for (i in 1:length(kVals)) {
        dataCluster <-
          kmeans(
            dfToCluster[, cols],
            kVals[i],
            nstart = 20,
            iter.max = 10000
          )
        clusterResults[[i]] <- dataCluster$cluster

        clusterCorrVals <- split(corrVals, clusterResults[[i]])

        cM <- unlist(lapply(clusterCorrVals, mean))
        cL <- unlist(lapply(clusterCorrVals, length))
        cS <- names(sort(-cM[cL>2]))

        afID1 <- cS[1]
        afID2 <- cS[2]
        testVals[i] <- t.test(clusterCorrVals[[afID1]], clusterCorrVals[[afID2]])$stat
        afClustNum[i] <- afID1
      }

      x1 <- 2
      x2 <- k
      y1 <- testVals[1]
      y2 <- testVals[length(testVals)]

      d1 <- (y2 - y1) *  kVals
      d2 <- (x2 - x1) * testVals
      d3 <- x2 * y1 - y2 * x1

      d <- d1 - d2 + d3

     kBest <- which.max(d)-1
     clustAF <- clusterResults[[kBest]] == afClustNum[kBest]


    } else {
      dataCluster <-
        kmeans(
          dfToCluster[, cols],
          k,
          nstart = 20,
          iter.max = 10000,
          algorithm = 'Lloyd'
        )
      clusterResults <- dataCluster$cluster

      clusterCorrVals <- split(corrVals, clusterResults)

      afID1 <- order(-unlist(lapply(clusterCorrVals, mean)))[1]

      clustAF <- clusterResults == afID1
    }

  } else {
    clustAF <- TRUE
  }

  afIdx <- dfToCluster$No[corrAF & clustAF]

  afMask <- mask
  afMask[!(afMask %in% afIdx)] <- 0

  return(afMask)
}



#' Calculate a union of the intersection masks for two masks.
#' @param mask1 A matrix containing a mask with annotated ROIs.
#' @param mask2 A matrix containing a mask with annotated ROIs.
#' @param minSize A numeric specifying the minimum number of pixels for a ROI to be retained.
#' @export
#' @return A matrix of a new mask.

intMask <- function(mask1, mask2, minSize = 5) {
mask <- mask1
mask[] <-  as.numeric(as.factor(paste(mask1,mask2,'_')))-1
mask[!mask%in%names(which(table(mask)>=minSize))] = 0
mask
}