Zenn
👋

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

2023/02/16に公開・約7,900字
R
#
seurat
#
bioinfomatics
#
scrnaseq
#
clusterprofiler
tech

230216

使用するパッケージの読み込み
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 = sc, only.pos = T)


FindAllMarkersの出力

クラスターごとのDEGの数を集計

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

gene_list <- split(deg$gene, deg$cluster)


gene_list

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

# clusterprofilerのbitr機能でID変換。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



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

length(unlist(gene_list))/2

clusterProfiler用に成型

空のリスト作成して、そこにENTREZIDを入れていく。

genelist <- vector(mode = "list", length = length(gene_list))
names(genelist) <- names(gene_list)
for(i in names(gene_list)){
  genelist[[i]] <- gene_list[[i]]$ENTREZID
}


genelist

compareCluster

複数の群のDEGがある時に、fun引数に指定した機能を同時に実行して、enrichment結果を全群同時に可視化できる。geneClusters引数、fun引数以外の引数はfun=に指定した関数によって変わるので、各関数のhelpを参照すべし。

Gene Ontology

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

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

結果をplot

dotの大きさでそのpathwayに登録されている遺伝子セット (geneset)のうち、ヒットした遺伝子の割合。dotの色で調整済みp値が表現されている。

dotplot(GOBP)




clusterProfilerで用意されている描画用の関数は内部で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))

KEGG

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

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

Reactome

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

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

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

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

WikiPathway

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

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

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

ENTREZIDから遺伝子名へ変換

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

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



初めからreadble引数をTRUEにしておくことでも対応可能。(enrichKeggにはreadable引数が使えず。。。)

reactome <- compareCluster(geneClusters = genelist, fun = "enrichPathway", readable = TRUE)

data.frameへの変換

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

df <- data.frame(reactome)

エラー対処履歴

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

enrichKEGGが使えなくなっていた。

230314
enrichKEGGを行うと以下のメッセージが出るようになった。

--> No gene can be mapped....
--> Expected input gene ID:
--> return NULL...
NULL

clusterProfilerのGitHub issuesを見るとBioconductorのversionを最新にすると直るとか。
私が使っていたversionは3.15で、最新版は3.16だったのでそちらにアップグレード。

BiocManagerのversionを確認
BiocManager::version()

[1] ‘3.15’

BiocManagerのversion 3.16をインストール
BiocManager::install(version = "3.16")

既にインストール済みのBioconductor管理パッケージのアップデートが求められるのでyを入力。
すんなり全てがエラー無くアップデートされれば良いがおそらくうまくいかないパッケージも出てくる。
今回はBiocManagerをアップデート後clusterProfilerを読み込むと以下のエラーが出た。

library(clusterProfiler)
Error: package or namespace load failed for ‘clusterProfiler’ in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]):
名前空間 ‘GOSemSim’ 2.22.0 はロードされていますが、>= 2.23.1 が要求されています

ということでGOSemSimを個別にアップデートしたら、clusterProfilerの読み込みおよびenrichKEGGに成功した。

BiocManager::install("GOSemSim")

そういえば3.15が配布されたのに3.14を使っていた時も同じようなエラーに出くわしたな。。。

enrichWPがタイムアウトする

230317
enrichWPを行うとタイムアウトでうまくいかなかったり、すんなりうまくいったりとよくわからない。

wp <- compareCluster(geneClusters = genelist, fun = "enrichWP", organism = "Homo sapiens")

Warning: URL 'https://data.wikipathways.org/current/gmt/': Timeout of 300 seconds was reachedError in file(con, "r") :
cannot open the connection to 'https://data.wikipathways.org/current/gmt/'

ソースコードを見ると https://data.wikipathways.org/current/gmt/ からwikipathwayの指定した生物種のgmtファイルを取得、成型してenrichr機能に渡しているようである。
enrichWP単独で使うときはタイムアウトしないので、compareCluster内で複数クラスターを扱うときはクラスターごとに毎回gmtファイルのダウンロードがされていてタイムアウトに繋がるのかな??

そこでgmtファイルを1度だけ取得してenrichr機能をcompareClusterのfun引数に渡すようにするとうまくいくようになった。

prepare_WP_data <- function(organism) {
  # wikipathwayのrepogitoryのURL
  wpurl <- 'https://wikipathways-data.wmcloud.org/current/gmt/'
  
  # URLから文字列を読み込む
  x <- readLines(wpurl)
  # .gmtと記載のある要素だけにする
  y <- x[grep('\\.gmt',x)]
  # wikipathways-...gmtの箇所をキャプチャーしてそこだけにする。
  gmtfile <- sub(".*(wikipathways-.*\\.gmt).*", "\\1",  y[grep('File', y)])
  
  # organism引数の半角スペースをアンダーバーに変換。
  organism <- sub(" ", "_", organism)  
  # 指定した生物種のURLを作成
  url <- paste0(wpurl, gmtfile[grep(organism, gmtfile)])
  
  # URLからgmtファイルをダウンロード  
  f <- tempfile(fileext = ".gmt")
  dl <- utils::download.file(url, method = "libcurl", destfile = f)
  if (is.null(f)) {
    message("fail to download wikiPathways data...")
    return(NULL)
  }
  # gmtファイルをdata.frameに変換
  wp2gene <- read.gmt.wp(f, output = "data.frame")  
  
  ## enrichrの引数用にデータ成型  
  # wikipathay IDと遺伝子IDのdata.frame
  wpid2gene <- wp2gene %>% dplyr::select(.data$wpid, .data$gene) 
  # wikipathway IDとpathway nameのdata.frame
  wpid2name <- wp2gene %>% dplyr::select(.data$wpid, .data$name)
  
  return(list(WPID2GENE = wpid2gene,
         WPID2NAME = wpid2name))
}

# wikipathwayのgenesetを取得
wpdata <- prepare_WP_data(organism = "Homo sapiens")

# enrichrでcompareCluster
wp <- compareCluster(geneClusters = genelist,
                     fun = "enricher", 
                     TERM2GENE = wpdata$WPID2GENE,
                     TERM2NAME = wpdata$WPID2NAME)
wp <- setReadable(wp, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

Discussion

ログインするとコメントできます