R/Flyception.R

In tkatsuki/FlyceptionR: R scripts and utilities for Flyception data analysis

#' FlyceptionR main script
#'
#' @param dir path to the directory that contains the data
#' @param prefix prefix used for output files
#' @param autopos logical. Perform camera alignment using FNCC?
#' @param interaction logical. Perform interaction detection? Requires two flies.
#' @param reuse logical. Reuse .RDS files?
#' @param fmf2tif logical. Convert fly-view and arena-view fmf files into tif format?
#' @param zoom numeric. Zoom factor between fly-view and fluo-view cameras.
#' @param FOI a vector of two numbers indicating the first and last frames to be analyzed. If not specified, all frames will be analyzed.
#' @param binning integer. Binning of the fluo-view camera.
#' @param fluo_flash_thresh numeric. A threshold for detecting flashes in a fluo-view video.
#' @param fv_flash_thresh integer. A threshold for detecting flashes in a fly-view video.
#' @param av_flash_thresh integer. A threshold for detecting flashes in a arena-view video.
#' @param dist_thresh numeric. A distance threshold for detecting fly-fly interactions.
#' @param rotate_camera integer. Angle of the fluo-view camera.
#' @export
#' @examples
#' FlyceptionR()
#'

FlyceptionR <- function(dir, prefix, autopos=T, interaction=T, stimulation=F, reuse=T,
                        fmf2tif=T, zoom=0.85, FOI=F, binning=1, fluo_flash_thresh=10000,
                        fv_flash_thresh=135, av_flash_thresh=100, dist_thresh=4,
                        rotate_camera=-180){

  # TO DO
  # Output: frid, frida, registered images, good frames
  # Create frame
  # Parallelism

  # ensure that dir ends with a slash:
  dir = paste0(normalizePath(dir, winslash="/"), "/")

  ## Part 0. Initialization
  # Start logging
  rlogging::SetLogFile(base.file=paste0(prefix, "_log.txt"), folder=dir)
  message(dir)

  # Prepare filenames
  output_prefix <- paste0(dir, prefix)
  fluo_view_tif <- paste0(dir, list.files(dir, pattern="ome\\.tif$"))
  fly_view_fmf <- paste0(dir, list.files(dir, pattern="^fv.*fmf$"))
  arena_view_fmf <- paste0(dir, list.files(dir, pattern="^av.*fmf$"))

  # Create ouput directory if it doesn't exist
  dir.create(output_prefix, showWarnings=FALSE, recursive=TRUE)
  # Save files within the output directory
  output_prefix <- file.path(output_prefix, prefix)

  ## Part 1. Detect flash
  message("Detecting flash in fluo-view")
  fluo_flash <- detect_flash(input=fluo_view_tif,
                             type="fluo",
                             output=output_prefix,
                             flash_thresh=fluo_flash_thresh,
                             reuse=reuse)
  message("Detecting flash in fly-view")
  fly_flash <- detect_flash(input=fly_view_fmf,
                            type="fly",
                            output=output_prefix,
                            flash_thresh=fv_flash_thresh,
                            reuse=reuse)
  message("Detecting flash in arena-view")
  arena_flash <- detect_flash(input=arena_view_fmf,
                              type="arena",
                              output=output_prefix,
                              flash_thresh=av_flash_thresh,
                              reuse=reuse)

  ## Part 2. Syncing frames and generate frame IDs
  syncing <- sync_frames(dir=dir,
                         fluo_flash=fluo_flash,
                         fly_flash=fly_flash,
                         arena_flash=arena_flash,
                         output=output_prefix,
                         reuse=reuse)

  ## Part 3. Analyze trajectories
  trj_res <- analyze_trajectories(dir=dir,
                                  output=output_prefix,
                                  fpsfv=syncing$fpsfv,
                                  interaction=interaction)

  ## Part 4. Detect stimulus
  if(stimulation==T){
    message("Detecting stimulus")
    fvtrj <- read.table(paste0(dir, list.files(dir, pattern="fv-traj-")))
    stimulus <- which(fvtrj[,10]==1)
    if(length(stimulus)==0){
      fridstim <- NA
      message(paste0("No stimulus was detected."))
    } else {
      stimfr <- sapply(stimulus, function(x) which.min(abs(syncing$frid-x)))
      message(paste0("Stimuli were given at the following frames:"))
      message(stimfr)
      dfstim <- data.frame(flview=stimfr, flyview=syncing$frid[stimfr], arenaview=syncing$frida[stimfr])
      write.table(dfstim, paste0(dir, prefix, "_fridstim.txt"))
    }
  }

  ## Part 5. Detect interaction
  if(interaction==T){
    message("Detecting interaction")
    closefr <- which(trj_res$flydist < dist_thresh)
    closefrid <- sapply(closefr, function(x) which.min(abs(syncing$frida-x)))
    write.table(closefrid, paste0(dir, prefix, "_closefrid.txt"))
  }

  ## Part 6. Load images
  message(sprintf("Reading %s", fluo_view_tif))

  # Analyze only part of the movie?
  if(FOI!=F && length(FOI)==2){
    flimg <- dipr::readTIFF2(fluo_view_tif, start=FOI[1], end=FOI[2])
    message(sprintf("Fluo-view frames from %d to %d will be analyzed.", FOI[1], FOI[2]))
    frid <- syncing$frid[FOI[1]:FOI[2]]
    frida <- syncing$frida[FOI[1]:FOI[2]]
  }else{
    message("All frames will be analyzed.")
    flimg <- dipr::readTIFF2(fluo_view_tif)
    frid <- syncing$frid
    frida <- syncing$frida
  }
  # Crop fluo-view movie for speed
  if(dim(flimg)[1] > 130){
    flimg <- flimg[round((dim(flimg)[1] - 128)/2):(round((dim(flimg)[1]/2+128/2))-1),
                   round((dim(flimg)[2] - 128)/2):(round((dim(flimg)[2]/2+128/2))-1),]
  }
  flimgrt <- EBImage::rotate(EBImage::flip(flimg), rotate_camera)
  # Load fly-view camera images
  fvimgl <- dipr::readFMF(fly_view_fmf, frames=frid)
  # Load arena-view camera images
  avimgl <- dipr::readFMF(arena_view_fmf, frames=frida)
  EBImage::writeImage(avimgl/255, file=paste0(output_prefix, "_avimgl_fr_", frida[1], "-", tail(frida, n=1), ".tif"))
  rm(avimgl)

  ## Part 7. Detect window on the head
  fvimgbwbrfh <- detect_window(fvimgl=fvimgl, output=output_prefix, reuse=reuse)

  ## Part 8. Position calibration
  flref <- dipr::readTIFF2(fluo_view_tif, frames=fluo_flash$flflashes[1])
  flref <- EBImage::normalize(EBImage::rotate(EBImage::flip(flref), rotate_camera))
  fvref <- dipr::readFMF(filepath=fly_view_fmf,
                         frames=fly_flash$fvflashes[1])[,,1]/255

  center <- align_cameras(source=fvref,
                          template=flref,
                          output=output_prefix,
                          center=c(0, 0),
                          zoom=zoom,
                          autopos=autopos)

  ## Part 9. Image registration
  registered_images <- register_images(frid=frid,
                                       fvimgl=fvimgl,
                                       flimgrt=flimgrt,
                                       fvimgbwbrfh=fvimgbwbrfh,
                                       angles=trj_res$angles,
                                       zoom=zoom,
                                       center=center,
                                       output=output_prefix,
                                       reuse=reuse)

  ## Part 10. Look for good frames based on size, position, motion, and focus
  goodfr <- find_goodframes(window_mask=fvimgbwbrfh,
                            fvimgl=fvimgl,
                            output=output_prefix,
                            reuse=reuse)

  ## Part 11. Calculate fluorescence intensity in the brain window
  message("Measuring fluorescence intensity...")
  if(file.exists(paste0(output_prefix, "_intensity_br.RDS"))==T & reuse==T){
    message("Using RDS file")
    intensity_br <- readRDS(paste0(output_prefix, "_intensity_br.RDS"))
  }else{
    intensity_br <- colSums(registered_images$fvimgbwbrfhregimg*registered_images$flimgreg, dims=2)/as.integer(goodfr$objsize)
    saveRDS(intensity_br, paste0(output_prefix, "_intensity_br.RDS"))
  }
  intensity <- zoo::na.approx(intensity_br)

  ## Part 12. Plot delta F over F0
  message("Creating dF/F0 plots")
  F0int <- mean(intensity[1:5])
  deltaFint <- intensity - F0int
  dFF0int <- deltaFint/F0int * 100
  dat <- data.frame(x=(1:length(dFF0int)), y=dFF0int, d=trj_res$flydist[frida])
  p <- ggplot2::ggplot(data=dat, ggplot2::aes(x=x, y=y)) +
    ggplot2::geom_smooth(method="loess", span = 0.4, level=0.95) +
    ggplot2::ylim(-5, 10)
  # add distance plot if appropriate
  if (!rlang::is_na(trj_res$flydist)) {
    P = P + ggplot2::geom_line(data=dat, ggplot2::aes(x=x, y=d))
  }
  ggplot2::ggsave(filename = paste0(output_prefix, "_dFF0int.pdf"), width = 8, height = 8)

  ## Part 13. Create delta F over F0 pseudocolor representation
  message("Calculating dF/F0 images...")
  dF_F0_image(flimgreg=registered_images$flimgreg,
              fvimgbwbrfhregimg=registered_images$fvimgbwbrfhregimg,
              regimgi=registered_images$regimgi,
              colmax=100, cmin=30, cmax=200,
              goodfr=goodfr$goodfr,
              output=output_prefix)

  ## Part 14. ROI measurement
  # Create ROI mask
  # Rectangle ROI example
  ROI_mask <- array(0, dim=dim(registered_images$flimgreg)[1:2])
  ROI_mask[120:(120+10-1),120:(120+10-1)] <- 1
  # Circular ROI example
  # ROI_mask <- array(0, dim=dim(registered_images$flimgreg)[1:2])
  # EBImage::drawCircle(img=ROI_mask, x=120, y=120, radius=5, col=1, fill=T)

  ROI_dFF0 <- measureROI(img=registered_images$flimgreg,
                         mask=ROI_mask,
                         output=output_prefix,
                         goodfr=goodfr$goodfr)

  ## Part 15. Create trajectory of the flies
  message("Creating trajectory of the flies...")
  pdf(file= paste0(output_prefix, "_trackResult.pdf"), width = 4.4, height = 4, bg = "white")
  par(plt = c(0, 1, 0, 1), xaxs = "i", yaxs = "i")
  plot(trj_res$trja[frida,1]*10, -trj_res$trja[frida,2]*10,
       type = "l", lty = 1, col="red",
       axes = F, xlim = c(-240, 240), ylim = c(-220, 220))
  par(new=T)
  plotrix::draw.ellipse(0,0,11.0795*20,10*20)
  dev.off()

  ## Part 16. Convert fmf to tif format, if tif doesn't exist yet
  if(fmf2tif==T){
    for (fvFile in list.files(dir, pattern="^fv.*fmf$")) {
      tiffFile = paste0(substr(fvFile, 1, nchar(fvFile) - 4), ".tif")
      if (!file.exists(paste0(dir, tiffFile))) {
        message(paste0("fv.tiff doesn't exist yet: ", tiffFile))
        dipr::fmf2tif(paste0(dir, fvFile), skip=10)
      }
    }
    for (avFile in list.files(dir, pattern="^av.*fmf$")) {
      tiffFile = paste0(substr(avFile, 1, nchar(avFile) - 4), ".tif")
      if (!file.exists(paste0(dir, tiffFile))) {
        message(paste0("av.tiff doesn't exist yet: ", tiffFile))
        dipr::fmf2tif(paste0(dir, avFile), skip=2)
      }
    }
  }

  message("Finished processing!")
  gc()
}