FlowType.R

 #!/usr/bin/env Rscript

require(docopt)
require(methods)

"
Usage:
FlowType.R (-h | --help | --version)
FlowType.R DIR

Description:   This script applies the flowType algorthim to high parameter cytometry data
Options:
--version       Show the current version.

Arguments:

DIR    Provide directory for cytools.args.Rdata to be found, automatically generated by invoking cyttools.R
" -> doc


args <- docopt(doc)

ARGS_DIR <- args$DIR

cat("\nLoading arguments from", ARGS_DIR, "\n")

load(paste(ARGS_DIR, "cyttools.args.Rdata", sep = ""))

RESULTS_DIR <- args$OUT

source("cyttoolsFunctions.R")

dir <- args$DIR # grabs directory from initial cyttools call
file <- list.files(dir ,pattern='.fcs$', full=TRUE) # captures all FCS files in the directory

targets <- read.delim(args$PANEL)
colsToCheck <- c("Ignore", "TransformCofactor", "Lineage", "Functional", "NRS", "PartitionsPerMarker")
if(checkDesignCols(targets, colsToCheck)){
  missingCols <- colsToCheck[which(colsToCheck %in% colnames(targets) == F)]
  cat("\n\nERROR: PANEL file does not include required columns.
      \n\nMissing Columns:", missingCols,
      "\n\nPlease run cyttools.R --makePanelBlank and cyttools.R --computeNRS to generate compatible panel file.\n\nStopping cyttools.R\n\n")
  q()
}


lineage_markers <- targets$name[targets$Lineage == 1]
functional_markers <- targets$name[targets$Functional == 1]

if(args$transform == "logicle"){
  # read in fcs files
  ncfs <- read.ncdfFlowSet(file)
  
  chnls <- colnames(ncfs)[grep("SSC|FSC|Time|\\-H", colnames(ncfs), invert = T)]
  safe_estimate_logicle <- safely(estimateLogicle)
  transFuncts <- fsApply(ncfs, safe_estimate_logicle, channels = paste0("^", chnls)) %>%
    modify_depth(1, 1) %>%
    discard(is_null)
  
  safe_transform <- safely(transform)
  
  for ( i in 1:length(transFuncts)){
    ncfs_trans <- safe_transform(ncfs, transFuncts[[i]])
    if(is.null(ncfs_trans$error)){
      flowSet.trans <- as.flowSet(ncfs_trans$result)
      break
    }else if(i == length(transFuncts)){
      cat("\nERROR: No transform can be estimated, exiting now\n")
      q()
    }
  }
}else if(args$transform == "arcsinh"){
  flowSet.trans <- read.flowSet.transVS(targets, file)
}else if(args$transform == "none"){
  flowSet.trans <- read.flowSet(file, transformation = F, truncate_max_range = F)
}else{
  cat("\nNo transform specified, exiting now\n")
  q()
}

colsToUse <- targets$name[targets$Lineage == 1 & targets$Ignore == 0]

if(length(colsToUse) %% 2  == 1){
  batches <- c(1, seq(4, length(colsToUse), by = 2))
}else{
  batches <- seq(1, length(colsToUse), by = 2)
}

flowTypeResults <- lapply(seq_along(flowSet.trans), function(y){
  PartitionsList <- lapply(seq_along(batches), function(x){
    start <- batches[x]
    if(is.na(batches[x+1])){
      end <- length(colsToUse)
    }else{
      end <- batches[x+1]-1
    }
    test <- flowType(flowSet.trans[[y]],
            PropMarkers = colsToUse[start:end],
            MFIMarkers = colsToUse[start:end],
            Methods = "kmeans",
            PartitionsPerMarker = targets$PartitionsPerMarker[targets$Ignore == 0][start:end],
            MarkerNames = targets$desc[colsToUse[start:end]])
    return(test@Partitions)
  })
  PartitionsList %>%
    do.call(cbind, .) %>%
  return(.)
})

names(flowTypeResults) <- file

dir.create(paste0(RESULTS_DIR, "PHENOTYPED_FCS/"),
           showWarnings = F)
for( files in file){
  rawFCS <- read.FCS(files, transformation = F)
  phenotype_data <- flowTypeResults[[files]]
  colnames(phenotype_data) <- paste0("Phenotype_", targets$desc[targets$name %in% colsToUse])
  clusterFCS <- flowCore::cbind2(rawFCS, phenotype_data)
  row.names(pData(parameters(clusterFCS))) <- paste0("$P", c(1:nrow(pData(parameters(clusterFCS)))))
  out.fcs.file <- paste0(RESULTS_DIR, "PHENOTYPED_FCS/phenotyped_", basename(files))
  write.FCS(clusterFCS, out.fcs.file)
}