cutoff <- 10 # ここで指定した数値(in bp)より離れたピークは同一ピークとみなさない infile <- "data.txt" # 入力ファイル名を指定 outfile <- "data_matrix.txt" # 出力ファイル名を指定 ############################################# ### データの読み込み ### ############################################# data <- read.table(infile, header=TRUE, row.names=NULL, sep="\t") # データファイルの読み込み data <- cbind(data, data[,3]/(sqrt(2*pi)*data[,4])) # peak area(data[,3])とheight(data[,4])情報から正規分布の標準偏差情報を算出 ############################################# ### Complete-linkage clustering (step 1) ### ### まずは全く同じ数字のものをマージ ### ############################################# data_s <- data[order(data[,2]),] # M.W.値が昇順になるように並び替え cluster_num <- 1 # クラスターのシリアル番号、1からはじめる cluster <- rep(0, nrow(data_s)) # 初期値として0を与えたベクトルを作成 TDF <- rep(0, nrow(data_s)) # 初期値として0を与えたベクトルを作成 TDF_num <- 1 # TDFのシリアル番号、1からはじめる tmp_mw <- 0 i <- 1 while(i <= nrow(data_s)){ # 初期値1のiがdata_sの行数分に達するまで繰り返される(基本的には最初の行から最後の行に向かっていく感じ) if((data_s[i,2] > tmp_mw) && (length(data_s[data_s[,2] == data_s[i,2],2]) >= 2)){ # i行目のM.W.がtmp_mwよりも大きく、かつi行目のM.W.と同じ分子量をもつものが二つ以上ある場合 cluster[i:(i+length(data_s[data_s[,2] == data_s[i,2],2])-1)] <- cluster_num # i行目のM.W.と同じ分子量をもつ全ての行について同じクラスター番号cluster_numを付与 TDF[i:(i+length(data_s[data_s[,2] == data_s[i,2],2])-1)] <- TDF_num # i行目のM.W.と同じ分子量をもつ全ての行について同じTDF番号TDF_numを付与 i <- i + length(data_s[data_s[,2] == data_s[i,2],2]) # iの値を「現在のiの値 + i行目のM.W.と同じ分子量をもつ行数」にする cluster_num <- cluster_num + 1 # クラスター番号を1足す TDF_num <- TDF_num + 1 # TDF番号を1足す } else if(i < nrow(data_s)){ # 最初のifを満たさず、且つ一番最後の行でない場合 cluster[i] <- cluster_num # i行目のクラスター番号にcluster_numを付与 cluster_num <- cluster_num + 1 # クラスター番号を1足す i <- i + 1 # iを1足す } else{ # 最初のifを満たさず、且つ一番最後の行である場合 cluster[i] <- cluster_num # i行目のクラスター番号にcluster_numを付与 i <- i + 1 # iを1足す } } data_s <- cbind(data_s, cluster, TDF) ############################################# ### Complete-linkage clustering (step 2) ### ### それ以外について ### ############################################# tmp_cluster <- NULL tmp_TDF <- NULL while(length(unique(cluster)) > 1){ ############################ ### 全通りの距離を計算 ### ############################ tmp <- NULL tmp1 <- NULL for(j in 1:(length(unique(cluster))-1)){ tmp[j] <- max(data_s[cluster == j+1, 2]) - min(data_s[cluster == j, 2]) # Complete-linkageなので「クラスター間の距離が最も遠いもの同士」の距離として定義 tmp1[j] <- j # シリアル番号も保持 } kouji <- cbind(tmp, tmp1) # 列同士を結合 tmp_posi <- kouji[order(kouji[,1]),2][1] # 最も距離が近いクラスターのシリアル番号をtmp_posiに格納 ######################################################################################## ### 最も距離が近い二つのクラスターについて同じサンプル由来ピークがあるかどうか判定 ### ######################################################################################## if(length(intersect(as.vector(data_s[cluster == tmp_posi,1]), as.vector(data_s[cluster == tmp_posi+1,1]))) == 0){ # 同じサンプル由来ピークがない場合 if((TDF[cluster == tmp_posi][1] == 0) && (TDF[cluster == tmp_posi+1][1] == 0)){ # 二つのクラスターともにTDFのシリアル番号が割り当てられていない(not >=1)場合 if(max(data_s[cluster == tmp_posi+1,2]) - min(data_s[cluster == tmp_posi,2]) <= cutoff){ # クラスター間の距離がcutoff以下の場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi] <- TDF_num # M.W.が低い方のクラスターに同じシリアル番号(TDF_num)を付与 TDF[cluster == tmp_posi+1] <- TDF_num # M.W.が高い方のクラスターに同じシリアル番号(TDF_num)を付与 } else{ # クラスター間の距離がcutoffより大きい場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi] <- TDF_num # M.W.が低い方のクラスターにTDF_numのシリアル番号を付与 TDF_num <- TDF_num + 1 # M.W.が高い方のクラスターに(1足した)別のシリアル番号を用意 TDF[cluster == tmp_posi+1] <- TDF_num # M.W.が高い方のクラスターにTDF_numのシリアル番号を付与 } } else if(TDF[cluster == tmp_posi][1] == 0){ # M.W.が低い方のクラスターにTDFのシリアル番号が割り当てられていない(not >=1)場合 if(max(data_s[cluster == tmp_posi+1,2]) - min(data_s[cluster == tmp_posi,2]) <= cutoff){ # クラスター間の距離がcutoff以下の場合 TDF[cluster == tmp_posi] <- TDF[cluster == tmp_posi+1][1] # M.W.が高い方のクラスターに既に割り当てられているTDFのシリアル番号をM.W.が低い方のクラスターに付与 } else{ # クラスター間の距離がcutoffより大きい場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi] <- TDF_num # M.W.が低い方のクラスターに新たなTDFのシリアル番号を付与 } } else if(TDF[cluster == tmp_posi+1][1] == 0){ # M.W.が高い方のクラスターにTDFのシリアル番号が割り当てられていない(not >=1)場合 if(max(data_s[cluster == tmp_posi+1,2]) - min(data_s[cluster == tmp_posi,2]) <= cutoff){ # クラスター間の距離がcutoff以下の場合 TDF[cluster == tmp_posi+1] <- TDF[cluster == tmp_posi][1] # M.W.が低い方のクラスターに既に割り当てられているTDFのシリアル番号をM.W.が高い方のクラスターに付与 } else{ # クラスター間の距離がcutoffより大きい場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi+1] <- TDF_num # M.W.が高い方のクラスターに新たなTDFのシリアル番号を付与 } } else{ # 二つのクラスター両方に既にTDFのシリアル番号が割り当てられている(>=1)場合 if(max(data_s[cluster == tmp_posi+1,2]) - min(data_s[cluster == tmp_posi,2]) <= cutoff){ # クラスター間の距離がcutoff以下の場合 TDF[cluster == tmp_posi+1] <- TDF[cluster == tmp_posi][1] # M.W.が低い方のクラスターに既に割り当てられているTDFのシリアル番号をM.W.が高い方のクラスターに付与 } else{ # クラスター間の距離がcutoffより大きい場合(何もしない) } } cluster[cluster > tmp_posi] <- cluster[cluster > tmp_posi] - 1 # tmp_posiより大きいクラスター番号から1引く } else{ # 同じサンプル由来ピークがある場合(マージしない) if((TDF[cluster == tmp_posi][1] == 0) && (TDF[cluster == tmp_posi+1][1] == 0)){ # 二つのクラスターともにTDFのシリアル番号が割り当てられていない(not >=1)場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi] <- TDF_num # M.W.が低い方のクラスターにTDF_numのシリアル番号を付与 TDF_num <- TDF_num + 1 # M.W.が高い方のクラスターに(1足した)別のシリアル番号を用意 TDF[cluster == tmp_posi+1] <- TDF_num # M.W.が高い方のクラスターにTDF_numのシリアル番号を付与 } else if(TDF[cluster == tmp_posi][1] == 0){ # M.W.が低い方のクラスターにTDFのシリアル番号が割り当てられていない(not >=1)場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi] <- TDF_num # M.W.が低い方のクラスターにTDF_numのシリアル番号を付与 } else if(TDF[cluster == tmp_posi+1][1] == 0){ # M.W.が高い方のクラスターにTDFのシリアル番号が割り当てられていない(not >=1)場合 TDF_num <- TDF_num + 1 # 新たなTDFのシリアル番号として現在のシリアル番号の最大値に1を足した数値を用意 TDF[cluster == tmp_posi+1] <- TDF_num # M.W.が高い方のクラスターに新たなTDFのシリアル番号を付与 } else{ # 二つのクラスター両方に既にTDFのシリアル番号が割り当てられている(>=1)場合(何もしない) } cluster[cluster > tmp_posi] <- cluster[cluster > tmp_posi] - 1 # tmp_posiより大きいクラスター番号から1引く } tmp_cluster <- cbind(tmp_cluster, cluster) # 一つ一つクラスターにしていったクラスター番号の変遷を保存 tmp_TDF <- cbind(tmp_TDF, TDF) # TDFのシリアル番号の変遷を保存 } TDF <- TDF - 1 tmp <- cbind(data[order(data[,2]),] , TDF) # 割り振られたTDF番号情報を入力データの一番右列に結合してtmpに保存 #write.table(tmp, "data_result.txt", sep="\t", row.names=FALSE, col.names=TRUE, quote=FALSE) # tmpを出力ファイルに保存 ############################################# ### ピークアラインメント結果を行列にする ### ############################################# sample <- as.vector(unique(data[,1])) mat <- matrix(rep(0, length(unique(TDF))*length(sample)), nrow=length(unique(TDF))) colnames(mat) <- sample tmp_mw <- 0 i <- 1 tmp_num <- 1 while(i <= nrow(tmp)){ if((tmp[i,6] != tmp_mw) && (length(tmp[tmp[,6] == tmp[i,6],6]) >= 2)){ # i行目のM.W.がtmp_mwよりも大きく、かつi行目のM.W.と同じ分子量をもつものが二つ以上ある場合 mat[tmp_num, as.vector(tmp[tmp[,6] == tmp[i,6],1])] <- as.numeric(as.vector(tmp[tmp[,6] == tmp[i,6],2])) i <- i + length(tmp[tmp[,6] == tmp[i,6],6]) # iの値を「現在のiの値 + i行目のM.W.と同じ分子量をもつ行数」にする tmp_num <- tmp_num + 1 } else if(i < nrow(tmp)){ # 最初のifを満たさず、且つ一番最後の行でない場合 mat[tmp_num, as.vector(tmp[tmp[,6] == tmp[i,6],1])] <- as.numeric(as.vector(tmp[tmp[,6] == tmp[i,6],2])) i <- i + 1 tmp_num <- tmp_num + 1 } else{ # 最初のifを満たさず、且つ一番最後の行である場合 mat[tmp_num, as.vector(tmp[tmp[,6] == tmp[i,6],1])] <- as.numeric(as.vector(tmp[tmp[,6] == tmp[i,6],2])) i <- i + 1 } } tmp2 <- cbind(1:nrow(mat), mat) colnames(tmp2) <- c("TDF", sample) # 列名を与えている write.table(tmp2, outfile, sep="\t", row.names=FALSE, col.names=TRUE, quote=FALSE)