😸

【シングルセル解析】シングルセル単位でgene set scoring (1) ~ AddModuleScore / UCell ~

2023/03/27に公開

bioinformatics

seurat

scrnaseq

tech

特定のpathway geneや細胞種markerなどの遺伝子群 (gene set)を使って、gene set scoringを行うことでactivity単位の評価ができる。
これを1細胞単位で行う手法について紹介する。

AddModuleScore
UCell
scGSEA ※ 別記事
scGSVA ※ 別記事
AUCell ※ 別記事
スコアの有意差検定 ※ 別記事

データの用意

Seuratのチュートリアルページから取得できるPBMC 2700細胞のデータを使用。

チュートリアルと同様に処理を進める。

デモデータの前処理

使用するパッケージの読み込み

library(Seurat)
library(dplyr)

データ処理

pbmc.data <- Read10X(data.dir = "pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

pbmc <- pbmc %>% 
  NormalizeData() %>% 
  FindVariableFeatures() %>% 
  ScaleData() %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:10) %>% 
  FindNeighbors(dims = 1:10) %>% 
  FindClusters(resolution = 0.5)
  
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
    "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
pbmc$Cell_type <- pbmc@active.ident
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

AddModuleScore

Seuratに搭載されている機能。
指定したgene setの平均発現量を基にスコアリングする。

features=引数にスコアリングしたいgene setを指定するのだが、遺伝子名のベクトルではなく、遺伝子名のベクトルをlistにしてから指定する点には注意。

pbmc <- AddModuleScore(object = pbmc,
                       features = list(c("CD8A","GZMA")),
                       name = "CD8_Tcell_score")

name引数に指定した文字列の後ろに1が付いた列名でModule scoreがmeta.dataに保存される。

scoreをFeaturePlot

FeaturePlot(pbmc,"CD8_Tcell_score1")

あくまで平均発現量ベースなので、特異性の低い遺伝子を入れるとスコアの差が出にくくなる。
例えば以下のように、CD8+ TcellにCD45やCD3が出ているからと言って追加すると、他のリンパ球にも十分に発現している遺伝子のためスコアのメリハリが無くなる。

特異性の低い遺伝子をgene setに入れた場合

pbmc <- AddModuleScore(object = pbmc,
                       features = list(c("PTPRC","CD3E","CD3D","CD8A","GZMA")),
                       name = "CD8_Tcell_score")
FeaturePlot(pbmc,"CD8_Tcell_score1")

複数のgenesetでスコアリングする場合は、features=に指定するリストの要素を増やす。
name=はlistの長さ分用意する必要は無いが、以下のように対応する文字列ベクトルを用意してもよい。

pbmc <- AddModuleScore(object = pbmc,
                        features = list(c("CD8A","GZMA"),
                                        c("MS4A1","CD19"),
                                        c("CD14","CD68")),
                        name = c("CD8_Tcell_score","Bcell_score","Mono_score"))

FeaturePlot(pbmc, features = c("CD8_Tcell_score1","Bcell_score2","Mono_score3"), ncol = 3)

アルゴリズム

Source codeを見たところ以下のような計算のようである。

遺伝子ごとの平均発現値を取得
(特定のクラスターに限らず、全細胞での平均である。)
平均発現量でnbin=引数の数に遺伝子を区切る。デフォルトでは24段階に遺伝子が区切られる。
(遺伝子が12000個あれば、発現量順に並び替えて、発現量レベル1から発現量レベル24まで500個ずつ遺伝子を含む集合に分けるイメージ)
features=引数に指定した各遺伝子が属するbinからctr数の遺伝子をランダムサンプリングする。
(全体からランダムに遺伝子をサンプリングするのではなく、同程度の発現レンジ内からランダムサンプリングする。)
featuresにCD3EとCD8Aを指定した場合、
まずCD3Eが属する発現帯の中からctr数の遺伝子をランダムサンプリング。
次にCD8Aが属する発現帯の中からctr数の遺伝子をランダムサンプリング。
これをfeaturesの遺伝子の数だけ繰り返す。
最終的にctr遺伝子の重複を無くしてctr遺伝子セットとする。
featuresに指定した遺伝子の細胞ごとの平均発現量を取得。
ctr遺伝子の細胞ごとの平均発現量を取得
features平均発現値からcontrol平均発現値を除してModule scoreとする。

AddModuleScoreを再現

# 乱数の固定
set.seed(1)
# gene setをリストでまとめたものを用意
gene.sets.list <- list(CD8_Tcell_score=c("CD8A","GZMA"),
                       Bcell_score=c("MS4A1","CD19"),
                       Mono_score=c("CD14","CD68"))
		       
# 発現マトリクスを取得
mat <- GetAssayData(object = pbmc)
# 遺伝子ごとの平均値を取得
genes.avg <- Matrix::rowMeans(x = mat)
# 昇順に並び替える
genes.avg <- sort(genes.avg)
# 発現レベルで区切る
nbin <- 24
genes.cut <- ggplot2::cut_number(x = genes.avg + rnorm(n = length(genes.avg))/1e+30,
                                 n = nbin)
names(genes.cut) <- names(genes.avg)

## control geneの選出
# gene setの数だけfor文を回す
ctrl <- 100
set.seed(1)
ctrl.genes.list <- vector(mode = "list", length = length(gene.sets.list))
for (i in 1:length(gene.sets.list)) {
  # gene setを1つ取り出す。
  gene_set <- gene.sets.list[[i]]
  # 1つのgene setの遺伝子の数だけfor文を回す
  for (j in 1:length(gene_set)) {
    # gene set遺伝子jが属する発現帯からctrl数の遺伝子をランダムサンプリング
    ctrl.genes <- names(sample(genes.cut[which(genes.cut==genes.cut[gene_set[j]])],
                               size = ctrl, 
                               replace = FALSE))
    # ctrl.genesリストに追加
    ctrl.genes.list[[i]] <- c(ctrl.genes.list[[i]], ctrl.genes)
  }
}
# control遺伝子の重複を無くす
ctrl.genes.list <- lapply(ctrl.genes.list, unique)

## control遺伝子の平均発現値を取得
# 計算結果の保存先matrixの用意
ctrl.scores <- matrix(data = 1,
                      nrow = length(ctrl.genes.list),
                      ncol = ncol(mat))
# gene setsの数だけfor文
for (i in 1:length(ctrl.genes.list)){
  # control遺伝子を取り出す
  ctrl.genes <- ctrl.genes.list[[i]]
  # control遺伝子の平均発現量を細胞ごとに算出して記録
  ctrl.scores[i, ] <- Matrix::colMeans(mat[ctrl.genes,])
}

## gene set遺伝子の平均発現値を算出
gene.set.scores <- matrix(data = 1, 
                          nrow = length(gene.sets.list), 
                          ncol = ncol(mat))
for (i in 1:length(gene.sets.list)) {
  # 1つのgene setを取り出す
  gene_set <- gene.sets.list[[i]]  
  # gene set遺伝子の平均発現量を細胞ごとに算出して記録
  gene.set.scores[i, ] <- Matrix::colMeans(mat[gene_set, ,drop = FALSE])
}

## gene set遺伝子平均発現値からcontrol遺伝子平均発現値を引く
gene.set.scores.final <- gene.set.scores - ctrl.scores

## 最終調整
rownames(gene.set.scores.final) <- names(gene.sets.list)
colnames(gene.set.scores.final) <- colnames(mat)
gene.set.scores.final <- data.frame(t(gene.set.scores.final))

AddModuleScore機能の出力と同じ結果が得られた。

## Seuratオブジェクトに追加
pbmc <- AddMetaData(pbmc, metadata = gene.set.scores.final)
FeaturePlot(pbmc, features = names(gene.sets.list), ncol = 3)

UCell

https://github.com/carmonalab/UCell

マンホイットニーのU統計量を使って、指定したgene setの発現順位に基づいたスコアリングを行う。negative markerも使用できる点は秀逸である。

使用するパッケージの読み込み

library(Seurat)
library(UCell) # BiocManager::install("UCell")

AddModuleScore_UCell()機能にSeuratオブジェクトを指定可能である。戻り値もSeuratオブジェクト。
features=引数にはgene setのベクトルをlistにしてから渡す。リストの名前属性がスコアの保存名になるので明示しておくとよい。

CD8 T cellマーカーでCD8.T.cellスコアを算出

pbmc <- AddModuleScore_UCell(obj = pbmc, 
			     features = list(CD8.T.cell = c("CD8A","GZMA")))

gene setの名前に_UCellが末尾についた列名でmeta.dataスロットに保存される。

SeuratのAddModuleScore()よりメリハリがついているように見える。

FeaturePlot(pbmc,"CD8.T.cell_UCell")

複数のgene setを使うときはlistに追加すればよい。

pbmc <- AddModuleScore_UCell(obj = pbmc,
                             features = list(CD8.T.cell = c("CD8A","GZMA"),
                                             B.cell = c("MS4A1","CD19"),
                                             Monocyte = c("CD14","CD68"))
                             )
FeaturePlot(pbmc,c("CD8.T.cell_UCell","B.cell_UCell","Monocyte_UCell"), ncol = 3)

アルゴリズム

遺伝子の発現順位を使用するのだが、シングルセルRNAseqデータでは1細胞内の多くの遺伝子が発現値0であり、ランキング下位は同値になってしまう。また各細胞で発現値が0である遺伝子の数は違うため計算バイアスが生じる。
そこでUCellでは上位1500個までを評価の対象にする。(ランク外の遺伝子は順位を1501として計算に用いる)。maxRank=引数でこの値は変更可能である。

1.指定したgene setのU統計量を細胞ごとに算出する。

U値は gene set遺伝子の発現順位の和 - n(n+1)/2で求められる。(n: gene set遺伝子数)
例えばある細胞において、CD8Aが発現順位5位、GZMAが20位なら順位和は25。U値は25 - 2(2+1)/2 = 22である。

ランク外の遺伝子はmaxRank+1の値が順位である。上記例にランク外の遺伝子がもう一つ追加されると、U値は (5 + 20 + 1501) - 3(3+1)/2 = 169 となる。

2.U値の補正

指定した遺伝子数などによってU値が取り得る範囲が変わってしまうため、これを「gene set遺伝子数 * maxRank」の値で割って補正する。さらにその値を1から引いてUCell scoreとする。

上記の例だと、以下のようになる。
U値 22 / (遺伝子数 2 * makRank 1500) = 0.007333..
1 - 0.007333.. = 0.99266..

U値 169 / (遺伝子数 3 * makRank 1500) = 0.037555..　
1 - 0.037555.. = 0.9624..

この計算によりUCell scoreは0-1の範囲を取り、1に近いほど高い値となる。

Ucellを再現

UCellの再現

# 発現マトリクス取得
mat <- Seurat::GetAssayData(object = pbmc, slot = "data")
# 昇順にランク付け
rank_matrix <- apply(-as.matrix(mat), 2, rank, ties.method="average")
# maxRankより大きい箇所はmaxRank+1の値を代入
maxRank <- 1500
rank_matrix[rank_matrix>maxRank] <- maxRank + 1

# UCell scoreの保存先matrixを用意
UCell.scores <- matrix(data = 1,
                       nrow = length(gene.set.list),
                       ncol = ncol(mat))

# gene setごとにfor文
for(i in 1:length(gene.set.list)){
  # 1つのgene setを取り出す
  gene.set <- gene.set.list[[i]]
  # gene set遺伝子だけのランク値matrixを取得
  gene.set.rank_matrix <- rank_matrix[gene.set,]
  # 細胞ごとに繰り返し
  UCell.scores[i,] <- apply(gene.set.rank_matrix, MARGIN = 2, function(x){
    # gene set遺伝子の合計ランク値
    rank_sum <- sum(x)
    # gene set遺伝子の数
    len_sig <- length(x)
    # マンホイットニーのU統計量を算出
    u_value <- rank_sum - (len_sig * (len_sig+1))/2
    # 補正してUCell scoreに
    ucell_score <- 1 - u_value/(len_sig * maxRank)
  })
}

# UCell score matrixの成型
rownames(UCell.scores) <- names(gene.set.list)
colnames(UCell.scores) <- colnames(mat)
UCell.scores <- data.frame(t(UCell.scores))

小数点以下の小さい値で違いはあるが凡そ再現できた。

## Seuratオブジェクトに追加
pbmc <- AddMetaData(pbmc, metadata = UCell.scores)
FeaturePlot(pbmc, features = names(gene.set.list), ncol = 3)

negative marker

UCellではnegative markerも使用することができる。
遺伝子名ベクトルの各遺伝子の末尾に+か-をつけることでpositiveかnegativeかを指定する。

pbmc2 <- AddModuleScore_UCell(obj = pbmc,
	features = list(CD8.T.cell = c("CD8A+","GZMA+",
	"CD4-","FOXP3-","MS4A1-","CD19-","CD14-","CD68-"))
	)

positive gene set、negative gene setそれぞれで上記のUCell scoreを計算し、postive UCell scoreからnegative UCell scoreを引く。

CD8陽性T細胞にはCD4 T細胞のマーカーや単球系のマーカーは出ないと考えられるので、そのランク値は低いはずである。negative markerのランク値が低ければUCell scoreとしては0に近いのでpositive markerのUCell scoreは影響を受けない。
一方で、negative markerがランク上位に入ってくるようであれば、negative marker UCellの値は高くなり、最終的なUCell scoreは低く補正される。

これを使えば、スコアのメリハリがつけやすくなるし、擬陽性のようなスコアを軽減することができる。

この例では、CD8+ T細胞集団でのCD8.T.cellはスコアは下がっていないが、単球系の集団にいたスコアが高い細胞が減っている。

左: negative marker無し　右: negative markerを使用

negative markerを沢山指定するほど、UCellスコアは下がっていく。
UCellスコアは最小が0になるようになっているので、それ以下の差は表現できない。
w_neg=引数でnegative marker UCellの重み係数を指定することも可能。デフォルトは1なのでnegative marker UCellスコアがそのまま使用されている。

データの用意

AddModuleScore

アルゴリズム

UCell

アルゴリズム

1.指定したgene setのU統計量を細胞ごとに算出する。

2.U値の補正

negative marker

Discussion