[R][Seurat] SCTransformで遺伝子数が減るのを防ぐ

SeuratのSCTransformを使用している例が増えてきたが、実際にやってみると、RNA assayに加えてSCT assayが保存されるためデータ容量も増えるし、RNA assayよりもSCT assayでは遺伝子数が減少するのが嫌で、見たい遺伝子がSCT assayには残ってこないことがあった。

遺伝子数が減少しないようにできたのでメモしておく。

事前準備

library(Seurat)
library(sctransform)

Seuratのチュートリアルページから取得できるPBMC 2700細胞のデータを使用。
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

データ読み込み

mat <- Read10X(data.dir = "pbmc3k_filtered_gene_bc_matrices\filtered_gene_bc_matrices\hg19")

この時点では32738遺伝子×2700細胞のmatrix

Seurat objectへの変換

典型例に従い、min.cells = 3, min.features = 200のフィルタリング条件を使用。

sc <- CreateSeuratObject(mat, min.cells = 3, min.features = 200)

この時点で13714遺伝子×2700細胞になった。

min.cells = 3は、その遺伝子を発現している細胞が3細胞以上の遺伝子を残すという条件。
min.features = 200は、1つの細胞の中で検出された遺伝子数が200以上の細胞を残すという条件。

上記の例では、1,2細胞でしか検出されていない遺伝子が2万個近くあったということ。
各細胞はちゃんと200種類以上の遺伝子が検出されていたので全細胞が条件をパスしている。

---

ちなみに各細胞で検出された遺伝子数の分布はこのような感じ。カウント値が0より大きいかどうかbool値に変換して、列方向に集計すると、最も少ない遺伝子数を持つ細胞でも212の遺伝子が検出されている。

遺伝子のフィルタリング。カウント値が0より大きいかどうかbool値に変換して、行方向に集計。

SCTransformの実行

sc <- SCTransform(object = sc)

次元数を確認すると12572遺伝子×2700細胞となり、遺伝子がさらにフィルタリングされてしまっている。

対策

SeuratのSCTransform関数は内部でsctransformパッケージのvst関数を使用している。vst関数のヘルプを見てみるとmin_cells引数のデフォルトが5となっている。

---

vst関数のhelp

---

SCTransformのhelp

---

SeuratのSCTransform関数のヘルプを見てもmin_cellsオプションは無いが、引数の説明にはAdditional parameters passed to sctransform::vstとあり、vst関数のオプションを認識するようである。

# min_cells引数を入れてSCTransform
sc <- SCTransform(object = sc, min_cells = 3)

※ min.cellsではなくmin_cells

これにより遺伝子数が元のRNA assayの遺伝子数のままSCT assayが作られる。

--> つまるところ、RNA assay作成時とSCT assay作成時のfiltering条件が同じになるようにしておけば良い。

SCTransform version.2でも確認

10X GENOMICSが提供しているPBMC 1万細胞のデータを使用
https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_10k_v3?

上記ページのFeature/cell matrix HDF5(filtered)のリンクからh5ファイルをダウンロード。

SeuratのRead10X_h5機能でmatrixを読み込む。

mat <- Read10X_h5(filename = "pbmc_10k_v3_filtered_feature_bc_matrix.h5")

33538遺伝子×11769細胞

Seurat objectへ変換

sc <- CreateSeuratObject(mat, min.cells = 3, min.features = 200)

20292遺伝子×11537細胞。

---

SCTransform (version 2)

sc <- SCTransform(sc, vst.flavor = "v2")

19109遺伝子×11537細胞

---

SCTransform (version 2 & min_cells = 3)

sc <- SCTransform(sc, min_cells = 3, vst.flavor = "v2")

20292遺伝子×11537細胞

--> verssion2でもmin_cells引数は同様に動く。

それでも遺伝子数が減る場合

ncells

vstのsource codeを見ると以下のようにmatrixからn_cells数をサンプリングした後にmin_cells閾値のフィルタリングを設けているので、細胞数が多いデータに対してサンプリングするn_cells数が少ないと、たまたま選ばれた集団でmin_cells条件に満たない遺伝子が生じ得る。

vst sorce code 86-99行目

  if (!is.null(n_cells) && n_cells < ncol(umi)) {
    cells_step1 <- sample(x = colnames(umi), size = n_cells)
    
    {中略}
    
    genes_cell_count_step1 <- rowSums(umi[, cells_step1] > 
      0)
    genes_step1 <- rownames(umi)[genes_cell_count_step1 >= 
      min_cells]

umiは発現マトリクス

SCTransformのsource codeの62行目を見るとSCTransformのncells引数がvstのn_cells引数に渡されるとのこと。

SCTransform sorce code 56-62行目

  vst.args[["umi"]] <- umi
  vst.args[["cell_attr"]] <- cell.attr
  vst.args[["verbosity"]] <- as.numeric(x = verbose) * 2
  vst.args[["return_cell_attr"]] <- TRUE
  vst.args[["return_gene_attr"]] <- TRUE
  vst.args[["return_corrected_umi"]] <- do.correct.umi
  vst.args[["n_cells"]] <- min(ncells, ncol(x = umi))

全細胞を使ってSCTransform。vstのn_cells引数ではなく、SCTransformのncells引数で指定。

sc <- SCTransform(sc, min_cells = 3, ncells = ncol(sc))