背景
一細胞ごとの蛍光強度を測定する機器としてフローサイトメーターがあります. これを使って細胞を測定すると測定結果としてFCSファイルが出力され, 当研究室ではCytExpertというソフトウェアを用いて解析を行っています. CytExpertはFCSファイルに特化しているだけあって, 散布図やヒストグラムを簡単に描くことができ,
任意の値の範囲内に含まれるデータのみを取り出すゲーティングも直感的に行うことができます. なので, 普通に使う分には不満はないそんなにないと思います. しかし, 自分の研究の中で, 異なる蛍光タンパク質(GFPとmCherryとか)が細胞ごとにどれくらい相関しているのかを調べる必要があり, これがCytExpertではできませんでした. どうしたものかと色々調べた結果,
Rのパッケージを使えば実現できそうだったので, その時の試行錯誤を書きたいと思います.
不満は無いと言ったけど...
まず, CytExpertはWindows専用でMac OSでは動かない. これは地味に面倒ですね. Boot Camp アシスタントやVirtualBoxなどの仮想化ソフトを使う方法もありますが, Windowsのライセンスのお金はどこから出るんですかね？ そしてグラフがイマイチ. 文字がギザギザしているし, 軸のメモリの数など細く設定できません.
発表用データとして使うにはIllustratorやPowerPointで上から線を引き直したり, 軸ラベルを書き直したりすることが必須なのではないでしょうか？にもかかわらず, グラフの出力はグラフ単位ではおこなえず(サイズをいじればできないことはないですが...), どちらかというとプリントアウトを前提としているため, 画像編集には向いていないという仕様... これはRを使うしかない！！
flowCoreを使ってみる
RのパッケージであるflowCoreはバイオ関係のデータを扱うためのパッケージを包括的にまとめたBioconductorで提供されています. インストールはBioconductorのflowCoreのページのInstllationと書かれたところに載っているコードをRに打ち込めばOKです. 同じページにマニュアルが載っているので詳しく知りたい方はそちらを参照してください（もちろん英語です）.
インストールを終えたら, ひとまずFCSファイルを読み込んでみましょう. FCSファイルの読み込みにはread.FCS()を使います. ()の中にはFCSファイルのパスを指定します.
R
library(flowCore) #flowCoreの呼び出し
sample = read.FCS("XXXXX.fcs") #XXXXXはファイル名
さっそくsampleの中身を見ていきたいと思います. まず素直にsampleを出力してみます.
R
sample
flowFrame object 'BY4741.fcs'
with 10671 cells and 22 observables:
name desc range minRange maxRange
$P1 FSC-H FSC-H 16777215 0 16777214
$P2 FSC-A FSC-A 16777215 0 16777214
$P3 SSC-H SSC-H 16777215 0 16777214
$P4 SSC-A SSC-A 16777215 0 16777214
$P5 FL1-H FITC-H 16777215 0 16777214
$P6 FL1-A FITC-A 16777215 -111 16777214
$P7 FL2-H PerCP-H 16777215 0 16777214
$P8 FL2-A PerCP-A 16777215 -111 16777214
$P9 FL3-H PE-H 16777215 0 16777214
$P10 FL3-A PE-A 16777215 -111 16777214
$P11 FL4-H ECD-H 16777215 0 16777214
$P12 FL4-A ECD-A 16777215 -111 16777214
$P13 FL5-H PC5.5-H 16777215 0 16777214
$P14 FL5-A PC5.5-A 16777215 -111 16777214
$P15 FL6-H PC7-H 16777215 0 16777214
$P16 FL6-A PC7-A 16777215 -111 16777214
$P17 FL7-H DAPI-H 16777215 0 16777214
$P18 FL7-A DAPI-A 16777215 -111 16777214
$P19 FL8-H HoechstRed-H 16777215 0 16777214
$P20 FL8-A HoechstRed-A 16777215 -111 16777214
$P21 SSC-Width SSC-Width 10000 0 9999
$P22 Time Time 900000000 0 899999999
440 keywords are stored in the 'description' slot
例として, BY4741.fcsというファイルを読み込んでいます. 結果を見てみると, 細胞数や蛍光の種類などが出力されています. 次にsummary()を使って基本統計量を見てみます.
R
summary(sample)
FSC-H FSC-A SSC-H SSC-A FL1-H FL1-A
Min. 56218.9 45745.8 182.2 99.8 198.000 -570.30
1st Qu. 244505.6 246539.8 176870.4 178636.2 1315.800 1186.80
Median 277527.7 296975.0 224108.5 240329.6 1677.400 1649.60
Mean 288193.9 320831.8 249568.7 275865.8 1849.468 1902.72
3rd Qu. 318508.1 366266.7 287560.1 321290.0 2142.850 2269.50
Max. 896167.2 2333654.2 2332447.5 9163624.0 11889.000 15146.20
FL2-H FL2-A FL3-H FL3-A FL4-H
Min. 187.0000 -314.4000 129.1000 -544.40002 47.200
1st Qu. 490.0500 359.3000 323.2000 -45.90000 1018.850
Median 621.0000 535.2000 375.9000 58.70000 1398.600
Mean 706.7995 644.0544 385.8352 62.75811 1667.769
3rd Qu. 787.1500 770.2000 437.4000 164.60001 1900.000
Max. 22218.6992 25638.1992 1501.2000 1667.09998 38896.102
FL4-A FL5-H FL5-A FL6-H FL6-A
Min. -1492.000 65.2000 -655.4000 48.7000 -533.79999
1st Qu. 367.300 251.4000 -31.4000 232.2000 -78.20000
Median 971.900 299.1000 85.6000 275.2000 24.00000
Mean 1248.423 317.5798 102.8028 282.7895 29.12037
3rd Qu. 1680.150 352.0000 202.9500 320.7500 126.30000
Max. 45497.398 5861.7002 7058.7998 3844.7000 4545.00000
FL7-H FL7-A FL8-H FL8-A SSC-Width
Min. 298.300 -494.500 22.5000 -1032.4000 578.9091
1st Qu. 787.500 709.200 272.8000 29.0500 1027.5735
Median 988.200 1022.400 325.3000 158.9000 1082.1372
Mean 1089.395 1194.914 338.2478 171.1674 1100.4071
3rd Qu. 1236.350 1413.650 385.5000 290.5000 1154.9559
Max. 11108.000 42786.699 3159.1001 3637.2000 4852.2441
Time
Min. 0.0
1st Qu. 463139.5
Median 904585.0
Mean 943792.8
3rd Qu. 1383440.0
Max. 1825126.0
蛍光ごとに結果が出力されます. 平均値を知りたいだけならここまでで十分ですね(もっともこれくらいはCytExpertでもできますけど). さらに詳しくsampleの中身を見ていきます.
R
str(sample)
Formal class 'flowFrame' [package "flowCore"] with 3 slots
..@ exprs : num [1:10671, 1:22] 296358 291005 284861 216855 294244 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : NULL
.. .. ..$ : Named chr [1:22] "FSC-H" "FSC-A" "SSC-H" "SSC-A" ...
.. .. .. ..- attr(*, "names")= chr [1:22] "$P1N" "$P2N" "$P3N" "$P4N" ...
.. ..- attr(*, "ranges")= num [1:22] 16777214 16777214 16777214 16777214 16777214 ...
..@ 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 Transforamtion" ...
.. .. ..@ data :'data.frame': 22 obs. of 5 variables:
.. .. .. ..$ name : 'AsIs' Named chr [1:22] "FSC-H" "FSC-A" "SSC-H" "SSC-A" ...
.. .. .. .. ..- attr(*, "names")= chr [1:22] "$P1N" "$P2N" "$P3N" "$P4N" ...
.. .. .. ..$ desc : 'AsIs' Named chr [1:22] "FSC-H" "FSC-A" "SSC-H" "SSC-A" ...
.. .. .. .. ..- attr(*, "names")= chr [1:22] "$P1S" "$P2S" "$P3S" "$P4S" ...
.. .. .. ..$ range : num [1:22] 16777215 16777215 16777215 16777215 16777215 ...
.. .. .. ..$ minRange: num [1:22] 0 0 0 0 0 -111 0 -111 0 -111 ...
.. .. .. ..$ maxRange: num [1:22] 16777214 16777214 16777214 16777214 16777214 ...
.. .. ..@ dimLabels : chr [1:2] "rowNames" "columnNames"
.. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
.. .. .. .. ..@ .Data:List of 1
.. .. .. .. .. ..$ : int [1:3] 1 1 0
..@ description:List of 440
以下省略
長々と結果が表示されますが, 最初の方だけ見てみるとsampleがflowFrameというS4クラスの構造を持っていることが読み取れます. sampleはexprs, parameters, descriptionの3つのスロットを持ち, このうちexprsに生の蛍光値のデータが含まれています. スロットへのアクセスは@演算子を使います.
そのまま出力すると細胞数x蛍光の種類分だけデータが出力されるのでhead()で先頭のデータだけを表示させます.
R
head(sample@exprs)
FSC-H FSC-A SSC-H SSC-A FL1-H FL1-A FL2-H FL2-A
[1,] 296357.7 338604.3 302511.8 323423.6 2076.3 2329.7 686.9 622.9
[2,] 291004.6 309358.6 230181.7 245796.2 1106.1 972.6 532.2 240.2
[3,] 284861.2 293847.6 216869.6 219008.9 1655.3 1276.1 497.7 267.2
[4,] 216855.4 228669.5 184996.2 192478.2 1784.9 1628.3 495.5 703.8
[5,] 294243.5 338679.2 233823.5 260693.1 1809.5 1944.5 863.6 478.4
[6,] 382794.6 421627.0 299554.3 323549.8 3040.0 2855.5 903.9 606.6
FL3-H FL3-A FL4-H FL4-A FL5-H FL5-A FL6-H FL6-A FL7-H
[1,] 272.3 -104.7 820.1 -85.5 288.8 261.7 259.4 31.4 1296.8
[2,] 365.6 127.8 1476.9 1685.6 316.4 251.4 161.9 -175.5 1172.9
[3,] 336.7 15.5 693.4 597.3 149.2 -105.5 209.4 -41.9 886.6
[4,] 368.8 -35.9 667.4 275.3 248.1 -49.1 214.0 -120.7 854.2
[5,] 420.8 -91.5 1909.3 1096.9 334.9 87.5 265.4 190.2 1039.1
[6,] 526.4 8.3 1791.4 1697.1 322.1 93.9 382.4 -168.5 1232.9
FL7-A FL8-H FL8-A SSC-Width Time
[1,] 1407.0 522.4 775.6 1094.812 0
[2,] 1487.8 426.5 869.5 1093.490 48
[3,] 1125.5 258.7 168.2 1034.124 99
[4,] 832.4 386.2 293.3 1065.435 327
[5,] 1109.8 371.1 281.8 1141.708 425
[6,] 1390.8 631.5 416.0 1106.034 603
蛍光の種類ごとに細胞一つ一つの蛍光値が行列として出力されました. 一応FSC-Aの平均値をmean()を使って計算してみます.
R
mean(sample@exprs[,"FSC-A"])
[1] 320831.8
先ほどのsummary()の結果と一致していますね.
ゲーティングを視覚的に行う
フローサイトメトリーではゲーティングを行うことによって死細胞を取り除いたり, 一定以上の蛍光値を持つ細胞の割合を調べたりします. flowCoreパッケージではゲーティングを行うための関数も用意されています. ゲーティングの形状によっていくつかの関数が用意されていますが, 主に使うのはrectangleGate()とpolygonGate()でしょうか.
とりあえずFSC-AとSSC-Aを軸にとってplotしてみます.
R
plot(sample@exprs[, c("FSC-A", "SSC-A")], xlim = c(0, 1e+06), ylim = c(0, 2e+06))
ここで, 以下に示す図において赤線で囲まれている長方形の部分のデータをゲーティングしてみます. ゲーティングにはrectangleGate()を使います.
R
abline(v=c(1e+05, 6e+05), h=c(0, 7e+05), col="red") #赤線を描画
R
rectGate = rectangleGate(filterId="Rect", "FSC-A"=c(1e+05, 6e+05), "SSC-A"=c(0, 7e+05)) #ゲーティングを作成してRectという名前をつける result = flowCore::filter(sample, rectGate) #flowCoreパッケージのfilter関数であることを明示的に記述 summary(result)
Rect+: 10197 of 10671 events (95.56%)
流れとしては, ゲーティングの枠を作って(1行目), それをfilter()でsampleに適応させる(2行目)イメージです. 結果を見てみると, "Rect"というゲーティングには全10671細胞のうち95.56%にあたる 10197細胞が含まれていることがわかります. rectangleGate()は長方形(2次元の場合. 多次元での指定も可能)のゲーティングですが,
もう少し複雑な領域の指定を行いたい場合はpolygonGate()を用います. 例えば, 下図の赤線で囲まれた領域をゲーティングすることを考えます.
R
loci = data.frame(x = c(122599, 180626.3, 435623.5, 576054, 514038.5, 324890.5, 135691.1), y = c(64095.739, 3031.304, 123772.963, 490907.576,
778359.222, 535815.096, 172651.711))
plot(sample@exprs[, c("FSC-A", "SSC-A")], xlim = c(0, 1e+06), ylim = c(0, 2e+06))
points(loci, col = "red", pch = 20)
lines(loci, col = "red")
R
colnames(loci) = c("FSC-A", "SSC-A") #x軸とy軸に名前をつける
polyGate = polygonGate(filterId = "Poly", loci) #ゲーティングを作成してPolyという名前をつける
result = flowCore::filter(sample, polyGate)
summary(result)
Poly+: 10090 of 10671 events (94.56%)
polygonGate()は座標の行列あるいはデータフレームを受け取り, その各座標で囲まれた部分のゲーティングとなります. この座標をグラフの見た目から決めるのは困難であるため, 視覚的に行えるように工夫します. 具体的には, グラフ上でクリックした位置の座標を返すlocator()を使います.
R
get.loci = function() {
loci = matrix(NA, 1, 2)
for (i in 1:20) {
if (i == 1) { #一番最初の座標の場合
loci = as.data.frame(locator(1))
points(loci[i, ])
} else {
locus = as.data.frame(locator(1))
if (sqrt(((locus - loci[i - 1, ])[1, 1])^2 + ((locus - loci[i - 1, ])[1, 2])^2) < 10) {
lines(x = c(loci[i - 1, 1], loci[1, 1]), y = c(loci[i - 1, 2], loci[1, 2]))
break #前回の座標と今回の座標の距離が10未満だった場合は座標の取得を終了する
}
loci = rbind(loci, locus) #lociに座標を追加していく
points(loci[i, ])
lines(x = c(loci[i - 1, 1], loci[i, 1]), y = c(loci[i - 1, 2], loci[i, 2]))
}
}
return(loci)
}
get.loci関数を定義したので, さっそく使ってみましょう. get.loci()は, あらかじめplotされた図に対してクリックしたところの座標を取得しつつ, 前回取得した座標と今回取得した座標を結ぶ線分を描画していきます. 座標取得を終了するには, 同じ箇所(厳密には非常に近い位置)で2回クリックします.
R
plot(sample@exprs[, c(2, 4)], xlim = c(0, 1e+06), ylim = c(0, 2e+06)) loci = get.loci() loci
x y
1 104843.52 -4040.746
2 294512.18 -27324.112
3 535055.54 55704.479
4 682037.74 319369.555
5 828802.49 688952.186
6 798527.15 1152457.084
7 627169.04 1172626.028
8 459621.81 1005263.237
9 331184.77 761971.105
10 187744.52 471649.839
11 96186.06 315425.527
あとはcolnames()で名前をつけてploygonGate()に放り込めばゲーティングが作成できます. filter()は与えられたゲーティングに含まれる細胞数と割合を返してくれますが, 実際にデータが抽出されているわけではありません. データの抽出にはSubset()を使います. 使い方はfilter()と同じです. 返り値はsampleと同じflowFrameの形です.
R
colnames(loci) = c("FSC-A", "SSC-A")
polyGate = polygonGate(filterId = "Poly", loci)
result = Subset(sample, polyGate)
plot(result@exprs[, c("FSC-A", "SSC-A")], xlim = c(0, 1e+06), ylim = c(0, 2e+06))
Subset()で得られたresultをplotしてみると, ゲーティングした領域のみになっていることがわかります.
ここで得られたresultをさらに別のゲーティングに持ち込むことで複数回のゲーティングを行うことも可能です. 例としてresultをDAPIの蛍光強度によってさらにゲーティングしてみたいと思います.
R
plot(result@exprs[, c("FSC-A", "FL7-A")], xlim = c(0, 7e+05), ylim = c(-500, 8000)) #"FL7-A"はDAPI-Aに相当
loci2 = get.loci()
colnames(loci2) = c("FSC-A", "FL7-A")
polyGate2 = polygonGate(filterId = "Poly2", loci2)
result2 = Subset(result, polyGate2)
plot(result2@exprs[, c("FSC-A", "FL7-A")], xlim = c(0, 7e+05), ylim = c(-500, 8000))
コメントをお書きください