julia で、flowcytometry のデータを plot する


julia で、flowcytometry のデータを plot してみましょう。 julia の速さがあれば、flowcytometry の解析も、基本的な解析ができます。

読み込みは、julia native であれば、FileIO.jl + FCSFiles が julia 1.0 でも動きます。 この方法では、julia の dictionary 形式で、データを読み込んでくれます。 dictionary を Array に変換する方法は、 別のブログ に記載してありますので、ご参照下さい。 生の FSC データは log や biexpornential への変換が必要になります。 読み込みから、このあたりの操作は、別の機会に記載させて頂きます。

今回は、簡便性のため、正規化された疑似データを plot することにします。 plot の方法は、業界標準の FlowJo の Dot, Pseudocolor, Contour の3種の plot を表現することにします。

まず、パッケージの読み込みと、正規化された疑似データを作ります。 速度の面等では、Plots ですが、Contour が不安定です。 Contour には、PyPlot を使うことにします。

  using Plots, PyCall , PyPlot
  using KernelDensity
  using Distributions
  using StatsBase
  @pyimport seaborn as sns

  data = vcat(
    rand(Normal(),10000,2),
    rand(Normal(3,0.8),10000,2),
    rand(Normal(-2,2),10000,2).+[0 2],
    rand(TDist(10),1000,2)
    )[rand(1:31000,100000) ,:];

  xdata = data[:,1];
  ydata = data[:,2];

散布図 (Dot) の描画

まずは Dot を描きましょう。 

  Plots.plot(data[:,1],data[:,2], t=:scatter, markersize = 0.1, c="black")

ざっと、こんな具合です。初回は時間がかかります。

Pseudocolor

Plots に関する限りは、恐らく FlowJo の Pseudocolor に相当する histogram2d が最も使いやすそうです。 

  clibrary(:colorcet)
  Plots.histogram2d(data[:,1],data[:,2],nbins=200, c=:rainbow)


さすがに FlowJo は売り物なので、もう少し見栄えする気がします。

デフォルトの色も捨て難いです。


Contour

いわゆる等高線グラフであるところの Contour は、かなり骨が折れます。 FlowJo 流だと、一番外の等高線から外れた点を plot しないといけません。 等高線だけなら、seaborn の kdeplot() で描くことができます。 ここに、一番外の等高線から外れた点を選んで scatter しないといけません。 一番外の等高線は、KernelDensity.jl で、おおよそ計算できます。 

  ik = InterpKDE(kde(data))

  function xy_ex(xdata,ydata, ik)
    xy_pdf = zeros(Float64, length(xdata))
      for i in 1:length(xdata)
        xy_pdf[i] = pdf(ik, xdata[i], ydata[i])
      end
    return xy_pdf
  end

  xy_ex_data = xy_ex(xdata,ydata, ik)
  xy_bool = xy_ex_data.< maximum(xy_ex_data)/30 - 0.0005

下は、描画のための関数です。 

  # サンプリング数を減らして描画する関数
  function counterplot(xdata::Array{Float64,1}, ydata::Array{Float64,1}, xy_bool::BitArray{1}, n_sample::Int64)
    select_rows = rand(1:length(xdata), n_sample)
    sns.kdeplot(xdata[select_rows], ydata[select_rows], n_levels=30 , cbar=true , cmap=:inferno)
    sns.scatterplot(xdata[xy_bool], ydata[xy_bool], s=10, color="black")
  end

  # 全部描く関数
  function counterplot_all(xdata, ydata, xy_bool)
    sns.kdeplot(xdata, ydata, n_levels=30 , cbar=true , cmap=:inferno)
    sns.scatterplot(xdata[xy_bool], ydata[xy_bool], s=10, color="black")
  end

ここまでが下準備です。 10^5 個描画すると、かなり時間がかかります。 

  julia> @time counterplot_all(xdata, ydata, xy_bool)
    32.639432 seconds (190 allocations: 76.484 KiB)
  PyObject

さすがに30秒は、つらいので、1万個くらい選んで描画にします。 

  julia> @time counterplot(xdata, ydata, xy_bool, 10000)
    2.076948 seconds (198 allocations: 311.125 KiB)
  PyObject

この位なら実用域ですね。

せっかくだから Makie.jl で Flowcytometry の疑似データを3次元プロットしてみよう。に3次元疑似データのプロットを載せました。

おしまい。

B! LINE