Process Fwave recorded by Neuropack

process_fwave.R

#%% function
load_raw <- function(fname_raw){
    samp_time <- 0.05 # ms
    line_ini <-0;line_end <-0
    linesRaw <- readLines(fname_raw)
    for (ii in seq(1, length(linesRaw))){
        if (length(grep(pattern='"1-1"', linesRaw[ii]))>0){
            line_ini <- ii
        }
        if (length(grep(pattern="Wave End", linesRaw[ii]))>0){
            line_end <- ii
            break
        }
    }
    raw_csv <- read.table(fname_raw, sep=',', header=F, skip=line_ini, nrows=line_end-line_ini-1)
    len_data <- dim(raw_csv)[1]
    n_trial <- dim(raw_csv)[2]
    times <- seq(0, length=len_data, by=samp_time)
    raw_data <- raw_csv[,2:n_trial]
    out <- list()
    out$times <- times
    out$raw <- matrix(as.matrix(raw_data), nrow(raw_data), ncol(raw_data))
    return (out)
}
get_nls_baseline <- function(epochs_norm, t_f_norm){
    windows();par(las=1, bty='l')
    matplot(t_f_norm, epochs_norm, type='l', lty=1, 
        main='click to set f wave segment\n :remove to calculate baseline by nls')
    t_lim_del <- c()
    for (ii in c(1,2)){
        loc <- locator(1)
        abline(v=loc$x)
        t_lim_del <- c(t_lim_del, loc$x)
    }
    epochs_norm_seg <- rbind(epochs_norm[(t_f_norm<=min(t_lim_del)),],
                        epochs_norm[(t_f_norm>=max(t_lim_del)),])
    t_f_norm_seg <- c(t_f_norm[(t_f_norm<=min(t_lim_del))], t_f_norm[t_f_norm>=max(t_lim_del)])
    list_rss_nls <- c()
    list_rss_poly <- c()
	list_rss_decay <- c()
    x <- t_f_norm_seg
    len_data <- dim(epochs_norm)[2]
	rss <- function(y,mdl){sum((y-fitted(mdl))^2)}
	mdl_nls <- function(x,y){nls(y~a - a/exp(b*x), start=c(a=max(y), b=max(x)/max(y)))}
	mdl_nls_decay <- function(x,y){nls(y~c/exp(a*x)*cos(b*x) + d, start=c(a=max(x)/max(y), b=pi/max(x), c=-tail(y,1), d=tail(y,1)))}
	mdl_poly <- function(x,y){lm(y~x+I(x^2)-1)}
    # try_get_rss <- function(x, y, mdl, mdl_name){
    #     rss <- function(y,mdl){sum((y-fitted(mdl))^2)}
    #     rss_try <- try(rss(y, mdl(x,y), silent=T))
    #     if (class(rss_try)=="try-error"){
    #         rss_try <- Inf
    #         print(paste("error in ", mdl_name, sep=''))
    #     }
    #     return(rss_try)
    # }
    for (ii in seq(1, len_data)){
        y <- epochs_norm_seg[, ii];
        # nls
        rss_nls <- try(rss(y, mdl_nls(x,y)), silent=T)
        if (class(rss_nls)=="try-error"){
            rss_nls <- Inf
            print("error in mdl_nls")
        }
        # rss_nls <- try_get_rss(x, y, mdl_nls, "mdl_nls")
        list_rss_nls <- c(list_rss_nls, rss_nls)
		# nls: decay curve
        rss_decay <- try(rss(y, mdl_nls_decay(x,y)), silent=T)
        if (class(rss_decay)=="try-error"){
            rss_decay <- Inf
            print("error in mdl_nls_decay")
        }
        # rss_decay <- try_get_rss(x, y, mdl_nls_decay, "mdl_nls_decay")
        list_rss_decay <- c(list_rss_decay, rss_decay)
        # polynomial2
        rss_poly <- try(rss(y, mdl_poly(x,y)), silent=T)
        if (class(rss_poly)=="try-error"){
            rss_poly <- Inf
            print("error in mdl_poly")
        }
        # rss_poly <- try_get_rss(x, y, mdl_poly, "mdl_poly")
        list_rss_poly <- c(list_rss_poly, rss_poly)
    }
	min_idx <- function(x){which(x==min(x))}
    list_rsss <- list(list_rss_nls, list_rss_decay, list_rss_poly)
    mdls <- c(mdl_nls, mdl_nls_decay, mdl_poly)
    min_rsss <- c(min(list_rss_nls), min(list_rss_decay), min(list_rss_poly))
    list_rss <- list_rsss[min_rsss == min(min_rsss)]
    mdl_min_rss <- mdls[min_rsss == min(min_rsss)][[1]]
    idx_min_rss <- min_idx(list_rss[[1]])
    y <- epochs_norm_seg[, idx_min_rss]
    x <- t_f_norm_seg
    mdl <- mdl_min_rss(x,y)
    y_predict <- predict(mdl, list(x=t_f_norm))
    ttl <- paste("click somewhere when ok: fitted-index",idx_min_rss)
    lines(t_f_norm, y_predict, col=2, lwd=2)
    loc_end <- locator(1)
    if (length(loc_end)>0){
        return(y_predict)
    } else {
        return(0)
    }
}
nonlinear_fit_fwave <- function(epochs_fwave, times_fwave){
    # len_data <- dim(raw_data)[2]
	#t1 <- 18;t2 <- 50
	t_f_norm <- times_fwave - times_fwave[1]
    len_epoch <- dim(epochs_fwave)[1]
	epochs_norm <- epochs_fwave - matrix(rep(1, len_epoch), nrow=len_epoch)%*%epochs_fwave[1,]
    continue <- T
    while (continue){
        y_predict <- get_nls_baseline(epochs_norm, t_f_norm)
        if (length(y_predict)>1){continue<-F} else {dev.off()}
    }
    return(y_predict)
}

detrend_and_measure_famp <- function(raw_data, times, tlim_f){
    # run after plot_m_f
    # params:
    #   tlim_f: time range for fwave segment-min, max
    epochs_fwave <- raw_data[(times >= min(tlim_f)) & (times <= max(tlim_f)), ]
    times_fwave <- times[(times >= min(tlim_f)) & (times <= max(tlim_f))]
    trend_nl_line <- nonlinear_fit_fwave(epochs_fwave, times_fwave)
    epochs_fwave_nls_detrend <- epochs_fwave - 
		matrix(trend_nl_line, nrow=length(trend_nl_line))%*%matrix(rep(1,ncol(epochs_fwave)),nrow=1)
    # add lines to f wave
	dev.set(dev.list()[1])
	lines(times_fwave, trend_nl_line + mean(epochs_fwave[1,]), col=2)
    # plot detrended fwave
    windows()
    par(bty='l',las=1, mfrow=c(2,1))
    matplot(times_fwave, epochs_fwave_nls_detrend, type='l', lty=1, lwd=1,
        xlab='time (ms)', ylab='amplitude (μV)', main='F wave detrended by nls')
    # locator
    xlims <- c()
    for (ii in c(1:2)){
        loc <- locator(1)
        abline(v=loc$x)
        xlims <- c(xlims, loc$x)
    }
    fwaves_seg <- epochs_fwave_nls_detrend[(times_fwave >= min(xlims)) & (times_fwave<=max(xlims)),]
    times_seg <- times_fwave[(times_fwave >= min(xlims)) & (times_fwave<=max(xlims))]
    fwaves_seg_detrend <- apply(fwaves_seg, 2, function(y){y-seq(y[1], tail(y,1), length=length(y))})
    famps <- apply(fwaves_seg_detrend, 2, function(y){max(y)-min(y)})
    matplot(times_seg, fwaves_seg_detrend, type='l', lty=1, 
        xlab='time (ms)', ylab='amplitude (μV)', main='F wave - linear detrended')
    return(famps)
}

plot_m_f <- function(raw_data, times, tlim_f, ttl){
    # info 
    xlim <- c(0, 60) # ms
    tmin_m <- 3 # ms
    step_y_f <- 200 # microV
    step_y_m <- 500 # microV
    # determine ylim for m wave
    raw_data_seg_m <- raw_data[(times>tmin_m),]
    min_seg_m <- floor(min(raw_data_seg_m/step_y_m))*step_y_m
    max_seg_m <- ceiling(max(raw_data_seg_m/step_y_m))*step_y_m
    ylim_m <- c(min_seg_m, max_seg_m)
    # determine ylim for fwave
    raw_data_seg_f <- raw_data[(times > min(tlim_f)) & (times< max(tlim_f)),]
    min_seg_f <- floor(min(raw_data_seg_f)/step_y_f)*step_y_f
    max_seg_f <- ceiling(max(raw_data_seg_f)/step_y_f)*step_y_f
    ylim_f <- c(min_seg_f, max_seg_f)
    windows()
    par(mfrow=c(1,2))
    par(bty='l',las=1)
    par(oma=c(0,0,4,0))
    matplot(times, raw_data, lty=1, lwd=2, ylab='amplitude (μV)', xlab='time (ms)',
        col=1, type='l', xlim=xlim, ylim=ylim_m, main='M wave')
    matplot(times, raw_data, lty=1, lwd=2, ylab='amplitude (μV)', xlab='time (ms)',
        col=1, type='l',ylim=ylim_f, xlim=xlim, main='F wave')
    mtext(text=ttl, outer=T)
    # abline
    t_seg <- c(0,0)
    for (ii in c(1,2)){
        loc <- locator(1)
        abline(v=loc$x, col = "blue")
        t_seg[ii]<-loc$x
    }
    out <- list()
    out$t_seg <- t_seg
    return(out)
}
#%% main
# load raw
path_rawcsv <- file.choose()
dir_raw <- dirname(path_rawcsv)
setwd(dir_raw)
raw_med_rest <- load_raw(path_rawcsv)
#%% plot
# info
tlim_f <- c(18,60) # ms to determine fwave plot ylim
# plot
times <- raw_med_rest$times
raw_data <- raw_med_rest$raw
# if title is required use below:
# print("Title?");stdin_title <- readLines(stdin(), n=1)
stdin_title="M wave and F wave"
latencies <- plot_m_f(raw_data, times, tlim_f, ttl=stdin_title)
# detrend and measure fwave amps
tlim_f <- latencies$t_seg
famps <- detrend_and_measure_famp(raw_data, times, tlim_f)
# save plot CMAP+F wave
dir_desktop <- file.path(Sys.getenv("USERPROFILE"),"Desktop")
# dir_fig <- paste(dir_desktop, 'botox_fwave','pdf',sep='/')
dir_fig <- paste(dir_raw,'pdf',sep='/')
fname_cond = strsplit(basename(path_rawcsv),'_')[[1]]
fname_cond_header = paste(fname_cond[1], fname_cond[2],sep="_")
fname_fig_m_f <- paste("fig_m_f_", fname_cond_header, format(Sys.time(), "_%y%m%d%H%M%S"), ".pdf", sep="") 
path_fig_m_f <- paste(dir_fig, fname_fig_m_f, sep='/')
len_dev <- length(dev.list())
dev.set(dev.list()[len_dev-2])
dev.print(device=pdf, file=path_fig_m_f, width=12, height=8)
# save detrended F wave
fname_detrend <- paste("fig_fwave_detrend_",fname_cond_header, format(Sys.time(), "_%y%m%d%H%M%S"), ".pdf", sep="")
path_fig_fdetrend <- paste(dir_fig, fname_detrend, sep='/')
dev.set(dev.list()[len_dev])
dev.print(device=pdf, file=path_fig_fdetrend, width=8, height=10)
print(famps)
fname_csv_save=paste('fwave_amps_', fname_cond_header, format(Sys.time(), "_%y%m%d%H%M%S"),'.csv',sep='')
path_csv_save = paste(dir_fig, fname_csv_save,sep='/')
write.table(famps, file=path_csv_save, sep=',',row.names=F, col.names=F)
#%% stats

日本光電で記録したF波を処理するRスクリプト

• 非線形ぽい基線を補正する処理を組み込みました
• 時間を区切ってその間の基線を補正してF波の振幅を測定します
• 処理する際の波形などpdfにして出力

注意点

• csvと同じディレクトリに"pdf"というフォルダがあるのを前提にそこへF波振幅のcsvやpdfを保存します。"pdf"のフォルダがないとエラーで止まります。
• 基線の補正は時々きれいにおさまりません。3種類(2次関数、指数関数、指数関数のdecay曲線)で補正してroot sum squareが最小になる方法で補正をするような段取りにしています。