UVデータを一発でキレイにする

UVデータを一発でキレイにした話

核酸のUV melting実験を行なっています。
機械はShimadzu UV1700 spectrophotometer (Shimadzu, Kyoto, Japan) です。
txtで、下画像のようなファイルが1スペクトルにつき1つ出てきますが、これを一つ一つエクセルで形式指定して開いて、それぞれのサンプルデータからバッファーをひいて、そして核酸のアニーリングとメルティングを１つのプログラムでやっているので、メルティング部分だけにカットして、規格化して、、、というのがめちゃくちゃめんどくさい！
f:id:xkwnc284:20191105205432p:plain
これをボタン一つでキレイにして一つのcsvファイルにまとめるプログラムを書こうとおもいます。

最終的にフィッティングもやりたかったので、Rで書き始めました。
途中でpythonに書き換えたかったのですが、うまいことできず挫折。
とりあえずRで動かしている。

bufferファイル1つとサンプルデータ２〜７つを同じフォルダに移す。

#グラフ一つにつき一つのファイル
#ディレクトリのファイル名をリスト化
  files <- list.files()
  for(file.name in files){
    if (grepl(".txt", file.name)==TRUE){
      if (grepl("buffer", file.name) == TRUE){
      buffer <- read.table(file.name, skip=95, header=F, sep=",")   
      colnames(buffer) <- c("Temp", "buffer")
      #print(head(buffer))
      #print(file.name)
      next
    }
    }
  }

ファイルネームを一つずつとってきてリスト化してfor文をまわす。ファイル名にbufferとあるファイルをbufferファイルとして読み込み、まずはバッファーのみのデータフレームを作成。

もういっちょfor文をまわして、次はbufferがファイル名に入ってないものを読み込んでいく。

#次のfor文で、bufferの文字のないファイルをデータファイルとしてよみこむ。
  options(digits=8) #出てくる結果の数字の桁数を設定
  df <- data.frame(buffer[,1])
  #temperatureだけ格納したデータフレーム、dfを準備。ここにデータを足していく
  for (file.name in files){
    if (grepl(".txt", file.name)==TRUE){
      if (grepl("buffer", file.name) == F){
        row <- read.table(file.name, skip=95, header=F, sep=",")
        sub <- row[,2]-buffer[,2]
        norm1 <- sub - min(sub)
        norm2 <- norm1/max(norm1)    #データの規格化
        
        name <- sub("-295.txt","",file.name) #ファイル名から.txtを除いたものをnameに入れる→グラフ表記やfigureのファイル名などに使う
       
        
        df <- cbind(df, norm2)
        #df <- data.frame(buffer[,1], norm2)    #x軸とy軸のデータを統合してdfに格納
        
        names(df)[ which( names(df)=="norm2" ) ] <- name
        #列名がnorm2のところを元のファイル名から取ってきたnameに変える
        #print(head(df))
   
      
        next
        }
  }
  }

ここまでで出来た１つのdfをresultファイルとして出力。

write.table(df, file="result.txt", quote=F, append=F, col.names=T, row.names=F)

ここまでのコードをプログラム名"before_plot2"として下記の{}にいれてやる

before_plot2 <- function(){ここに上記のコード}

これでコンソール上で

> source('~/R/before_plot2.R') #RStudioで、Source on Saveにチェックいれとけば、コードを書いたファイルを保存したときに自動的に実行される。
> before_plot2()

とすれば同じディレクトリにバッファーが引かれて、規格化されたデータが全て人つのファイルにまとまったものが出てくる。（下図）
シアワセ。

f:id:xkwnc284:20191106212045j:plain

まだまだ効率化の余地はありそうですが、とりあえずはここで満足。
ディレクトリをセットするのを忘れがちなので注意。

(おまけ）出てくるデータを非線形回帰によってフィッティングする

一気にやろうとすると、途中でエラーでひっかかったときに結局止まったところからやりなおしたり、初期値をいじいじしたりと逆に手間なので、引数にbufferのファイル名とサンプルのファイル名をとって、指定してやっていくことにした。

#bufferのファイル名とデータのファイル名を引数に
fitting <- function(buffer.file.name, data.file.name){
  buffer <- read.table(buffer.file.name, skip=95, header=F, sep=",")    

row <- read.table(data.file.name, skip=95, header=F, sep=",")
sub <- row[,2]-buffer[,2]
norm1 <- sub - min(sub)
#write.table(norm1, file="output_norm1.txt", quote=F, append=T)
norm2 <- norm1/max(norm1)    #データの規格化
#write.table(norm2, file="output_norm2.txt", quote=F, append=T)
name <- sub("-295.txt","",data.file.name)    #ファイル名から.txtを除いたものをnameに入れる→グラフ表記やfigureのファイル名などに使う

pdf(paste(name,".pdf"))  #プロット用デバイス作製

df <- data.frame(buffer[,1], norm2)  #x軸とy軸のデータを統合してdfに格納

options(digits=8)
par(mar=c(5.5, 6.0, 4.1, 2)) #余白の広さを行数で指定
par(mgp = c(4, 1.2, 0)) #余白の使い方、説明、ラベル、軸の位置を行で指定
par(ps=20)

#write.table(df, file=paste(name,"_df.txt", quote=F, append=T))
plot(df, main= name, xlab="Temperature (ºC)", ylab="Normalized absorbance at 295 nm", cex.lab = 1.2)

fitting <- nls(norm2 ~ (m1*(buffer[,1]+273.15)+m2)+((m3-m1)*(buffer[,1]+273.15)+(m4-m2))*(1/(1+1/(exp((m5/1.987)-(m6)/(1.987*(buffer[,1]+273.15)))))), df, start=list(m1=-0.01, m2=1, m3=0.0001, m4=0.34, m5=91, m6=30440))
lines(buffer[,1], fitted(fitting), col=2) #フィッティングラインを書く
df_2 <- data.frame(buffer[,1], norm2, fitted(fitting)) #フィッティングカーブのデータもデーターフレームに格納
write.table(df_2, file=paste(name,"_result.txt"), quote=F, append=F, col.names=T) #それをresultファイルに出力。カレイダなどに持って行って描画したいときのため
print(fitting)
print(summary(fitting))
dev.off()
}