R/CytoTron.R

In qbarbier/CytoTron: Cytometry Perceptron annotation

#` @export
CytoTron_run <- function(){
  appDir <- system.file("shinyApp", "app", package = "CytoTron")
  if (appDir == "")
  {
    stop("Could not find app directory. Try re-installing `CytoTron`.", call. = FALSE)
  }
  shiny::runApp(appDir, display.mode = "normal", launch.browser = T)
}

#` @export
overClusteringData <- function(fcs,markers,methode="kmeans",args=list(200,1000)){

  data <- fcs@exprs[,markers]
  table <- NULL
  clusters <- NULL
  if(methode == "kmeans" && length(args == 2)){
    km <- kmeans(data, centers=args[[1]],iter.max = args[[2]])
    table <- km$centers
    clusters <- km$cluster
  }
  if(methode == "FlowSOM" && length(args == 2)){
    fcs <- flowFrame(data)
    fs <- FlowSOM::BuildSOM(ReadInput(fcs,compensate = F,transform = F),
        xdim=args[[1]], ydim=args[[2]], colsToUse=c(1:dim(fcs)[2]))
    mst <- BuildMST(fs)
    clusters <- mst$map$mapping[,1]
    table <- mst$map$medianValues
  }
  if(methode == "CLARA" && length(args==3)){
    cl <- clara(data, k = args[[1]], samples=args[[2]], sampsize = args[[3]])
    table <- cl$medoids
    clusters <- cl$clustering
  }
  return(list(table, clusters))
}

#` @export
reduDimData <- function(fcs, markers, methode="PCA",args=NULL){

  data <- fcs@exprs[,markers]
  new.dim <- list(NULL)
  if(methode=="PCA"){
    new.dim <- prcomp(data)$x[,1:2]
  }
  if(methode == "UMAP"){
    new.dim <- uwot::umap(data)
    colnames(new.dim) <- c("UMAP1","UMAP2")
  }
  if(methode == "tSNE"){
    new.dim <- Rtsne(data,perplexity = args[[1]],theta = args[[2]], max_iter = args[[3]])$Y
    colnames(new.dim) <- c("TSNE1","TSNE2")
  }
  if(methode=="EmbedSOM"){
    fcs <- flowFrame(fcs@exprs[,markers])
    print(args[[3]])
    fs <- FlowSOM::BuildSOM(ReadInput(fcs,compensate = F,transform = F),rlen=args[[3]],
                            xdim=args[[1]], ydim=args[[2]], colsToUse=c(1:dim(fcs)[2]))
    new.dim <- data.frame(EmbedSOM::EmbedSOM(data=fs$data,map=fs$map))
    colnames(new.dim) <- c("EmbedSOM1","EmbedSOM2")
  }
  return(new.dim)
}

#` @export
overClusteringTrainData <- function(flow.frames,markers,methode="kmeans",args=list(200)){

  th <- args[[1]]
  if(methode == "FlowSOM"){th <- args[[1]]*args[[2]]}
  cluster <- c()
  res <- lapply(c(1:length(flow.frames)),function(fcs){
    data <- flow.frames[[fcs]]@exprs[,markers]
    print(fcs)
    if(dim(data)[1]<th || is.null(dim(data)[1])){
      if(is.null(dim(data)[1])){
        popName <- names(flow.frames)[[fcs]]
        table <- matrix(data,nrow=1)
        colnames(table) <- markers
        table<-data.frame(cbind(table,popName),check.names = FALSE)
      } else {
        x <- dim(data)[1]
        popName <- rep(names(flow.frames)[[fcs]],x)
        table <- cbind(data,popName)
      }
      return(table)
    }
    if(methode == "kmeans"){
      km <- kmeans(data,centers = args[[1]],iter.max = 1000)
      table <- km$centers
      cluster <<- c(cluster, km$cluster)
    }
    if(methode == "FlowSOM"){
      ff <- flow.frames[[fcs]][,markers]
      fs <- FlowSOM::BuildSOM(ReadInput(ff,compensate=F,transform=F),rlen=args[[3]],
                              xdim=args[[1]], ydim=args[[2]], colsToUse=c(1:dim(ff)[2]))
      mst <- BuildMST(fs)
      table <- mst$map$medianValues
    }
    if(methode == "CLARA"){
      clara <- clara(data, samples=50,k=args[[1]])
      table <- clara$medoids
    }
    popName <- rep(names(flow.frames)[[fcs]],args[[1]])
    table <- cbind(table,popName)
    return(table)
  })

  trainData <- do.call(rbind, res)
  return(trainData)
}

#` @export
ceilTrainData <- function(flow.frames, markers,replace=T,size=1000){

  res <- lapply(c(1:length(flow.frames)), function(i){
    mat <- flow.frames[[i]]@exprs[,markers]
    if(is.null(dim(mat))){
      table <- flow.frames[[i]]@exprs[rep(1,1000),markers]
    } else {
      table <- mat[sample(1:nrow(mat),replace=replace,size=size),]
    }
    popName <- rep(names(flow.frames)[[i]],size)
    table <- cbind(table,popName)
    return(table)
  })

  trainData <- do.call(rbind, res)
  return(trainData)

}

#` @export
createModel <- function(train,nameDim,hidden=c(5),func=c("relu"),epochs =50,
                        compile_function = "categorical_crossentropy",batch_size=50){

  other <- colnames(train)[which(colnames(train)!= nameDim)]
  #df <- train
  df <- data.frame(train, stringsAsFactors = F,row.names=c(1:dim(train)[1]),check.names = FALSE)
  for(i in other){
    print(i)
    v <- df[,i]
    v <- as.numeric(as.character(v))
    v[which(is.na(v))] <- 0
    df[,i] <- v
  }

  Train_Features <- data.matrix(df[,-dim(df)[2]])
  Train_Labels <- df[,nameDim]
  Train_Labels <- keras::to_categorical(as.numeric(as.factor(Train_Labels)))
  Train_Labels <- Train_Labels[,2:dim(Train_Labels)[2]]


  model <- keras::keras_model_sequential()
  for(i in c(1:(length(hidden)+1))){
    print(i)
    if(i == 1){
      print("a")
      keras::layer_dense(model,input_shape=length(other),
                         units=hidden[[i]],activation=func[[i]])
    }else if((i)>length(hidden)){
      print("c")
      keras::layer_dense(model,units=length(unique(train[,nameDim])),
                         activation="softmax")
    }else{
      print("b")
      keras::layer_dense(model,units=hidden[[i]],activation=func[[i]])
    }
  }

  keras::compile(model,loss=compile_function,optimizer="adam",metrics="accuracy")

  history <- keras::fit(model, x=Train_Features, y=Train_Labels,
                     epochs = epochs, batch_size = batch_size ,verbose=0)


  res <- evaluate(model, Train_Features, Train_Labels, verbose=0)
  list <- list("model"=model,"history"=history,"evaluate"=res)

  return(list)
}

#` @export
annotClusterFCS <- function(fcs, cluster, annotation){
  annotRaw <- annotation[cluster]
  fcs <- enrich.FCS(fcs, annotRaw, "CytoTron")
  return(fcs)
}

#` @export
getParamsCytoTron <- function(i){
  if(is.null(i)) return(NULL)
  data <- as.matrix(pData(i@parameters),stringsAsFactors = F)
  labels <- as.vector(unlist(data[,2]))
  params <- as.vector(unlist(data[,1]))
  if(length(which(is.na(labels)))>0){
    labels[which(is.na(labels))] <- params[c(which(is.na(labels)))]
  }
  if(length(which(labels=="<NA>"))>0){
    labels[which(labels=="<NA>")] <-params[c(which(labels=="<NA>"))]
  }
  if(length(which(labels=="NA"))>0){
    labels[which(labels=="NA")] <- params[c(which(labels=="<NA>"))]
  }
  names(params) <- labels
  return(params)
}

#` @export
found.spill.FCS <- function(fcs){
  comp <- "NULL"
  dim <- "NULL"
  desc <- fcs@description
  id1 <- unlist(lapply(desc, function(i){return(dim(i)[1])}))
  if(!is.null(id1)){
    comp <- names(id1)
    s <- fcs@description[[names(id1)]]
    if(dim(s)[1]==dim(s)[2]){
      dim <- dim(s)[1]
    }
  }
  return(c(comp,dim))
}

#` @export
compensate.FCS <- function(flow.frame, spill=NULL){
  if(is.null(spill) || spill == "NULL"){
    spill <- found.spill.FCS(flow.frame)[[1]]
    if(spill=="NULL"){
      warning("No compensation apply and/or found")
      return(flow.frame)
    }
  }
  fcs <- flowCore::compensate(flow.frame, flow.frame@description[[spill]])
  return(fcs)
}

#` @export
delete.column.FCS <- function(fcs, marker, spill=NULL){
  id <- which(colnames(fcs@exprs)==marker)
  print(id)
  mat <- exprs(fcs)
  mat <- mat[,-id]
  exprs(fcs) <- data.matrix(mat)
  desc <- fcs@description
  new.desc <- list()
  new.name <- list()
  for(i in c(1:length(desc))){
    if(length(grep(paste0("P",id),names(desc[i])))>0){
      print(names(desc)[i])
    } else {
      new.name <- c(new.name,names(desc)[i])
      new.desc <- c(new.desc,desc[i])
    }
  }
  names(new.desc) <- new.name
  fcs@description <- new.desc
  if(!is.null(spill)){
    new.spill <- fcs@description[[spill]]
    id.spill <- which(colnames(new.spill)==marker)
    if(length(id.spill)>0){
      new.spill <- new.spill[-id.spill,-id.spill]
      fcs@description[[spill]] <- new.spill
    }
  }
  return(fcs)
}

#` @export
enrich.FCS <- function(original, new.column, nw.names=NULL){
  new_p <- parameters(original)[1,]
  ## Now, let's change it's name from $P1 to $P26 (or whatever the next new number is)
  new_p_number <- as.integer(dim(original)[2]+1)
  rownames(new_p) <- c(paste0("$P", new_p_number))

  ## Now, let's combine the original parameter with the new parameter
  library('BiocGenerics') ## for the combine function
  allPars <-  BiocGenerics::combine(parameters(original), new_p)

  ## Fix the name and description of the newly added parameter, say we want to be calling it cluster_id

  if(is.null(nw.names)){
    new_p_name <- "cluster"
  } else {
    new_p_name <- nw.names
  }

  allPars@data$name[new_p_number] <- new_p_name
  allPars@data$desc[new_p_number] <- new_p_name

  new_exprs <- cbind(original@exprs, new.column)
  colnames(new_exprs) <- c(colnames(original@exprs),new_p_name)

  new_kw <- original@description
  new_kw["$PAR"] <- as.character(new_p_number)
  new_kw[paste0("$P",as.character(new_p_number),"N")] <- new_p_name
  new_kw[paste0("$P",as.character(new_p_number),"S")] <- new_p_name
  new_kw[paste0("$P",as.character(new_p_number),"E")] <- "0,0"
  new_kw[paste0("$P",as.character(new_p_number),"G")] <- "1"
  new_kw[paste0("$P",as.character(new_p_number),"B")] <- new_kw["$P1B"]
  new_kw[paste0("$P",as.character(new_p_number),"R")] <- new_kw["$P1R"]
  new_kw[paste0("flowCore_$P",as.character(new_p_number),"Rmin")] <- new_kw["flowCore_$P1Rmin"]
  new_kw[paste0("flowCore_$P",as.character(new_p_number),"Rmax")] <- new_kw["flowCore_$P1Rmax"]

  ## Now, let's just combine it into a new flowFrame
  new_fcs <- new("flowFrame", exprs=new_exprs, parameters=allPars, description=new_kw)

  return(new_fcs)
}

#` @export
logicle.FCS <- function(flow.frame, value=NULL, markers=NULL){

  if(is.null(markers)){
    if(is.null(flow.frame@description[[found.spill.FCS(flow.frame)[[1]]]])){
      markers.transform <- colnames(flow.frame)
    } else {
      markers.transform <- colnames(flow.frame@description[[found.spill.FCS(flow.frame)[[1]]]])
    }
  } else {
    markers.transform <- markers
  }

  list.index <- names(unlist(lapply(markers.transform, function(x) return(which(flow.frame@description==x)))))
  list.index <- gsub("N","", list.index)
  list.index <- gsub("\\$P","", list.index)

  if(is.null(value)){
    if(!is.null(flow.frame@description[[paste0("$P",list.index[1],"MS")]])){
      r.values <- unlist(lapply(list.index, function(x)
        as.integer(flow.frame@description[[paste0("$P",x,"MS")]]))
      )
    } else if(!is.null(flow.frame@description[[paste0("P",list.index[1],"MS")]])) {
      r.values <- unlist(lapply(list.index, function(x)
        as.integer(flow.frame@description[[paste0("P",x,"MS")]]))
      )
    } else {
      r.values <- rep(90, length(list.index))
    }
  }
  else {
    r.values <- rep(value, length(list.index))
  }

  w.values <- (4.5-log10(262143/abs(r.values)))/2
  w.values[which(w.values<0)] <- 0.5
  w.values[which(is.infinite(w.values))] <- 0.5
  for(t in 1:length(markers.transform)){
    lgcl <- flowCore::logicleTransform(w=w.values[t])
    flow.frame <- flowCore::transform(flow.frame, transformList(markers.transform[t],lgcl))
  }
  return(flow.frame)
}

#` @export
invers.logicle.FCS <- function(flow.frame, value=NULL, markers=NULL){
  if(is.null(markers)){
    if(is.null(flow.frame@description[[found.spill.FCS(flow.frame)[[1]]]])){
      markers.transform <- colnames(flow.frame)
    } else {
      markers.transform <- colnames(flow.frame@description[[found.spill.FCS(flow.frame)[[1]]]])
    }
  } else {
    markers.transform <- markers
  }

  list.index <- names(unlist(lapply(markers.transform, function(x) return(which(flow.frame@description==x)))))
  list.index <- gsub("N","", list.index)
  list.index <- gsub("\\$P","", list.index)

  if(is.null(value)){
    if(!is.null(flow.frame@description[[paste0("$P",list.index[1],"MS")]])) {
      r.values <- unlist(lapply(list.index, function(x)
        as.integer(flow.frame@description[[paste0("$P",x,"MS")]]))
      )
    } else if(!is.null(flow.frame@description[[paste0("P",list.index[1],"MS")]])) {
      r.values <- unlist(lapply(list.index, function(x)
        as.integer(flow.frame@description[[paste0("P",x,"MS")]]))
      )
    } else {
      r.values <- rep(90, length(list.index))
    }
  }
  else {
    r.values <- rep(value, length(list.index))
  }

  w.values <- (4.5-log10(262144/abs(r.values)))/2
  w.values[which(w.values<0)] <- 0.5
  w.values[which(is.infinite(w.values))] <- 0.5

  flow.frame.inv <- flow.frame

  for(t in 1:length(markers.transform)){
    invLgcl <- flowCore::inverseLogicleTransform(trans = logicleTransform(w=w.values[t]))
    flow.frame.inv <- flowCore::transform(flow.frame.inv, transformList(markers.transform[t],invLgcl))
  }

  return(flow.frame.inv)
}

#` @export
arcsinh.FCS <- function(flow.frame, marker=NULL, arg){
  raw <- flow.frame@exprs
  mat <- flow.frame@exprs
  if(is.null(marker) || length(marker)<1){
    marker <- colnames(flow.frame)
  }
  mat <- mat[,marker]
  colnames(mat) <- marker
  mat <- asinh(mat/arg)
  raw[,marker] <- mat[,marker]
  flow.frame@exprs <- raw
  return(flow.frame)
}

#` @export
invers.arcsinh.FCS <- function(flow.frame, marker_untrans, arg){
  raw <- flow.frame@exprs
  mat <- flow.frame@exprs
  marker_untrans_index <- which(colnames(flow.frame)%in%marker_untrans)
  mat <- mat[,-marker_untrans_index]
  marker <- colnames(mat)
  mat <- sinh(mat)*arg
  raw[,marker] <- mat[,marker]
  flow.frame@exprs <- raw
  return(flow.frame)
}

#` @export
decompensate.FCS <- function(x, spillover=NULL) {
  if(!is.null(spillover)){
    cols <- colnames(spillover)
    sel <- cols %in% colnames(x)
    e <- exprs(x)
    e[, cols] <- e[, cols] %*% spillover
    exprs(x) = e
    return(x)
  } else {
    return(x)
  }
}

#` @export
concatenate.FCS <- function(flow.frames, params="Flag"){
  ff.concat <- NULL
  n <- length(flow.frames)
  for(i in 1:n){
    ff.raw <- flow.frames[[i]]
    p <- matrix(i, nrow = nrow(ff.raw), ncol=1, dimnames = list(NULL, params))
    new.col <- as.vector(p)
    ff.raw <- enrich.FCS(ff.raw, new.col, nw.names=params)
    if(is.null(ff.concat)){
      ff.concat  <- ff.raw
    } else {
      exprs(ff.concat) <- rbind(exprs(ff.concat),exprs(ff.raw))
    }
  }
  return(ff.concat)
}

#` @export
deconcatenate.FCS <- function(data, params){
  print(params%in%colnames(data@exprs))
  if(params%in%colnames(data@exprs)){
    flow.frames <- lapply(sort(unique(unlist(data@exprs[,params]))), function(i){
      fcs <- data
      fcs@exprs <- data@exprs[which(data@exprs[,params]==i),]
      fcs <- delete.column.FCS(fcs,marker=params,spill=NULL)
      return(fcs)
    })
  } else {
    warning("Params does'nt exist in flowFrame files")
    return(fcs)
  }
  return(flow.frames)
}

#` @export
write.Label.Enrich.FCS <- function(fcs, annotation.column, populations.dataframe){
  #populations : >id1;pop1> ... >idN;popN

  #populations.dataframe: noms de colonne = ID, Name

  #Number of populations
  keyword.data <- ""

  #Populations
  for(i in 1:nrow(populations.dataframe))
  {
    pop <- populations.dataframe[i,]
    tmp.data <- paste(pop$Name, pop$ID, sep=";")
    keyword.data <- paste0(keyword.data, tmp.data, ">")
  }

  fcs@description[[paste0("P",annotation.column,"PopN")]] <- keyword.data

  return(fcs)
}

#` @export
read.Label.Enrich.FCS <- function(fcs, annotation.column, add.pop.size=T){
  pop.table <- NULL
  pop.keyword <- fcs@description[[paste0("P",annotation.column,"PopN")]]
  if(!is.na(pop.keyword) && length(pop.keyword)>0)
  {
    tmp.populations <- unlist(strsplit(pop.keyword, ">", fixed = T))[-1]
    pop.table <- matrix(0, ncol=3, nrow=length(tmp.populations))
    colnames(pop.table) <- c("ID", "Name", "Events")
    for(i in 1:length(tmp.populations))
    {
      tmp.pop <- tmp.populations[[i]]
      pop.table[i,c(1,2)] <- as.character(unlist(strsplit(tmp.pop, ";", fixed = T)))
      if(add.pop.size)
      {
        pop.table[i,3] <- as.integer(sum(fcs@exprs[,annotation.column]==as.integer(pop.table[i,1])))
      }
    }
    pop.table <- data.frame(pop.table, stringsAsFactors = F)
    pop.table$Events <- as.integer(pop.table$Events)
    if(!add.pop.size)
    {
      pop.table <- pop.table[,c(1,2)]
    }
  }
  return(pop.table)
}

#` @export
norm.percentile.FCS <- function(fcs, min.value=0.5, max.value=4.5, markers=NULL){
  if(is.null(markers)){
    markers <- colnames(fcs@description[[found.spill.FCS(fcs)[1]]])
  }
  if(is.null(markers)){
    warning("No markers found !")
    return(fcs)
  }
  for(i in markers){
    fcs@exprs[,i] <- (fcs@exprs[,i]-min.value)/(max.value-min.value)
  }
  return(fcs)

}

#` @export
unNorm.percentile.FCS <- function(fcs, min.value=0.5, max.value=4.5, markers=NULL){
  if(is.null(markers)){
    markers <- colnames(fcs@description[[found.spill.FCS(fcs)[1]]])
  }
  if(is.null(markers)){
    warning("No markers found !")
    return(fcs)
  }
  for(i in markers){
    fcs@exprs[,i] <- (fcs@exprs[,i]*(max.value-min.value))+min.value
  }
  return(fcs)
}

#` @export
clean.tails.FCS <- function(fcs, markers=NULL, zero=TRUE,max=TRUE){
  for(i in markers){
    if(zero){
      fcs <- fcs[which(fcs@exprs[,i]>0),]
    }
    if(max){
      m <- max(fcs@exprs[,i])
      fcs <- fcs[which(fcs@exprs[,i]<m),]
    }
  }
  return(fcs)
}