[R][Seurat] シングルセル単位でgene set scoring (2)
特定のpathway geneや細胞種markerなどの遺伝子群 (gene set)を使って、gene set scoringを行うことでactivity単位の評価ができる。
これを1細胞単位で行う手法について紹介する。
使用したデータの用意に関しては[R][Seurat] シングルセル単位でgene set scoring (1)を参照して欲しい。
ssGSEA
GSEA (gene set enrichment analysis)は群間比較によって得られた発現変動値やp値を使用して遺伝子をランク付けし、特定のgene setがランキング上位に集まっているのかを指標にenrichmentを解析する手法である。
GSEAはまず群間比較を必要とするが、ssGSEA (single sample GSEA)は1検体単位でGSEA解析ができる拡張版である。
これをシングルセルデータに使えば、1細胞単位でenrichment解析が行える。
scRNAseqデータにssGSEAを行う方法は幾つかのRパッケージで実装されている。
escapeパッケージ
SeuratオブジェクトやSingleCellExperimentオブジェクト、発現マトリクスに対してssGSEA解析を行うことができる。加えて、escapeパッケージは先述のUCellを選択することもできる。
ライブラリ読み込み
BiocManager::install("escape")
library(escape)
gene setの用意
1) custom gene setの場合
遺伝子名のベクトルをリストでまとめる。
gene.sets1 <- list(CD8.T.cell = c("CD8A","GZMA"),
B.cell = c("MS4A1","CD19"),
Monocyte = c("CD14","CD68"))
2) escapeパッケージ提供のgene setの場合
細胞のマーカーや細胞周期、炎症、代謝などに関する20個のgene setがescapeパッケージで提供されている。
data("escape.gene.sets", package="escape")
gene.sets2 <- escape.gene.sets
3) MSigDBの場合
Molecular Signatures Database (MSigDB)のgene setをgetGeneSets機能で取得することができる。
getGeneSetsの引数
- species: デフォルトで"Homo sapiens"になっている。
- library: MSigDBでは9つの大きい分類 (collection)でgene setをまとめている。この引数では各collectionの頭文字を指定する。collectionの内容は上記リンクから確認して欲しい。
library = c("H","C5") - subcategory: MSigDB C5にはBP/MF/CCがあるように各collectionにはsubcollectionがある。collectionの中で取得するgene setを制限する場合はここでsubcollection名を指定する。
subcaterogy=c("BP") - gene.sets: 特定のgene set名が分かっているのなら、直接ここに指定しても良い。
gene.sets = c("HALLMARK_TNFA_SIGNALING_VIA_NFKB")
# 必ずlibrary引数は指定する。
gene.sets3 <- getGeneSets(library = "H")
取得結果
Scoringの実行
enrichIt機能を使用する。
enrichItの引数
- obj: 発現マトリクス or Seuratオブジェクト or SingleCellExperimentオブジェクトを指定
- gene.sets: 上述のように用意したgene.setsオブジェクトを指定
- method: "ssGSEA"か"UCell"を指定。デフォルトは"ssGSEA"
- groups: 一度に解析する細胞数。多くすると早いがPC負荷がかかる。
- cores: 並列実行コア数
- min.size: gene setに含まれる最小遺伝子数。デフォルトではmin.size=5になっているので、少ない遺伝子数のgene setを使うときはmin.sizeを下げる必要があるので注意。
- ssGSEA.norm: ssGSEA結果を正規化するかどうか。正規化すると0~1の値になる。
ソースコードを見ると、内部でGSVAパッケージを使っているとのこと。GSVA::gsvaやUCell::ScoreSignatures_UCellのオプションも受け付けるようである。
T cell, Myeloid gene setを使ってスコアリングしてみる。遺伝子が3つずつしかないのでmin.size引数を下げる必要がある。
ES.seurat1 <- enrichIt(obj = pbmc,
gene.sets = gene.sets1,
groups = 3000, cores = 2,
min.size = 1, ssGSEA.norm = TRUE)
出力はデータフレーム。
ssGSEAの結果はAddMetaData機能でSeuratオブジェクトに保存することができる。
pbmc <- AddMetaData(pbmc, metadata = ES.seurat1)
FeaturePlot(pbmc, features = colnames(ES.seurat1), ncol=3)
escapeパッケージでは結果のplotコマンドも用意されているので活用するよとい。vignette
scGSVAパッケージ
シングルセルレベルでGSVA解析を行うパッケージであるが、ssGSEAも選択できる。
devtools::install_github("guokai8/scGSVA")
library(scGSVA)
gene setの用意
1) pathway database
buildAnnot機能でGO, KEGG, Reatomeのgene setが取得できる。
go <- buildAnnot(species="human",keytype="SYMBOL",anntype="GO") # OP引数で"BP","CC","MF"の指定も可能
kegg <- buildAnnot(species="human",keytype="SYMBOL",anntype="KEGG")
reactome <- buildAnnot(species="human",keytype="SYMBOL",anntype="Reactome")
2) custom gene set
Gene ID列とAnnot列が必要。遺伝子名をGeneID列、その遺伝子が属するgene set名をAnnot列に記載したdata.frameを作成する。
2列目は決まりは無いが、Annot列と同じ内容のものを入れておくと良い。(3列でないとエラーが出るため。)
geneset <- data.frame(GeneID= c("CD8A","GZMA","MS4A1","CD19","CD14","CD68"),
PATH=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2)),
Annot=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2))
)
自作したgene setのdata.frame
Scoringの実行
scgsva機能を使用する。obj引数には発現マトリクス、Seuratオブジェクト、SingleCellExperimentオブジェクトが指定可能である。
method="ssgsea"にするとssGSEAが行える。内部でGSVAパッケージによるssGSEA解析が使われている。
res <- scgsva(obj = pbmc,
annot = geneset,
method = "ssgsea",
ssgsea.norm = TRUE
)
scgsvaの出力を使ったplot例もGitHubにある。
元のSeuratオブジェクトに追加するにはAddMetaDataを使用するとよい。
pbmc <- AddMetaData(object = pbmc, metadata = res@gsva)
FeaturePlot(pbmc, features = unique(geneset$Annot), ncol=3)
アルゴリズム
こちらにssGSEA機能を自作するRスクリプトがあった。
GSVAパッケージのssgsea機能と同様のスクリプトとなっている。
例として1サンプル、1 gene setの場合を考える。実際にはこの処理をgene setの数だけ、サンプル(細胞)の数だけ繰り返すだけである。
# Seuratオブジェクトから発現マトリクスを取り出す
set.seed(1)
cell <- sample(rownames(pbmc@meta.data[pbmc$Cell_type == "CD8 T",]),1)
mat <- GetAssayData(object = pbmc_CD8T, slot = "count")
# デモ用に1細胞分のデータを取り出す。
X <- mat[,cell] # 13714遺伝子*1細胞のベクトル
# デモ用のgene set。遺伝子名のベクトル。
gene_set <- gene_set <- c("PTPRC","CD3D","CD8A","GZMA","CCL5")
ちなみにこの細胞ではデモgene setの5つの遺伝子のうち4つが検出されていた。
1. 発現値をランク値に変換 (低い値からランク値がつけられる)
R <- rank(X)
この細胞は発現値0の遺伝子が12863個あり、それらが同位でランク値1となる。そしてその次のランク値が12864になる。ties.method=引数で同位の扱いを指定でき、ties.method=averageで12864/2がランク値になる。ties.method=minにすると12863個目までがランク値1になる。デフォルトはties.method=averageになっている。
ランク値と遺伝子発現値のランキングは以下のような関係になる。
ちなみにデモgene setのランク値は以下のようになっていた。
2. 発現ランキングを降順に並び替え
gene_ranks <- sort(R, decreasing = T)
この例ではCCL5は発現ランキングの上位に位置していた。
3. gene setの遺伝子の箇所、gene set以外の遺伝子の箇所を記録
# gene setの箇所はTRUE、それ以外はFALSEのbool値ベクトル
indicator_pos <- names(gene_ranks) %in% gene_set
# bool値を逆転
indicator_neg <- !indicator_pos
今回は5遺伝子なので5か所ヒットしている。
4. gene set遺伝子のランク値のみ取得
bool値を計算に使うとTRUEが1, FALSEが0としてふるまうので、発現ランクベクトルとgene set遺伝子か否かのbool値ベクトルをかけると、gene setの箇所は元のランク値、それ以外は0のランク値ベクトルが得られる。
その結果を0.25乗しているが、これはGSEAの決まり事のようである。(GSVAではalphaは1)
alpha <- 0.25
rank_pos <- (gene_ranks * indicator_pos) ^ alpha
gene setの遺伝子のランク値以外が0、gene setの遺伝子が0.25乗された値になっている。
5. gene set遺伝子のランク値累積和ベクトル
cumsum関数の例
発現ランキングの4番目と7番目にgene setの遺伝子があれば、以下のように合計が累積したベクトルが得られる。
cumsum関数を使用してgene set発現ランクベクトルの累積和ベクトル求める。さらにgene set遺伝子の合計ランク値で割ることで累積和の最大値を1にする。
step_cdf_pos <- cumsum(rank_pos)/sum(rank_pos)
gene set遺伝子の累積和ベクトルを図示すると以下のように左上に膨らんだ線となる。発現ランキングの上位にgene set遺伝子がある程早く累積ランク値が1に達するのが分かる。
逆に発現ランキングの下位にgene set遺伝子が集まっていると、右下にへこんだ線となるのが想像できる。
GSEA/ssGSEAではこの線の膨らみを数値化することでgene setのスコアリングをしている。
6. gene set遺伝子以外の累積和ベクトル
gene set遺伝子のランク値を0, gene set以外の遺伝子はランク値1として累積和ベクトルを求める。こちらも最大値を1に補正する。
step_cdf_neg <- cumsum(indicator_neg)/sum(indicator_neg)
可視化すると以下のようになる。
gene set遺伝子での処理と同じように処理しても凡そ結果は変わらない。
rank_neg <- (gene_ranks * indicator_neg) ^ alpha
step_cdf_neg <- cumsum(rank_neg)/sum(rank_neg)
7. 累積和ベクトルの差
gene set遺伝子の累積ランクベクトル、gene set遺伝子以外の累積ランクベクトルの各要素を減ずる。
step_cdf_diff <- step_cdf_pos - step_cdf_neg
感覚としては以下の網掛け部分を数値化していく。
8. 遺伝子数で補正
解析に使用された全遺伝子数で補正する。
num_genes <- length(X)
step_cdf_diff <- step_cdf_diff / num_genes
9. ssGSEAスコア化
上記までの累積和ベクトルを合計してenrichment scoreとする。
es <- sum(step_cdf_diff)
その後、複数サンプルでの値の範囲で補正する。(この例では1サンプルなのでそれができない。)
複数サンプル(細胞)、複数gene setの場合
上記は主にベクトルで行ってきたが、これを発現マトリクスに拡張して、サンプルごと、gene setごとにapplyで繰り返し処理をする。
Seuratでは最初に全ての細胞で発現が0のような遺伝子をフィルタリングしているが、特定の細胞をsubsetした後には再フィルタリングしないので全細胞で発現値が0のような遺伝子が残り続けてしまう。
GSEAではinputの遺伝子全てにランク値をつけたり、遺伝子数での補正をしたりと、使う遺伝子数の影響を受ける。シングルセル用にGSEAを実装しているパッケージではinputで全細胞が発現値0の遺伝子はまず除かれるようになっている。
gene_sets <- list(CD8.T.cell = c("CD8A","GZMA"),
B.cell = c("MS4A1","CD19"),
Monocyte = c("CD14","CD68"))
mat <- GetAssayData(object = pbmc, slot = "data")
# 10細胞を抽出
X <- as.matrix(mat[,1:10])
# 10細胞全てで発現が0の遺伝子を除く
X <- X[rowSums(X > 0) != 0, ]
alpha <- 0.25
num_genes <- nrow(X)
# 遺伝子のランク値マトリクスを作成 (列単位でランク値に変換)
R <- matrixStats::colRanks(x = X, preserveShape = TRUE, ties.method = "average", useNames = T)
## 1列ごと（1細胞ごと）に繰り返し
es <- apply(R, 2, function(R_col) {
# ランク値降順に並び替え
gene_ranks <- sort(R_col, decreasing = TRUE)
## gene setごとに繰り返し
es_sample <- sapply(gene_sets, function(gene_set) {
# gene setの遺伝子かどうかのbool値ベクトル作成
indicator_pos <- names(gene_ranks) %in% gene_set
indicator_neg <- !indicator_pos
# gene setの箇所のランク値を取得
rank_pos <- (gene_ranks * indicator_pos) ^ alpha
# gene setの遺伝子の累積ランク値ベクトル (合計を1に補正)
step_cdf_pos <- cumsum(rank_pos)/sum(rank_pos)
# gene set以外の遺伝子の累積ランク値ベクトル (合計を1に補正)
step_cdf_neg <- cumsum(indicator_neg)/sum(indicator_neg)
# 累積ランク値ベクトルの差分を取得
step_cdf_diff <- step_cdf_pos - step_cdf_neg
# 遺伝子数で補正
step_cdf_diff <- step_cdf_diff/num_genes
# 累積ランク値ベクトルを合計
sum(step_cdf_diff)
})
# sapplyの出力がlistなのでunlist
unlist(es_sample)
})
# 値の範囲の幅で補正
es_correct <- es/diff(range(es))
es_correct
GSVAパッケージのssGSEAと同様の結果が得られるのが確認できた。
escapeパッケージでは最後の補正を値の範囲で行わずに、0-1にscalingしている。
normalize <- function (x) {
(x - min(x))/(max(x) - min(x))
}
es_correct2 <- apply(es, 1, normalize)
es_correct2
escapeパッケージのenrichItと同様の結果が得られるのが確認できた。
scGSVA
ssGSEAよりenrichment scoreの差が出やすいらしいGSVA解析をシングルセルレベルで行うこともできる。上記のscGSVAパッケージで行える。
※ 処理時間が長い。一晩かかるつもりで臨んでほしい。
library(scGSVA)
geneset <- data.frame(GeneID= c("CD8A","GZMA","MS4A1","CD19","CD14","CD68"),
PATH=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2)),
Annot=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2))
)
res <- scgsva(obj = pbmc,
annot = geneset,
method = "gsva"
)
値は-1から1の範囲を取り、ssGSEAよりもメリハリがついたスコアリングができているようだ。処理時間をかけた割にssGSEAやUCellよりも格段に有能かというと何とも言えないが。
pbmc <- AddMetaData(object = pbmc, metadata = res@gsva)
FeaturePlot(pbmc, features = unique(geneset$Annot), ncol=3)
scGSVAパッケージのソースコードを見ると、遺伝子のフィルタリングしたマトリクスをただ単純にGSVAパッケージでのGSVA解析に渡しているだけのようだった。
アルゴリズム
ssGSEAは1サンプルごとの遺伝子発現量を使用するが、GSVAは全サンプルを使用して、遺伝子の発現量を発現密度に変換する点が異なる。その後、ランク値へ変換して解析するところの考え方は類似している。
ランク値変換後は、発現上位から遺伝子を見ていき、gene set遺伝子であれば加算、gene set遺伝子以外であれば減算を行う。
スコアはその最大値/最小値を使って計算される。scoreの算出方法は3通り用意されている。
GSVAパッケージは一部の機能がC言語で書かれている。以下は私がRコードに書き直したものだが、理解しきれていない部分もあるので信用しすぎないように。
1) 発現マトリクスを発現密度に変換
遺伝子発現の分布に応じて処理が異なる。マイクロアレイデータでは正規分布、RNAseqデータではポアソン分布に従うことを想定して、GSVA::gsva機能のkcdf引数では"Gaussian"か"Poisson"が選択できる。
以下はポアソン分布を使った発現密度変換を示す。
expr <- GetAssayData(pbmc, slot = "count")
# デモ用に10細胞取り出す
expr <- as.matrix(expr)[,1:10]
# どの細胞にも発現していない遺伝子は除く
expr <- expr[rowSums(expr > 0) != 0, ]
# 遺伝子ごとに繰り返す
gene.density <- apply(expr, 1, function(x){
n_sample <- length(x)
# 原著で定めたパラメーター
bw <- 0.5
# 結果の保存先
left_tail <- matrix(0.0, nrow = n_sample, ncol = 1)
## サンプルごとに繰り返す
# ポアソン分布の期待値がx[i]+0.5の時の、各サンプルの発現値の累積確率を求める。これをサンプルの数だけ加算していく。
for(i in 1:n_sample){
left_tail <- left_tail + ppois(x, x[i]+bw, lower.tail = TRUE, log.p = FALSE)
}
# 合計した累積確率をサンプル数で割る。0~1の範囲へ
left_tail <- left_tail/n_sample
# 値を変換。0に近いほどマイナスに、1に近いほどプラスの値になる。
return(-1.0 * log((1.0-left_tail)/left_tail))
}
)
gene.density <- t(gene.density)
Gene densityの算出について、ある遺伝子Xと細胞(サンプル)が5つの場合のイメージを作ってみた。灰色の塗り潰し箇所が累積確率である。
これを全遺伝子に対して行っている。
単純に対数変換するよりは発現の差が表現できているように見える。
| 発現マトリクスを対数変換 | 発現マトリクスを発現密度に変換 |
|---|---|
 |
 |
RNAseqではGC含有量や遺伝子長によるバイアスも発現量のスケール感に影響を与えているため、各遺伝子を同じスケール間の中でランク付けするために発現密度という変換を行っているとのことである。原著
2) ランク値に変換
# 発現密度を降順に並び替えた際のindexを取得
sort.sgn.idxs <- apply(gene.density, 2, order, decreasing=TRUE)
# ランク値に変換
rank.scores <- apply(sort.sgn.idxs, 2, function(x){
num_genes <- length(x)
tmp <- numeric(length = num_genes)
tmp[x] <- abs(seq(from=num_genes, to=1) - num_genes/2) # 中央値ほど低く、発現上位or下位になるほどランク値が高くなる。
return(tmp)
})
3) enrichment score
gene set遺伝子のエンリッチメントスコアを計算する。
mx_diff, abs_rnkの引数で結果が変わる。GSVA::gsvaのデフォルトはmx_diff=TRUE, abs_rnk=FALSEである。
gene_sets <- list(CD8.T.cell = c("CD8A","GZMA"),
B.cell = c("MS4A1","CD19"),
Monocyte = c("CD14","CD68"))
# 遺伝子名をindexに変換
gset.idx.list <- lapply(gene_sets, function(genes) {which(row.names(expr) %in% genes)})
custom_gsva <- function(rank_matrix,
gene_density,
geneset_idx_list,
tau = 1, # <-- ssGSEAで0.25だった箇所
mx_diff = TRUE,
abs_rnk = FALSE){
n_genes <- nrow(rank_matrix)
n_samples <- ncol(rank_matrix)
# 各列をgene.densityの降順に並び替えた時のindexを取得
sidxs <- apply(gene.density, 2, order, decreasing=TRUE)
# gene setの数だけ繰り返し
res.list <- lapply(geneset_idx_list, function(geneset_idxs){
# gene set遺伝子の数
n_geneset <- length(geneset_idxs)
# gene set遺伝子以外の遺伝子で0から1まで累積する際の1 step辺りの値 (1/n値)を計算
dec <- 1 / (n_genes - n_geneset)
# gene set遺伝子かどうかのmaskベクトル
geneset_mask <- rep(0, n_genes)
geneset_mask[geneset_idxs] <- 1
# 結果の保存先
res <- numeric(length = n_samples)
# サンプルの数だけ繰り返し
for (i in 1:n_samples) {
# サンプルiのランク値を取り出す
rank_value <- rank_matrix[,i]
# サンプルiを発現密度降順に並べた際の遺伝子順を取り出す
x_sort_indxs <- sidxs[,i]
# gene set遺伝子のランク値合計を計算しておく
sum_gset <- sum(rank_value[geneset_idxs] ^ tau)
mx_value <- 0
mx_value_sign <- 0
cum_sum <- 0
mx_pos <- 0
mx_neg <- 0
## 発現上位から遺伝子を取り出していき、gene set遺伝子の場合はランク値を加算、gene set遺伝子以外の場合はランク値を減算していく
for (j in 1:n_genes) {
idx <- x_sort_indxs[j]
# geneset_maskは遺伝子一覧のうちgene set遺伝子は1, gene set遺伝子以外は0でマスクしている
if (geneset_mask[idx] == 1) {
# gene set遺伝子の場合、累積ランク値にランク値を加算
cum_sum <- cum_sum + (rank_value[idx] ^ tau / sum_gset)
} else {
# gene set遺伝子以外の場合、累積ランク値から1/n値を引く
cum_sum <- cum_sum - dec
}
# 最大値、最小値を更新
mx_pos <- max(mx_pos, cum_sum)
mx_neg <- min(mx_neg, cum_sum)
}
## enrichment scoreの計算
if (mx_diff != 0) {
# mx_diffがTRUEなら最大値+最小値
mx_value_sign <- mx_pos + mx_neg
# さらにabs_rnkがTRUEなら最大値-最小値
if (abs_rnk != 0)
mx_value_sign <- mx_pos - mx_neg
} else {
# mx_diffがFALSEなら絶対値が大きい方を採用
mx_value_sign <- ifelse(mx_pos > abs(mx_neg), mx_pos, mx_neg)
}
res[i] <- mx_value_sign
}
return(res)
})
return(t(data.frame(res.list)))
}
custom_gsva(rank_matrix = rank.scores, gene_density = gene.density,
geneset_idx_list = gset.idx.list,
tau = 1, mx_diff = TRUE, abs_rnk = FALSE)
scGSVAパッケージのscgsvaと同様の結果が得られるのが確認できた。
geneset <- data.frame(GeneID= c("CD8A","GZMA","MS4A1","CD19","CD14","CD68"),
PATH=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2)),
Annot=c(rep("CD8.T.cell", 2),
rep("B.cell", 2),
rep("Monocyte", 2))
)
res <- scGSVA::scgsva(obj = expr, annot = geneset, method = "gsva")
GSVAパッケージのgsvaも同じ結果。
GSVA::gsva(expr = expr, gset.idx.list = gene_sets, method = "gsva", kcdf = "Poisson")
Discussion