Characteristic Direction Robustnes.R

#' @param datain expression data with gene id in the first column
#' @param ncontrol number of samples from control
#' @param ntreatment number of samples from treatment
#' @param control indices of control samples
#' @param treatment indices of treatment samples
#'
random_chdir <- function(datain, ncontrol, ntreatment, control, treatment, randomize_gamma=FALSE) {
    # Take ncontrol samples form control
    control <- sample(control, ncontrol)
    # and ntreatment from treatment
    treatment <- sample(treatment, ntreatment)

    # Create sampleclass vector
    sampleclass <- factor(c(rep('1', ncontrol), rep('2', ntreatment)))

    # Choose gamma
    gamma <- ifelse(randomize_gamma, runif(1), 1)
    
    # Run chdir and extract chdir vector
    png(tempfile())
    chdir <- GeoDE::chdirAnalysis(
        datain=datain[, c(colnames(datain)[1], control, treatment)],
        sampleclass=sampleclass,
        gammas=gamma
    )$chdirprops$chdir[[1]]
    dev.off()    
    
    # Results
    list(
        chdir=chdir,
        control=control,
        treatment=treatment,
        gamma=gamma
    )
}



#' Compare ranks of genes in two chdir vectors
#'
#' @param chdir_a list
#' @param chdir_b list
#' @param nbreaks number of breaks for cut function
#'
chdir_check <- function(chdir_a, chdir_b, nbreaks=50, show_xlab=TRUE, show_ylab=TRUE, show_gammas=FALSE) {
    
    br <- seq(0, length(chdir_a$chdir), length.out = nbreaks)
    
    im <- table(
        cut(rank(chdir_a$chdir), breaks = br),
        cut(rank(chdir_b$chdir), breaks = br)
    )

    xlab <- paste(
        paste(chdir_a$control, collapse=' '),
        paste(chdir_a$treatment, collapse=' '),
        sep='\n'
    )
    
    ylab <- paste(
        paste(chdir_b$control, collapse=' '),
        paste(chdir_b$treatment, collapse=' '),
        sep='\n'
    )
  
    main <- ifelse(
        show_gammas,
        paste('Gammas: (', round(chdir_a$gamma, 3), ', ', round(chdir_b$gamma, 3), ')', sep=''),
        ''
    )

    image(im, col=gplots::rich.colors(64), axes=FALSE, main=main)
    if(show_xlab) { mtext(xlab, side=1, line=2, cex=0.50) }
    if(show_ylab) { mtext(ylab, side=2, line=1.1, cex=0.50) }
}