SPADEを使いやすくする

SPADEを昔やったけど、免疫システムの多次元データ解析という総説があって、「FCSの高次元データ、とりあえずSPADE使ってみなよ」的なことが書いてあったのだが、BioconductorのSPADEはFCSデータ専用感が満載なので、行列データが生で存在していても使えるようにした。というのも、FCSデータに限らず、N*Mのサイズの定量値が収まっている行列データであれば、SPADEを使って分化ネットワークっぽいものが作成できるからである。
というわけで、二次元キャラのプロポーションがどう変化していくかやってみる。
データを読み込んで行列にしておく。データはすべて連続値としておく。

library(spade)
library(flowViz)
dat <- read.delim("clipboard")
dat1 <- subset(dat, select=-c(name, type))

関数を探すと、FCSフォーマットでデータを保存する関数があるようなので、使ってみると

write.FCS(dat1, "test.fcs")

 以下にエラー write.FCS(dat1, "test.fcs") : 
  Argument 'x' has to be a 'flowFrame'.

と言われ、flowFrameというオブジェクトでなければならないらしい。flowFrameを調べると

 Objects can be created using
     ‘ new("flowFrame",’
     ‘ exprs = ...., Object of class matrix’
     ‘ parameters = ...., Object of class AnnotatedDataFrame’
     ‘ description = ...., Object of class list’
     ‘ )’
     or the constructor ‘flowFrame’, with mandatory arguments
     ‘exprs’ and optional arguments ‘parameters’ and
     ‘description’.

     ‘flowFrame(exprs, parameters, description=list())’

‘exprs’: Object of class ‘matrix’ containing the measured
          intensities. Rows correspond to cells, columns to the
          different measurement channels. The ‘colnames’ attribute of
          the matrix is supposed to hold the names or identifiers for
          the channels. The ‘rownames’ attribute would usually not be
          set.

     ‘parameters’: An ‘AnnotatedDataFrame’ containing information
          about each column of the ‘flowFrame’. This will generally
          be filled in by ‘read.FCS’ or similar functions using data
          from the FCS keywords describing the parameters.

     ‘description’: A list containing the meta data included in the
          FCS file.

らしい。というわけで実行すると

fcs <- new("flowFrame", as.matrix(dat1))
write.FCS(fcs, "animedata.fcs")

となり、成功する。確認すると、S4クラスのオブジェクトで、@exprs, @parameter, @description という項目ができる。exprs で入力した行列の列名が勝手にparameter で設定される。
GvHD データと比較してみてほしい。

str(fcs)

Formal class 'flowFrame' [package "flowCore"] with 3 slots
  ..@ exprs      : int [1:309, 1:6] 17 18 17 16 16 17 17 17 18 18 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:309] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:6] "age" "height" "weight" "B" ...
  ..@ parameters :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
  .. .. ..@ varMetadata      :'data.frame':	5 obs. of  1 variable:
  .. .. .. ..$ labelDescription: chr [1:5] "Name of Parameter" "Description of Parameter" "Range of Parameter" "Minimum Parameter Value after Transformation" ...
  .. .. ..@ data             :'data.frame':	6 obs. of  5 variables:
  .. .. .. ..$ name    : Factor w/ 6 levels "B","H","W","age",..: 4 5 6 1 3 2
  .. .. .. ..$ desc    : Factor w/ 6 levels "B","H","W","age",..: 4 5 6 1 3 2
  .. .. .. ..$ range   : int [1:6] 31 182 60 105 65 93
  .. .. .. ..$ minRange: int [1:6] 9 127 28 60 47 65
  .. .. .. ..$ maxRange: int [1:6] 31 182 60 105 65 93
  .. .. ..@ dimLabels        : chr [1:2] "rowNames" "columnNames"
  .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slots
  .. .. .. .. ..@ .Data:List of 1
  .. .. .. .. .. ..$ : int [1:3] 1 1 0
  ..@ description:List of 1
  .. ..$ note: chr "empty"

というわけで二次元キャラの変化を追っていくわけだが、左下のほうは年齢も背もおっぱいも小さいのでロリ、右上になっていくと成長していく感じ。
途中の枝でロリグラマラスとか、貧乳とか、そういう分化になるっぽい。

# SPADE
library(spade)
library(flowViz)
dat <- read.delim("clipboard")
dat1 <- subset(dat, select=-c(name, type))
fcs <- new("flowFrame", as.matrix(dat1))
path <- "hoge"       # 保存パス
write.FCS(fcs, path) # FCS データの出力

data_file_path <- path # 上の fcs データを保存したパス
output_dir <- path     # クラスタリング結果やグラフの保存パス
SPADE.driver(data_file_path, out_dir=output_dir)

mst_graph <- igraph0:::read.graph(paste(output_dir,"mst.gml",sep=.Platform$file.sep),format="gml")
clust <- read.table(paste(output_dir, "/clusters.table", sep=""), sep=" ", header=TRUE)
lay0 <- read.table(paste(output_dir,"layout.table", sep=.Platform$file.sep))
library(gplots)
par(mfrow=c(2, 3), mar=c(0, 0, 3, 0))
for(m in head(seq(ncol(clust)), -1)){
	cut0 <- seq(min(clust[,m]), max(clust[,m]), length=99) # ノードの色付け
	g0 <- mst_graph
	V(g0)$label <- NA
	V(g0)$size <- clust[,m]+2
	V(g0)$frame.color <- "black"
	V(g0)$color <- bluered(100)[cut(clust[,m], cut0, include.lowest=TRUE)]
	plot(g0, layout=as.matrix(lay0))
	title(colnames(clust)[m], cex.main=3)
}

各データがどのクラスターに属するかを調べると、SPADEは計算量の削減にdown samplingをして、その後up samplingをして元データを復元、クラスタリングしている。

data_file_path = "animedata.fcs"     # 情報を付け足したい元データ
cells_file_path <- "clusters.fcs"    # クラスターの情報。上のスクリプトを実行していれば自動的に出る
upsample_file_path <- "upsamle.fcs"  # upsampling してクラスター情報を付加した fcs ファイルの出力名
SPADE.addClusterToFCS(data_file_path, upsample_file_path, cells_file_path)
up <- read.FCS("upsamle.fcs")
up@exprs[, "cluster"]                # クラスターの番号が出ている