Wanderlustという手法がある(Cell. 2014 Apr 24;157(3):714-25.)。Trajectory detectionともいう。どちらで呼べばいいかはようわからん。
FACSデータのような高次元データに対して、SPADEやviSNEなどのように、適当にデータ量を削減しながら、分化系統樹が作成でき、多次元データを2次元にプロットできるようになる。
アルゴリズムの特徴としては、viSNEのように非線形な関係の多次元データに対して適応を考えている。また、分化系統樹はSPADE(これは線形仮定)のように枝分かれを仮定せず、一直線になるものと思っている。たぶん拡張版がいつか出ると思う。一直線の木を考えるため、N個のマーカーはそれぞれなめらかに発現量が変化するものと仮定する。
アルゴリズムの概略は、k-NN graphを作成してひとつなぎにして骨格的なものを想定し、グラフ上の適当なノードを選んでwaypointsと呼ばれる個の監視点
を置く。分化の最初の地点(early cellとかinitiatorとかいう)だけはユーザーが設定する必要があり、この点からの相対的な分化度が最終的なtrajectory scoreとなる。論文では未熟B細胞のマーカーであるCD34が最大の細胞を設定している。[0, 1]補正することで、0、つまりearly cellだが、これに近いと未分化、1に近いと高分化となる。
k-NN graphにおいて、本のエッジから適当に
本のエッジに間引いて、subgraphであるl-k-NN graphを作成する。すべての細胞
について、
から
までと、
からl-k-NN graphのwaypointsすべて(
個)の間の、コサイン距離(非線形仮定なので使っている。仮定次第ではユークリッド距離でもよさそう)に基づいた重み付けのあるノード間距離を計算する。
から
の距離だけでは、各細胞の相対的な位置関係が分からないとかなんとかで、waypointsが導入されているっぽい。あるwaypointsの組を決めて
からの距離をそれっぽく収束するまで計算する。これがとあるl-k-NN graphのtrajectory scoreになる。
l-k-NN graphの作り方は何通りもあるが、適当に個のパターンからtrajectory scoreを計算して、平均でも取れば最終的なtrajectory scoreになる。
論文っぽいデータをlocatorでポチポチして作成した。ここから,
,
,
くらいでk-NN graphを作成する。
l-k-NN graphのひとつとしてこんなのを作った。赤が自分で決めるearly cell 、灰色が
個のwaypointsである。赤のノードから始まって、クルッと42番のノードまで行けば分化終了、というイメージ。
このl-k-NN graphにおいて、すべてのノードとwaypointsの間の距離を計算すると、こんな感じになる。waypoints自身との距離は当然、0になる。
あるノードと、
もしくはwaypoints
との距離
において、
か
かで条件分けする。そうしたらとりあえずのtrajectory scoreが出る。
l-k-NN graphを適当に変化させて、回くらい繰り返すとこんな感じになる。平均を取れば赤線。これが最終的なtrajectory scoreになる。
trajectory scoreと各マーカーの関係を見ると、x軸のマーカーは分化最中にピークを迎えて後半にダラっと低下する。y軸のマーカーは最初からクライマックスであとは低下しっぱなし、という様子がtrajectory scoreが大きくなる、つまり、分化が進むに従ってわかる。
というわけで、Wanderlust: Trajectory scoreができるようになったので、定番のアイドルデータセットでやってみよう。はもっともロリだと思われる横山千佳にした。
ラブライブのメンバーだけ取ってきたが、結果は以下の通り。かわいいかしこいエリーチカとスピリチュアルでぶがBBAの方にくるのは予想通りだった。2年生組はど真ん中、まきちゃんだけお姉さんぶってちょっとBBA寄りだった。りんちゃんはわかるが、かよちんが穂乃果、海未ちゃんよりもロリ側だったのはちょっと意外だった。
(l-k-NN graphの発生はランダムなので結果が変わってしまった。)
age height weight B W H score 矢澤にこ 17 154 41 74 57 79 0.3002550 星空凛 15 155 43 75 59 80 0.3133763 高坂穂乃果 16 157 45 78 58 82 0.3260504 小泉花陽 15 156 44 82 60 83 0.3324079 南ことり 16 159 42 80 58 80 0.3366060 園田海未 16 159 48 76 58 80 0.3502635 西木野真姫 15 161 43 78 56 83 0.3523379 東條希 17 156 43 90 60 82 0.3989719 絢瀬絵里 17 162 47 88 60 84 0.4103085
パラメータはそもそも線形関係が見てわかるので、コサイン距離で計算しなくてもよかったかもしれないけど、こんな感じになった。
ちなみに一番BBAになったのは諸星きらりである。
kNNの計算にcccgパッケージを使った。
V1,V2 1.2064275403,8.4123651097 2.0327395939,9.0011381683 2.7535650023,9.5899112269 1.1185220027,7.9805982001 1.7514418735,8.5104939528 2.507429497,8.6675001018 3.4919715183,8.9226350939 2.0854829165,7.2936962984 2.7711461098,7.9020951256 3.1403493678,7.1759416867 3.5798770559,7.7647147453 4.0369858515,6.9796840005 4.4413513245,6.7638005457 5.0566900878,7.058187075 4.7402301524,6.2142790243 5.5489610985,6.4301624791 4.1248913891,5.6451317343 5.0566900878,5.5862544285 6.1115565392,5.6451317343 6.6214086574,6.2339047929 6.3928542596,4.9778556013 7.043355238,4.683469072 7.6938562163,4.6049659975 7.3598151734,3.4862971861 8.1509650119,3.702180641 8.9069526354,4.0554444761 9.2937370009,2.7405179786 8.3443571947,2.6816406727 7.9399917216,2.1909964572 7.0257741304,1.2685853321 6.8499630552,2.1909964572 5.8302588189,2.2891253003 5.935745464,1.4059657124 5.0742711954,1.5433460928 4.019404744,1.6414749359 4.6699057223,1.1508307204 3.8084314537,1.2293337948 2.5250106045,1.739603779 3.7205259161,2.7993952844 4.3534457869,2.8778983589 3.5271337333,3.4862971861 2.4546861744,3.133033351
# trajectory_detection のスクリプトっぽいもの library(cccd) library(proxy) trajectory_detection <- function(dat, s=1, k=30, l=5, ng=20, nl=20, alpha=0.95, niter=20){ traj_ng <- NULL # cosine距離の使用 cos_dist <- as.matrix(dist(dat, method="cosine")) ng0 <- 0 g <- nng(dat, k=k, mutual=T) while(ng0 < ng){ # 隣接行列の作成 mat <- as.matrix(get.adjacency(g)) # l-out-of-k NNG の作成 for(i in seq(nrow(mat))){ idx1 <- which(mat[i, ] == 1) if(length(idx1) >= l){ idx <- sample(idx1, size=l) mat[i, setdiff(idx1, idx)] <- 0 mat[setdiff(idx1, idx), i] <- 0 } } mat[mat == 1] <- cos_dist[mat == 1] g1 <- graph.adjacency(mat, mode="undirected", weighted=TRUE) # l-k-NNG の最短経路 sh_path <- shortest.paths(g1, algorithm="dijkstra") if(all(sh_path[s,] < Inf)){ # 適当に100回くらいやればtrajencyは収束するっぽいので traj <- NULL for(line in seq(niter)){ # waypoints の導入 waypoints <- sort(sample(setdiff(seq(nrow(dat)), s), nl)) # 重み付け距離っぽいもの w <- matrix(0, nr=nrow(dat), nc=nl) for(i in seq(nrow(w))){ for(j in seq(ncol(w))){ w[i, j] <- (sh_path[i, waypoints[j]]^2)/sum(sh_path[i, waypoints]^2) } } idx <- ifelse(outer(sh_path[s,], sh_path[s, waypoints], ">"), 1, -1) tij <- sweep(sh_path[, waypoints]*idx, 2, sh_path[s, waypoints], "+") traj <- cbind(traj, rowSums(tij*w/nl)) } traj_ng <- cbind(traj_ng, traj[, ncol(traj)]) ng0 <- ng0 + 1 print(ng0) } } traj_ng_m <- rowMeans(traj_ng) matplot(traj_ng, pch=16, col=1) points(traj_ng_m, type="l", col=2, lwd=3) traj_ng_q <- quantile(rowMeans(traj_ng), c(1-alpha, alpha)) traj <- (traj_ng_m - traj_ng_q[1])/diff(traj_ng_q) traj <- (traj-min(traj))/diff(range(traj)) return(list(score=traj, s=s, k=k, l=l, nl=nl, ng=ng)) }
# 上のデータを読み込んで dat <- as.matrix(read.csv("wanderlust.txt")) k <- 6 l <- 3 g <- nng(dat, k=k, mutual=T) V(g)$size <- 0 par(mar=c(rep(0, 4))) plot(g) legend("topright", legend=paste("k =", k), bty="n", cex=2) # 適当に l-k-NN graphを作ってみる mat <- as.matrix(get.adjacency(g)) for(i in seq(nrow(mat))){ idx1 <- which(mat[i, ] == 1) if(length(idx1) >= l){ idx <- sample(idx1, size=l) mat[i, setdiff(idx1, idx)] <- 0 mat[setdiff(idx1, idx), i] <- 0 } } sample(which(mat[i, ] > 0), l), 1), seq(nrow(mat))) # cosine距離の使用 cos_dist <- diag(0, nrow(dat)) for(i in seq(nrow(dat))){ for(j in seq(nrow(dat))){ cos_dist[i, j] <- 1-t(dat[i,])%*%t(t(dat[j,]))/sqrt(t(dat[i,])%*%t(t(dat[i,])) * t(dat[j,])%*%t(t(dat[j,]))) } } mat[mat == 1] <- cos_dist[mat == 1] s <- 1 nl <- 6 waypoints <- sort(sample(setdiff(seq(nrow(dat)), s), nl)) g1 <- graph.adjacency(mat, mode="undirected", weighted=TRUE) V(g1)$size <- 0 V(g1)$label.color <- replace(rep(1, nrow(dat)), waypoints, 2) E(g1)$color <- "blue" idx <- mapply(function(x) which(apply(get.data.frame(g), 1, setequal, get.data.frame(g1)[x, 1:2])), seq(nrow(get.edgelist(g1)))) par(mar=c(rep(0, 4))) E(g)$color <- replace(rep(grey(0.8), nrow(get.edgelist(g))), idx, "blue") E(g)$width <- replace(rep(1, nrow(get.edgelist(g))), idx, 2) V(g)$label.color <- replace(rep(1, nrow(dat)), waypoints, 2) V(g)$size <- replace(rep(0, nrow(dat)), c(1, waypoints), 8) V(g)$color <- replace(rep(grey(0.8), nrow(dat)), 1, "red") plot(g) legend("topright", legend=paste("k =", k), bty="n", cex=2)
# dijkstra 法によって最短経路を計算する sh_path <- shortest.paths(g1, algorithm="dijkstra") # 各ノードから waypoints までの距離 matplot(t(sh_path[waypoints,]), type="l", xlab="# Node", ylab="distance", lty=1, lwd=2) abline(v=waypoints, lty=3) legend("topleft", legend=waypoints, lwd=2, col=seq(waypoints), bg="white")
# アイドルデータに対して実行 # 名前と所属グループの column は不要。数値のみにしておく。 dat <- read.delim("girl.txt") ll <- subset(dat, type=="LoveLive") rownames(dat) <- dat$name dat <- as.matrix(subset(dat, select=-c(name, type))) # s の設定 loli <- which(rownames(dat) == "横山千佳") res <- trajectory_detection(dat, s=loli, k=25, l=7, nl=20, ng=10) par(mfrow=c(2, 3)) for(i in seq(ncol(dat))){ plot(res$score, dat[, i], xlab="Trajectory score", ylab=colnames(dat)[i], main=colnames(dat)[i], cex.main=2) }
library(jpeg) # 適当なところから jpg を取ってきておく。 pics <- vector("list", 9) for(i in seq(9)){ pics[[i]] <- readJPEG(paste(i, "_add.jpg", sep="")) } ra <- 1 #原点に近いところが潰れるので拡大したかったけど、等倍でやった。 xy0 <- sapply(pics, dim)[1:2, ] #pixel rownames(xy0) <- c("height", "width") s1 <- 0.002 #拡大縮小率 s0 <- 0.6 #lay0 <- matrix(unlist(locator(9)), nc=2) x_idx <- sort(match(match(ll$name, rownames(dat)), order(res$score))) y_idx <- sort(cbind(dat, res$score)[order(res$score)[x_idx], 7]) par(mar=c(4, 4.5, 4, 2)) plot(sort(res$score), pch=16, xlab="アイドル人数", ylab="ロリ ← Trajectory score → ボイン", cex.lab=1.6, cex.axis=1.6, las=1) points(x_idx, y_idx, pch=16, col=2, cex=1.5) title("Wanderlust: Trajectory score による成長度", cex.main=1.5) for(j in seq(pics)){ arrows(lay0[j, 1], lay0[j, 2], x_idx[j], y_idx[j], lwd=2, length=0.12) xleft=lay0[j, 1]*ra - xy0[2, j]/2*s0 ybottom=lay0[j, 2]*ra - xy0[1, j]/2*s1 xright=lay0[j, 1]*ra + xy0[2, j]/2*s0 ytop=lay0[j, 2]*ra + xy0[1, j]/2*s1 rasterImage(image=pics[[j]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE) }
