R/UVfitting.R

In yanxianUCSB/yxhelper: Yx's Personal R packages

#!/usr/bin/Rscript

# devtools::install_github('yanxianUCSB/myRCodes')

UVfit <- function() {
    cat('------------UVfit V0.1------------\n')
    wd <- getwd()
    # baseline
    baseline <- get.abs(wd, 'baseline')
    # components
    protein <- get.abs(wd, 'protein')
    nucleicacid <- get.abs(wd, 'nucleic acid')
    if (is.null(protein) || is.null(nucleicacid))
        return(NULL)
    protein.c <-
        as.numeric(as.character(readline('protein conc [uM] = ')))
    nucleicacid.c <-
        as.numeric((as.character(
            readline('nucleic acid conc [ng/uL] = ')
        )))
    # baseline correction
    if (is.null(baseline))
        baseline <- rep(0, length(protein))
    protein <- protein - baseline
    nucleicacid <- nucleicacid - baseline

    for (i in 1:1000) {
        if (readline('Next? [y]/n\n') %in% c('', 'y')) {
            # spec to fit
            spc <- get.abs(wd, 'sample')
            spc <- spc - baseline
            # components normalization
            dat <- data.frame(
                spc = spc,
                protein = protein / protein.c,
                nucleicacid = nucleicacid / nucleicacid.c
            )
            # optimization
            result <- optim(
                par = c(protein.c, nucleicacid.c),
                two.component,
                data = dat
            )
            cat(
                sprintf(
                    "[Protein] = %.2f uM,\n[NA] = %.2f ng/uL.\n",
                    result$par[1],
                    result$par[2]
                )
            )
            cat(sprintf("RMSE = %.2f", result$value[1]))
        } else
            break
    }
}
two.component <- function(data, par) {
    x <- par[1]
    y <- par[2]
    spc <- data$spc
    protein <- data$protein
    nucleicacid <- data$nucleicacid
    return(sqrt(sum(abs(
        spc - x * protein - y * nucleicacid
    )) ^ 2 / length(spc)))
}
get.abs <- function(path, prompt = '') {
    cat(paste0('Select ', prompt, ' spec, Enter 0 to skip: '))
    spec <- select.list(dir(path, '.txt'), multiple = F)
    if (spec == "") {
        cat(paste0(prompt, ' spec skipped\n'))
        return(NULL)
    } else {
        spec <- read.delim(paste0(path, '\\', spec),
                           header = T,
                           as.is = T,
                           skip = 1)
        return(approx(
            x = spec$Wavelength.nm.,
            y = spec$Abs.,
            xout = seq(240, 320, 1)
        )$y)
    }
}
# UVfit()
# abs <- function(path) {
#     ds <-
#         read.delim(
#             select.list(dir(path, '.txt'), multiple = F),
#             header = T,
#             as.is = T,
#             skip = 1
#         )
#
#     abs <- c()
#     wl = c(260, 274, 320)
#     for(i in 1:3) {
#         abs[i] <- approx(ds$Wavelength.nm., ds$Abs., wl[i])$y
#     }
#     return(abs)
# }