Zenn
⛳

【シングルセル解析】シングルセル単位でgene set scoring (4) ~ スコアの有意差検定 ~

2023/09/25に公開
R
#
bioinformatics
#
seurat
#
scrnaseq
tech

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


変動パスウェイ解析

Geneset/Pathwayを同時にたくさんスコアリングした場合、どのGeneset/Pathwayが変動しているのかを検定したくなる。

scGSVAパッケージではscgsva()機能の出力を使ったpathwayの有意差検定を提供している。limmaパッケージでの線形モデルを使った findPathway()機能 と rstatixパッケージでWilcox/T test/ANOVAを行うsigPathway()機能が使用可能。
https://github.com/guokai8/scGSVA

デモデータ

デモ用にpbmc2700のデータでKEGG pathway databaseのgenesetをssGSEAアルゴリズムでスコアリングしておく。（226個のpathway）

library(scGSVA) 

# genesetの用意
kegg <- buildAnnot(species="human",keytype="SYMBOL",anntype="KEGG")

# scoring
res <- scgsva(obj = pbmc,
              annot = kegg,
              method = "ssgsea",
              ssgsea.norm = TRUE
)


objスロットにはSeuratオブジェクト、gsvaスロットには各genesetスコアのベクトル（長さは細胞数）、annotスロットには用いたgenesetオブジェクトの情報が保存されている。


res@obj


res@gsva


res@annot

findPathway()

findPathway()機能のref=引数にreferenceとしたい群名を指定すると、ref群 vs その他の群1, ref群 vs その他の群2 ,,,, の比較が行われる。
※ ref=引数を指定しない場合は、group=引数に指定した群情報からアルファベット順で最初の水準が使用される。

検定にはlimmaパッケージを用いた線形モデルへの当てはめが行われる。

例として、CD8 T細胞とその他の8つの細胞種で有意なpathwayを探索してみる。

dep1 <- findPathway(object = res, group = "Cell_type", ref = "CD8 T")

8つの比較パターン × 226個のpathway = 1808行の結果が得られる。

findPathway()の出力


計算スクリプト

内部では次のスクリプトが働いている。
線形モデルの切片項に reference群を設定し、説明変数にその他の群 というモデルデザインにデータを当てはめている。

inputの用意
# 行に細胞バーコード、列にpathway
gsva <- res@gsva
# 行にpathway、列に細胞バーコード
gsva <- t(gsva)

# 各細胞の群情報ベクトル
group <- pbmc$Cell_type

# reference群の設定
ref <- "CD8 T"
線形モデルへの当てはめ
# refがある場合は、reference群が水準の1番目になるように
level <- unique(group)
if (!is.null(ref)) {
  level <- c(ref, setdiff(level, ref))
} else {
  ref <- level[1]
}

# reference群とそれ以外の群の比較表作成
lev <- expand.grid(setdiff(level, ref), ref)
# 比較の表記をわかりやすくしておく。@@は後々削除される。ただの文字区切り用。
levn <- apply(lev, 1, function(x){
  paste0(x[1], "_vs_", x[2], "@@")
  })

# 群情報をfactor型にしておく
group <- factor(as.vector(group), levels = level, ordered = F)

# モデルデザインの作成
design <- model.matrix(~group)
colnames(design) <- c(ref, levn)

# 線形モデルへの当てはめ
fit <- limma::lmFit(object = gsva, design = design)
fit1 <- limma::eBayes(fit)
結果の抽出
# 比較パターンが2つ以上あれば、各比較パターンごとに結果を取り出して、最終的に1つにまとめる。
if (length(level) > 2) {
  # 各比較パターンの結果を取得
  tops <- sapply(levn, function(x){
    # topTable機能はlimmaの出力から統計量を含めた結果を取得する機能
    limma::topTable(fit = fit1,
                    adjust = "BH", 
                    coef = x,
                    number = Inf)},
    simplify = F)
  # それをdata.frameとしてまとめる
  res <- do.call(rbind, tops)
  
  res$term <- sub("^\\.", "", sub(".*@@", "", rownames(res)))　## @@以前を削除してterm列へ
  res$comparision <- sub("@@.*", "", rownames(res)) ## @@以降を削除してcomparison列へ
  rownames(res) <- NULL
} else {
  res <- topTable(fit1, adjust = "BH", coef = 2, number = Inf)
  res$comparision <- sub("@@.*", "", levn)
}


sigPathway()

群間の統計検定を実施する機能。検定はT検定（Welch）、Wilcox検定、ANOVAが使用できる。デフォルトはWilcox検定が使用される。test.use=引数で次の文字から指定する。
T検定： "t.test"/"t"　Wilcox検定："wilcox"/"w"　ANOVA："anova"/"aov"/"a"。

ref=引数でreference群を設定した場合は、ref群 vs その他群1, ref群 vs その他群2 ,,,,の検定になるが、NULLの場合はすべての組み合わせで検定が行われる。デフォルトはNULLとなっている。全組み合わせの場合は計算時間がかかるので注意。
※ ANOVAの場合はどの群とどの群で差があるかではなく、genesetエンリッチメントスコアが群によって異なっているか という検定になる。ref=引数は使用されない。


例として、9つの細胞種の全比較パターンで有意なpathwayを探索してみる。

dep_w <- sigPathway(res, group = "Cell_type")

36の比較パターン（ 対角成分を除いた三角行列：(9群 * 9群 - 9) / 2）× 226個のpathwayで8136行の結果が得らえる。

Wilcox (デフォルト)

dep_t <- sigPathway(res, group = "Cell_type", test.use = "t")


T検定

一方、ANOVAであれば226個のpathway分の結果が得られる。比較群が多いときはまずはANOVAで有意なpathwayを見つけるのもよいだろう。

dep_a <- sigPathway(res, group = "Cell_type", test.use = "a")


ANOVA


計算スクリプト

内部では次のスクリプトが働いている。

library(tidyr)
library(rstatix)
inputの用意
# 行に細胞バーコード、列にgenesetのdata.frame
x <- res@gsva
# group列に群情報を追記する
x$group <- pbmc@meta.data$Cell_type

# reference群の設定
ref <- NULL
検定用にデータ整形
# ロング型へ変換。type列がgeneset, val列がエンリッチメントスコア, group列が群情報
x <- gather(data = x, key = type, value = val, - group)
# genesetでグループ化したtibbleに変換
x <- group_by(x, type)
T検定の場合
# Welchのt検定
res <- t_test(data = x, formula = val ~ group, ref.group = ref)
# p値の多重検定補正
res$p.adj <- p.adjust(res$p, method = "BH")
# 必要な列を抽出
res <- res[, c("type", "group1", "group2", "statistic", 
               "p", "p.adj")]
colnames(res)[1] <- "Path"
Wilcoxの場合
# Wilcox
res <- wilcox_test(data = x, formula = val ~ group, ref.group = ref)
# p値の多重検定補正
res$p.adj <- p.adjust(res$p, method = "BH")  
# 必要な列を抽出
res <- res[, c("type", "group1", "group2", "statistic", "p", "p.adj")]
colnames(res)[1] <- "Path"
ANOVAの場合
res <- anova_test(data = x, formula = val ~ group)
res$p.adj <- p.adjust(res$p, method = method)
res <- res[, c("type", "F", "p", "p.adj")]
colnames(res)[1] <- "Path"

Discussion