R/02_logicle_transforms.R

In ebecht/infinityFlow: Augmenting Massively Parallel Cytometry Experiments Using Multivariate Non-Linear Regressions

#' Transform FCS files
#' @description Backbone measurements use a common transformation across the whole dataset. Exploratory measurements are transformed well/file-wise.
#' @param yvar name of the exploratory measurement
#' @param paths Character vector of paths to store intput, intermediary results, outputs...
#' @param xp Logicle-transformed backbone expression matrix
#' @param chans vector of backbone channels' names
#' @param events.code vector of length nrow(xp) specifying from which well each event originates
#' @param annot annotation table generated by infinityFlow:::intialize()
#' @param verbose Verbosity
#' @importFrom utils read.table
#' @importFrom stats quantile
#' @importFrom flowCore logicleTransform inverseLogicleTransform
#' @noRd
logicle_transform_input <- function(
                                 yvar,
                                 paths,
                                 xp=readRDS(file.path(paths["rds"],"xp.Rds")),
                                 chans=readRDS(file.path(paths["rds"],"chans.Rds")),
                                 events.code=readRDS(file.path(paths["rds"],"pe.Rds")),
                                 annot=read.table(paths["annotation"],sep=",",header=TRUE,stringsAsFactors=FALSE),
                                 verbose=TRUE
                                 ){
    
    ## ##################
    ## Computing parameters for each channel for each project using the code from flowCore's estimateLogicle
    ## ##################
    if(verbose){
        message("Logicle-transforming the data")
        message("\tBackbone data")
    }
    
    transforms_chan <- setNames(
        lapply(
            chans,
            function(x){
                data <- xp[,x]
                t <- max(data)
                m <- 4.5
                q <- 0.05
                r <- .Machine$double.eps + quantile(data, q)
                w <- max((m-log10(t/abs(r)))/2,0.1)
                w <- min(w,m/2)
                a <- 0
                logicleTransform(w=w,t=t,m=m,a=a) ##Just use summary() to retrive the parameters
            }
        ),
        chans
    )

    if(verbose){
        message("\tExploratory data")
    }
    transforms_pe <- lapply(
        split(xp[,yvar], events.code),
        function(x){
            data <- x
            t <- max(data)
            m <- 4.5
            q <- 0.05
            r <- .Machine$double.eps + quantile(data, q)
            w <- max((m-log10(t/abs(r)))/2,0.1)
            w <- min(w,m/2)
            a <- 0
            logicleTransform(w=w,t=t,m=m,a=a)
        }
    )
    
    if(verbose){
        message("\tWriting to disk")
    }
    saveRDS(transforms_chan,file=file.path(paths["rds"],"transforms_chan.Rds"))
    saveRDS(transforms_pe,file=file.path(paths["rds"],"transforms_pe.Rds"))

    ## ##################
    ## Exporting transformed expression matrices
    ## ##################
    if(verbose){
        message("\tTransforming expression matrix")
    }
    
    for(chan in chans){
        xp[,chan] <- transforms_chan[[chan]](xp[,chan])
    }

    d.e <- split(as.data.frame(xp),events.code)
    d.e <- lapply(d.e,as.matrix)
    for(chan in unique(events.code)){
        d.e[[chan]][,yvar] <- transforms_pe[[chan]](d.e[[chan]][,yvar])
    }

    xp <- do.call(rbind,d.e)

    if(verbose){
        message("\tWriting to disk")
    }
    
    saveRDS(xp,file=file.path(paths["rds"],"xp_transformed.Rds"))
    invisible()
}