MBCluster.Seqで遺伝子間クラスタリング

Clustering RNA-seq R

RNA-seqデータのクラスタリング

RNA-seqでは得られたリードを遺伝子やエキソンごとに数え上げることで発現解析を行うことは一般的なワークフローの一つです。得られたカウントデータは通常、行に遺伝子またはエキソン、列にサンプルが並ぶような形の数値行列として表現されます。つまり、行数は遺伝子またはエキソンの数、列数はサンプル数となります。

今回はこのデータのクラスタリングを考えていこうと思います。RNA-seqデータのクラスタリングでは、まずどちらの向きでデータ単位を考えるかで2つに分けられます。つまり列ごとに縦に並ぶデータを1つのデータ単位と見てサンプル間のクラスタリングを行うか、もしくは、行ごとに横に並ぶデータを1つのデータ単位とみて遺伝子間のクラスタリングを行うかです。

今回は遺伝子間のクラスタリングを行うことを考えてみます。

クラスタリング手法

クラスタリング手法には様々な方法があります。k-means やら階層的クラスタリングやらありますが、今回はModel-based Clusteringによるクラスタリングを行ってみたいと思います。Model-based Clusteringでは統計モデルに基づきクラスタリングを行う方法で、この方法の1つを実装したRパッケージとして MBCluster.Seq というものがあります。

以下実際に MBCluster.Seq を試した結果です。詳細は論文とマニュアルを参照してください。

# ライブラリ
library(MBCluster.Seq)
library(ggplot2)
library(tidyr)
library(dplyr)

# 今回用いるデータ(MBCluster.Seq同梱)
# 実験毎にコードを変える部分
data("Count")
head(Count)
##      N1.1 N2.1 N3.1 N4.1 N1.2 N2.2 N3.2 N4.2
## [1,]    2    0    0    0    4    0    0    0
## [2,]    4  357 2537 1295   19 1056 2690 4411
## [3,]    0    0    6    8    1    2    8   18
## [4,]    1    1    1    0    2    5    1    2
## [5,]    2   10  107   32    2   31   94   69
## [6,]   79    8   18    5  102   24   21   14
GeneID <- 1:nrow(Count)
Treatment <- c(1,2,3,4,1,2,3,4)
nstart <- 5

# 遺伝子数、クラスタ数など
G <- nrow(Count)
K <- 2:50
I <- 4

# logFC や dispersion を計算
d <- RNASeq.Data(Count, NULL,  Treatment)

# クラスタリング
res <- lapply(K, function(k) {
  # 各 k について nstart 回繰り返し対数尤度が最大となるクラスタリング結果を採用
  # （R の kmeans 関数の nstart 引数を参考にした）
  all_res <- lapply(seq_len(nstart), function(i) {
    # centers の初期化
    c0 <- KmeansPlus.RNASeq(d, k, "nbinom")
    # EMアルゴリズム
    res <- Cluster.RNASeq(d, "nbinom", c0$centers)
    # 対数尤度の計算
    res$logl <- lglk.cluster(d, "nbinom", res$cluster)
    res
  })
  all_logl <- sapply(all_res, "[[", "logl")
  res <- all_res[[which.max(all_logl)]]
  res
})

# 対数尤度から AIC を計算する
logl <- sapply(res, "[[", "logl")
aic <- -2 * (logl * G - G * (K + 1) - K * I + 1)
plot(K, aic, type = "o", col = ifelse(seq_along(aic) == which.min(aic), "red", "black"))

f:id:menseki:20201004231347p:plain

# AIC が最小となるクラスタ数 k を選択
res.obj <- res[[which.min(aic)]]

# 図の描画
logfc <-
  as.data.frame(d$logFC) %>% 
  mutate(cluster = res.obj$cluster) %>% 
  tibble::rowid_to_column() %>% 
  pivot_longer(-c(cluster, rowid), names_to = "group", values_to = "logFC")

center <- 
  unclass(res.obj$center) %>% 
  as.data.frame() %>% 
  mutate(cluster = row_number()) %>% 
  tibble::rowid_to_column() %>% 
  pivot_longer(-c(cluster, rowid), names_to = "group", values_to = "logFC")
  
ggplot(logfc, aes(group, logFC)) +
  geom_line(aes(group = rowid)) +
  facet_wrap(vars(cluster)) +
  geom_line(data = center, aes(group, logFC, group = rowid), col = "red")

f:id:menseki:20201004231424p:plain