Zenn
👋

【シングルセル解析】[Seurat][clusterProfiler] 各細胞集団のDEGを一度にエンリッチメント解析して比較

2023/02/16に公開
R
#
bioinformatics
#
seurat
#
scrnaseq
#
clusterprofiler
tech

Seuratで各クラスターの発現変動遺伝子（DEG）を検出し、clusterProfilerで各クラスターのDEGのpathway enrichment解析を行う流れを紹介する。

使用するパッケージの読み込み
library(Seurat)
library(dplyr)
library(clusterProfiler)
library(ggplot2)
library(org.Hs.eg.db)

デモデータ

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

チュートリアルページ通りにクラスタリングまで進める。

デモデータ前処理
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)

SeuratでDEG解析

deg <- FindAllMarkers(object = pbmc, only.pos = T)


FindAllMarkersの出力

クラスターごとのDEGの数を集計を集計するにはtable関数でcluster列を指定すればよい。

クラスターごとのDEG数
table(deg$cluster)


clusterごとのDEGをベクトルにしてリストにまとめる。

split関数に同じ長さのベクトルを渡すと、2つ目のベクトルの情報を元に集計した結果をリストで受け取れる。

gene列をcluster列の情報を元に集計
gene_list <- split(deg$gene, deg$cluster)


gene_list


遺伝子名から遺伝子IDへ変換する。

clusterProfilerでエンリッチメント解析する際には、遺伝子名ではなくEntrez IDが必要となる。clusterProfilerにはID変換のためのbitr機能が用意されている。

# for分でlistの各要素に実行
for(i in names(gene_list)){
  gene_list[[i]] <- bitr(gene_list[[i]], 
                         fromType = "SYMBOL", 
                         toType = "ENTREZID", 
                         OrgDb = "org.Hs.eg.db")
			 }


bitrで変換後のgene_list



遺伝子名とEntrez IDは1:1対応していないので、ID変換ができない遺伝子も出てくる。
変換が成功した遺伝子は以下のように調べることができる。

残った遺伝子数を調べる。
length(unlist(gene_list))/2

bitrの返り値はSYMBOL列, ENTREZID列の2列のデータフレームなので、unlistでベクトル化した後に÷2している。(bitrの際、drop=FALSEをつけると、ID変換ができなかったものはNAとして出力に残すこともできる。)

Entrez IDが無い遺伝子は機能アノテーションがされていないような遺伝子なので気にすることはない。変換遺伝子数が異常に少ない時は何か間違っているかもしれない。


clusterProfiler用に成型

ENTREZID列のみを取り出したものに更新しておく。

for(i in names(gene_list)){
  gene_list[[i]] <- gene_list[[i]]$ENTREZID
}


gene_list


compareCluster

clusterProfilerにはenrichGOやenrichKEGGのように遺伝子ベクトルに対してエンリッチメント解析を行う機能があるが、compareCluster()を使うと複数の遺伝子ベクトルに対して比較エンリッチメント解析を行うことができる。
scRNAseqではクラスターごとのDEGを求めることが多いので、compareCluster()は有用である。

geneClusters=引数に上記で作成した遺伝子ベクトルのリストを指定するまではどのデータベースを使う場合でも同じであるが、それ以外の引数はfun=に指定した関数 (enrichGOやenrichKEGG)によって変わるので、各関数のhelpを参照してほしい。


Gene Ontology

ont=引数で"MF", "BP", "CC"を選択。デフォルトは"MF" (molecular function)

GOBP <- compareCluster(geneClusters = gene_list, 
                       fun = "enrichGO", 
                       ont="BP",
                       OrgDb = "org.Hs.eg.db")


KEGG

organism=引数で生物種を記載。デフォルトで"hsa"なのでヒトの場合は無くても良い。

kegg <- compareCluster(geneClusters = gene_list,
		       fun = "enrichKEGG",
		       organism = "hsa")


Reactome

Reactomeデータベースを使用する場合は、ReactomePAパッケージが必要となる。

# BiocManager::install("ReactomePA")
library(ReactomePA)

organism=の指定名はkeggと異なり"human"。デフォルトで"human"なのでヒトの場合は指定しなくても良い。

reactome <- compareCluster(geneClusters = gene_list, 
                           fun = "enrichPathway",
			   organism ="human")


WikiPathway

wikipathwayの際はorganism=引数が必要。keggとは異なり"Homo sapience"と入れる。

WP_clusterplot <- compareCluster(geneClusters = gene_list,
                                 fun = "enrichWP",
				 organism = "Homo sapiens")

指定できる生物種はget_wp_organisms()で確認可能。


Disease Ontology

https://disease-ontology.org/

Disiase Ontologyを使用する場合は、DOSEパッケージが必要となる。ヒトの疾患のデータベースなので生物種の指定は無く、ヒト用。

# BiocManager::install("DOSE")
library(DOSE)
DO <- compareCluster(geneClusters = gene_list,
                     fun = "enrichDO")


DisGeNET

https://www.disgenet.org/

DisGeNETを使用する場合も、DOSEパッケージが必要となる。こちらもヒト用。

library(DOSE)
DGN <- compareCluster(geneClusters = gene_list,
                     fun = "enrichDGN")


結果のplot

clusterProfilerでは結果のplot機能も用意されている。試しにdotplot()の例を示す。

dotの大きさ： pathwayのgene setのうち、ヒットした遺伝子の割合。
dotの色：調整済みp値

dotplot(GOBP)


clusterProfilerで用意されているplot機能は内部でggplot2を使用しているため、ggplot2によるplotの編集が可能。特にpathway term名が長い場合は、潰れて見えなくなることもあるので文字サイズや行間の編集が必要になることもしばしば。

dotplot(GOBP)  + 
  theme(axis.text.y = element_text(size = 6, lineheight = 0.6),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 8))


plotコマンドはたくさん用意されているので公式ページの紹介を参照してほしい。


ENTREZIDから遺伝子名へ変換

上記までのcompareCluster()の返り値はENTREZIDを入力に使用しているため、出力もENTREZIDのままである。SetReadable()関数を使用すると、compareCluster()の出力オブジェクト内のENTREZIDを遺伝子名に戻してくれる。

GOBP <- setReadable(x = GOBP, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
kegg <- setReadable(x = kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")



初めからreadble=引数をTRUEにしておくことでも対応可能。

GOBP <- compareCluster(geneClusters = gene_list, fun = "enrichGO", ont="BP", OrgDb = "org.Hs.eg.db", readable = TRUE)
reactome <- compareCluster(geneClusters = gene_list, fun = "enrichPathway", readable = TRUE)

ただし、enrichKEGG()のようにreadble=引数が用意されていないものもある。それらはSetReadable()を追加で行う必要がある。


data.frameへの変換

clusterProfilerの出力は複雑なリストのような形状をしているが、data.frame()関数で変換するだけでtable形式になる。またwrite.csv等で書き出すだけでtable形式のものが出力される。

df <- data.frame(reactome)




エラー対処履歴

clusterProfilerはupdateの度に前まで使えていた機能がエラーをはくようになったりと不安定な印象があります。


Discussion