研究資料 Research Materials

参考

• 原著論文
• ソフトウェア

R / R studioの導入

以下を参考にした。（Windowsの場合）

• R言語インストール(2020年 Windows)
• RStudioインストール(2020年 Windows)

Splatterのインストール

1. 以下をR studio のconsoleで実行する。（BiocManagerのインストール）
   install.packages(“BiocManager”)
   BiocManager::install(c(“Biostrings”, “GenomicRanges”, “rtracklayer”))
2. 以下をR studio のconsoleで実行する。（splatterのインストール）
   BiocManager::install(“splatter”)

テスト

1. File → New File → R ScriptでR Scriptを新規作成する。
2. 以下をR scriptに記載する。# Load package
   suppressPackageStartupMessages({
   library(splatter)
   library(scater)
   })
   # Create mock data
   set.seed(1)
   sce <- mockSCE()
   # Estimate parameters from mock data
   params <- splatEstimate(sce)
   # Simulate data using estimated parameters
   sim <- splatSimulate(params)
3. Scriptを実行（Ctrl + Shift + Enter）する。
   sceが参照データ、paramsがsceから推定されたsplatetrのパラメータ集合、simがシミュレーションデータとなる。
   カウントマトリックスはcounts(sim)で参照することができる。

疑似scRNA-seqデータの構成手順

原著論文のFig. 1（下記引用）も参照。
コードフォントはSplatterのパラメータに対応している。

1. 各遺伝子の平均発現量の計算 [Original mean λ’_{i}]
   Shapeパラメータα (mean.shape）とRateパラメータβ（mean.rate）を用いて、ガンマ分布Gamma(α,β)に従う各遺伝子の平均発現量を計算する。
   rowData(sim)$GeneMeanに保存されている。
   バッチ数が1の場合、分散が0のデータ（各遺伝子で値が一定）となっている。
2. 外れ値を含んだ平均発現量の計算 [Gene mean λ_{i}]
   パラメータπ^O=out.probの割合の遺伝子を外れ値として扱い、外れ値の遺伝子をLocatrionパラメータμ^O=out.facLocとScaleパラメータσ^O=out.facScaleに基づく対数正規分布lnN(μ^O σ^O)で計算する。
   外れ値遺伝子iの平均発現量はψ’_{i}Med(λ’_{i})となる。Medはメジアン（中央値）。
   rowData(sim)$GeneMeanに保存されている。
   ここまでは細胞間の不均一性は存在しない（バッチやグループ（細胞種）の違いは除く）。
3. 総RNA数不均一性を含んだ発現量の計算 [Cell mean λ’_{i,j}]
   細胞jの総RNA数（総cDNAライブラリ数）L_jをLocatrionパラメータμ^L=lib.facLocとScaleパラメータσ^L=lib.facScaleに基づく対数正規分布lnN(μ^L, σ^L)で計算する。
   assays(sim)$BaseCellMeanに保存されている。
4. 生物学的不均一性を含んだ発現量の計算 [Trended cell mean λ_{i,j}]
   生物学的変動係数（BCV）B_{i,j}をCommon dispersionパラメータΦ=bcv.commonと自由度df_0=bcv.dfの逆カイ2乗分布に基づいて計算する（詳細はFig.1のBCVの式を参照）。
   その後、B_{i,j}を用いたガンマ分布Gamma(1/B_{i,j}^2,λ’_{i,j}B_{i,j}^2)でλ_{i,j}を取得する。ここで、ガンマ分布のパラメータはshape(k)=1/B_{i,j}^2, scale(θ)=λ’_{i,j}B_{i,j}^2となっており、「1. 各遺伝子の平均発現量の計算」のGamma(α,β)とは対応関係が異なることに注意（α=k, β=1/θ）。
   Gamma(1/B_{i,j}^2,λ’_{i,j}B_{i,j}^2)の平均はλ’_{i,j}、分散はλ’_{i,j}^2B_{i,j}^2である。
   B_{i,j}はassays(sim)$BCVに保存されている。
   λ_{i,j}はassays(sim)$CellMeansに保存されている。
5. カウントデータの計算 [True count Y’_{i,j}]
   パラメータ（平均）λ_{i,j}のポアソン分布でカウントY’_{i,j}を計算する。
   assays(sim)$TrueCountsに保存されている。
6. ドロップアウトを含んだカウントデータの計算 [ Count Y_{i,j} ]
   確率中心パラメータdropout.midおよび形状パラメータdropout.shapeに基づくロジスティック関数でドロップアウト確率π^D_{i,j}を計算し、確率(1-π^D_{i,j})のベルヌーイ分布で1^D_{i,j}（0か1）を計算し、ドロップアウトを含んだカウントデータY_{i,j}を計算する。（図のBer(π^D_{i,j})は誤りであり、Ber(1-π^D_{i,j})が正しいので注意）

パラメータの設定

参考：Splat simulation parameters
パラメータの初期化はparams <- newSplatParams()で行う。このとき、paramにはデフォルトパラメータが入る。
パラメータを更新するときは
params <- setParam(params, “batchCells”, 2000)
のように行う。調整パラメータの意味と設定方法は以下を参照。英名は原著論文のFig. 1（上記にも記載）に対応付けている。

1. batchCells
   細胞数とバッチ数（疑似実験の数）のパラメータ。nBatchesやnCellsは直接編集できないため、細胞数やバッチ数はここで設定する。
   1実験2000細胞 → params <- setParam(params, “batchCells”, 2000)
   2実験1000細胞/500細胞 → params <- setParam(params, “batchCells”, c(1000,500))
2. nGenes
   遺伝子数。
3. group.prob
   細胞のグループ（細胞種）の数とそれぞれの存在確率。２グループ以上はparams <- setParam(params, “group.prob”,c(0.1,0.2,….))のように設定する（c()の括弧内の合計が1になるようにする）。
   1グループ → params <- setParam(params, “group.prob”,1)
   2グループ等確率 → params <- setParam(params, “group.prob”,c(0.5,0.5))
4. mean.shape/mean.rate
   細胞内の平均遺伝子発現量（Original mean）が従うガンマ分布 Gamma(α, β) のパラメータ（α=mean.shape, β=mean.rate）。Wikipediaのガンマ分布の記事の場合、k=mean.shape, 1/θ=mean.rateが対応している。
5. lib.loc/lib.scale/lib.norm
   総RNA数（総cDNAライブラリ数）の細胞間不均一性の設定パラメータ。シミュレーションデータのRead数やUMI数の分布に関連する。
   lib.loc/lib.scale → 対数正規分布 lnN(μ^L, σ^L) のパラメータ（μ^L=lib.loc, σ^L=lib.scale）。
   lib.norm → lib.norm=TRUEとすると上の対数正規分布 lnN(μ^L, σ^L)が正規分布 N(μ^L, σ^L)となる。lib.normは基本的にデフォルト（FALSE）で良さそう。
6. out.prob/out.facLoc/out.facScale
   外れ値（Outlier）の設定。
   out.prob → 外れ値の存在確率。
   out.facLoc/out.facScale → 外れ値の細胞が従う対数正規分布lnN(μ^O, σ^O)のパラメータ（μ^O=out.facLoc, σ^O=out.facScale）。
7. de.prob/de.facLoc/de.facScale
   発現変動遺伝子（Differentially Expressed Genes, DEG）の設定パラメ―タ。
   de.prob → 遺伝子数に対する発現変動遺伝子数の割合。例えばnGenes=20,000、de.prob=0.1であれば、各グループ毎に約2000遺伝子がで発現変動遺伝子となる。ただし、異なるグループで共通の発現変動遺伝子になる場合もある。発現変動遺伝子の確認はrowData(sim)[“DEFacGroup1”]が1でない遺伝子がGroup1の発現変動遺伝子である。de.probが高いほどGroup間の違いが局在化する。
   de.downProb → 発現変動遺伝子が規定の発現値より低くなる確率。「発現変動遺伝子＝既定値より高い発現値を持つ」としたい場合はde.downProb=0とすれば良い。
   de.facLoc/de.facScale → 発現変動遺伝子が従う対数正規分布lnN(μ^D, σ^D)のパラメータ（μ^D=de.facLoc, σ^D=de.facScale）。
8. bcv.common/bcv.df
   生物変動係数（Biological Coefficient of Variance, BCV）の設定パラメータ。
9. dropout.type/dropout.mid/dropout.shape
   ドロップアウトの設定パラメータ。ここで、ドロップアウトとはシークエンス後のカウントマトリックスの非0値が強制的に0になることを表している（これが実際のドロップアウトの現象を表しているか不明）。
   dropout.type → “none”：ドロップアウト無し（デフォルト）。”experiment”：ドロップアウト有り。
   dropout.mid/dropout.shape → ドロップアウトされる確率を決めるパラメータ。dropout.mid=0, dropout.shape < 0とすると発現量に応じてドロップアウトの確率が小さくなる。

※調整できないパラメータ

1. nBatches
   バッチ数（疑似実験の数）。batchCellsで調整する。
2. nCells
   細胞数。batchCellsで調整する。

その他

1. 真の細胞の疑似データはassays(sim)$CellMeansが対応している。
2. RがMicrosoft R Openの状態だとBiocManagerのインストールでうまくいかなかった。