Rでフローサイトメトリー(FCM/FACS)

R Rpackage Rを使いこなす 医学

MikuHatsune2013-09-20

フローサイトメトリー(FCM)という実験がある。これは、細胞表面に抗原があるのを抗体反応で標識して、標識された細胞に光を当てて抗体反応が起きていたか、つまりその抗原が存在するかを識別する実験である。
実験機器が許す限り、抗原標識は可能で、たいてい2〜3の抗体、多くて何十できるとかできないとか。
正常な細胞でも白血病などの異常な細胞でも、ある種の抗原パターンがあって、それらの組み合わせで理論的には細胞を分離できる。
flowCoreパッケージでFCMデータを扱ってみる。

# Required once per R installation
source("http://bioconductor.org/biocLite.R")
biocLite("flowCore")

library(flowCore)
file.name <- system.file("extdata","0877408774.B08", package="flowCore")
x <- read.FCS(file.name, transformation=FALSE)

# 定量データの部分だけ別に取っておく。時系列も抜いておく。
dd <- x@exprs[, grep("-", colnames(x@exprs))]
time0 <- x@exprs[, "Time"]

サンプルファイルを読み込んでみると、7つの抗体で標識された50 sec ほどのデータであることが分かる。

x

flowFrame object '0877408774.B08'
with 10000 cells and 8 observables:
     name              desc range minRange maxRange
$P1 FSC-H             FSC-H  1024        0     1023
$P2 SSC-H             SSC-H  1024        0     1023
$P3 FL1-H                    1024        0     1023
$P4 FL2-H                    1024        0     1023
$P5 FL3-H                    1024        0     1023
$P6 FL1-A              <NA>  1024        0     1023
$P7 FL4-H                    1024        0     1023
$P8  Time Time (51.20 sec.)  1024        0     1023
147 keywords are stored in the 'description' slot

実験データは1細胞当たり、7つの抗体それぞれの光量の実数データを持つ。

head(x@exprs)

     FSC-H SSC-H FL1-H FL2-H FL3-H FL1-A FL4-H Time
[1,]   382    77   618     0   225    55   286    1
[2,]   628   280   245   431   259     0   371    1
[3,]  1023   735   699   448   215   143   638    1
[4,]   373   128   202   354    94     0   149    1
[5,]  1023  1023   618   742   408    61   866    1
[6,]   489   292   179   374   154     0   363    1

とりあえずかっこよく散布図を描く。なんとなくクラスターっぽいのがある。

library(GGally)
colnames(dd) <- gsub("-", "_", colnames(dd)) # - が入っているとなぜかエラー
ggpairs(dd,
	upper = list(continuous = "density", combo = "box"),
	lower = list(continuous = "points", combo = "dot"),
	params=c(cex=0.3, col=memb),
)


ある抗原があるかないかは、gating という閾値を決めてその値を超えるものはさらに次のgatingを満たすものは以下略というようにして絞っていく。
ここでmultiple experienceを読み込んでgatingしてみる。

# loading 5 experiment together into an object fs
fs <- read.flowSet(path=system.file("extdata", "compdata", "data", package="flowCore"), name.keyword="SAMPLE ID", phenoData = list(name="SAMPLE ID", Filename="$FIL"))
pData(phenoData(fs)) # meta data for 5 samples
summary(fs)
x1 <- fs[[1]]            # same structure as x above

Gatingを行う。

plot(x1@exprs[, c("FL1-H", "FL2-H")], cex=0.5)
# どこに垂直・水平の２直線を引きたい？ How do you draw vertical and holizontal lines?
abline(v=c(0, 12), h=c(0, 12), lty=2, col=2, lwd=3)
# 枠内に入った割合は? Fractio within the gated rectangle?
rectGate <- rectangleGate(filterId="Fluorescence Region", "FL1-H"=c(0, 12), "FL2-H"=c(0, 12))
result <- filter(x1, rectGate)
summary(result)

Fluorescence Region+: 9811 of 10000 events (98.11%)


98%が右下にgatingされている。この区切り方で本当にいいかと思ったりしたら、クラスタリングしてみる。4つくらいに分かれるかなと思ったら

dist.1.2 <- dist(x1@exprs[, c("FL1-H", "FL2-H")])
cl <- hclust(dist.1.2)
plot(cl)
memb <- cutree(cl, k = 4) # クラスタリングの枝葉切り
plot(x1@exprs[, c("FL1-H", "FL2-H")], col=memb, pch=16, cex=1)



この条件でgatingすると、陽性の98%の細胞クラスターと陰性の2%の細胞クラスターが抽出される。

# extracting a subset with the gate
fs2 <- Subset(fs, rectGate)

fs[[1]]

flowFrame object 'NA'
with 10000 cells and 7 observables:
     name       desc range minRange maxRange
$P1 FSC-H FSC-Height  1024        0     1023
$P2 SSC-H SSC-Height  1024        0     1023
$P3 FL1-H       <NA>  1024        1    10000
$P4 FL2-H       <NA>  1024        1    10000
$P5 FL3-H       <NA>  1024        1    10000
$P6 FL1-A       <NA>  1024        0     1023
$P7 FL4-H       <NA>  1024        1    10000
141 keywords are stored in the 'description' slot

fs2[[1]]

flowFrame object 'NA'
with 9811 cells and 7 observables:
     name       desc range minRange maxRange
$P1 FSC-H FSC-Height  1024        0     1023
$P2 SSC-H SSC-Height  1024        0     1023
$P3 FL1-H       <NA>  1024        1    10000
$P4 FL2-H       <NA>  1024        1    10000
$P5 FL3-H       <NA>  1024        1    10000
$P6 FL1-A       <NA>  1024        0     1023
$P7 FL4-H       <NA>  1024        1    10000
141 keywords are stored in the 'description' slot

x1split <- split(x1, rectGate)
show(x1split)      # split into two subsets with the gate

$`Fluorescence Region+`
flowFrame object 'NA (Fluorescence Region+)'
with 9811 cells and 7 observables:
     name       desc range minRange maxRange
$P1 FSC-H FSC-Height  1024        0     1023
$P2 SSC-H SSC-Height  1024        0     1023
$P3 FL1-H       <NA>  1024        1    10000
$P4 FL2-H       <NA>  1024        1    10000
$P5 FL3-H       <NA>  1024        1    10000
$P6 FL1-A       <NA>  1024        0     1023
$P7 FL4-H       <NA>  1024        1    10000
141 keywords are stored in the 'description' slot

$`Fluorescence Region-`
flowFrame object 'NA (Fluorescence Region-)'
with 189 cells and 7 observables:
     name       desc range minRange maxRange
$P1 FSC-H FSC-Height  1024        0     1023
$P2 SSC-H SSC-Height  1024        0     1023
$P3 FL1-H       <NA>  1024        1    10000
$P4 FL2-H       <NA>  1024        1    10000
$P5 FL3-H       <NA>  1024        1    10000
$P6 FL1-A       <NA>  1024        0     1023
$P7 FL4-H       <NA>  1024        1    10000
141 keywords are stored in the 'description' slot

kmeans法でクラスターするのもいいかも。

k0 <- kmeans(x1@exprs, 4)
memb <- k0$cluster
# ggpairs の色付けがうまくいかない
plot(as.data.frame(x1@exprs), col=memb, cex=0.5)


　
外れ値解析を混ぜてもいいかもしれない。

library(Rlof)
lof0 <- lof(dd, 20)
# 95% tile より外れ値スコアが大きい細胞は外れ値と考える。
plot(as.data.frame(dd), col=(lof0 > quantile(lof0, 0.95))+1, cex=0.2)