7.クラスタリング

• 階層的クラスタリング(Hierarchical Clustering):
  • 凝集型(agglomerative)：各データのみを含むクラスタを併合してゆくことでデータの階層構造を作成する。
    • 最短距離法
    • 最長距離砲
    • 群平均法
    • ウォード法
  • 分割型(divisive)：データ集合を分割することで階層を生成する。
• 非階層クラスタリング：分割の良さを評価する関数（同質性(homogeneity)と分離性(separation)）を定め、良い分割を探す。最適解はNP困難なので準最適解を探す。
  • クラスタの評価基準(Cluster validity Index)
    • 教師なし
      • クラスタ内誤差平方和（エルボー法で極度にSSEが落ちるクラスタ数が良い。） $$\operatorname{SSE} = \sum_k\sum_{x\in C_k}\|x-m_i\|^2$$
      • Calinski-Harabasz index $$\operatorname{VRC} = \frac{\sum_k n_k\|m_k-m_{all}\|^2}{\sum_j\sum_{x\in C_k}\|x-m_k\|^2}\times\frac{N-k}{k-1}$$
      • Dunn index $$\operatorname{DI} = \underset{i=1\cdots n_C}{\min}\left\{\underset{j=i+1\cdots n_C}{\min}\left\{\frac{\underset{x\in C_i,y\in C_j}{\min\{d(x,y)\}}}{\underset{k=1\cdots n_C}{\max}\left\{\underset{x,y\in C_k}{\max}\{d(x,y)\}\right\}}\right\}\right\}$$
      • Silhouette coefficient $$\operatorname{Si} = \frac{b_i-a_i}{\max\{a_i,b_i\}}$$
        • $a_i$ はデータ $i$ と同一クラスタ内のデータとの距離の平均
        • $b_i$ はデータ $i$ が含まれないクラスタのデータとの距離の平均のうち最小値
        • 高い値ほど含まれるクラスタに一致し、他のクラスタと離れていると判定。
        • シルエット分析：個々のデータのシルエット係数を計算し、各クラスタ別に高い順にプロットされたシルエット図を作成する。
        • Shilhouette index：$\operatorname{Sil} = 1/N\sum_i\operatorname{Si}$
    • アルゴリズムの例
      • $k$-means clustering：
        1. $k$個のクラスタの重心との距離から各点を分類。
        2. 更新されたクラスタから再度重心を求める。
        3. 重心が収束するまで1,2を繰り返す。
      • DBSCAN：データを隣接点の数から3種に分類しクラスタリングする。
      • Mean-shift：密度分布の極大を探す。一定半径の球の中の平均値へ移動。
    • 教師あり
      • クラスタリング結果を $C = \{C_1,C_2,\cdots,C_k\}$
      • 正解のクラスタを $A=\{A_1,A_2,\cdots,A_k\}$
      • $C_i$ に含まれる $A_j$ の要素の個数を $X_ij=|C_i\cap A_j|$ とする。
      • エントロピー：
        • $C_i$ のエントロピー $E_i$ は $$E_i = -\sum_{k}\operatorname{Prob}(A_k|C_i)\log \operatorname{Prob}(A_k|C_i)$$
        • 全体のエントロピーは以下で表される。低いほど良いと判定。 $$\sum_k\frac{|C_k|}{N}E_k$$
      • 純度：ある正解のクラスタをどの程度含むか？
        • $C_i$ の純度 $P_i$ は $$P_i = \frac{1}{|C_i|}\max_k|X_{ik}|$$
        • 全てのデータの純度は以下で表される。高いほど良いと判定。 $$\sum_k\frac{|C_k|}{N}P_k$$
      • F値：
        • 再現度 $R_{hk}$ と精度 $P_{hk}$ を以下のように定義 $$R_{hk} = \frac{\left|A_h\cap C_k\right|}{|A_h|}\quad P_{hk} = \frac{\left|A_h\cap C_k\right|}{|C_k|}$$
        • $A_h$ と $C_k$ のF値 $F_{hk}$ は $R_{hk}$ と $P_{hk}$ の調和平均 $$F_{hk} = \frac{2R_{hk}P_{hk}}{R_{hk} + P_{hk}}$$
        • 全体のF値は $A_h$ に対して最大になる $k$ を求めて $$F = \sum_{h=1}^K\frac{|A_h|}{N}\max_kF_{hk}$$

$k$-means法¶

　まずは、$k$-means clusteringの実装を行う。

• 長所：
  • 原理がわかりやすい。
  • 1ステップの計算量が $O(kN)$
• 短所：
  • 「次元の呪い」
    • 高次元において点が疎になりやすい。
    • 次元 $d$ が大きくなるほど、データは中心点から離れ、点の分布は疎になる。
      • 半径 $r$ の $d$ 次元ユークリッド球面の体積は $$V_d(r) = \frac{\pi^{d/2}}{\Gamma\left(\frac{d}{2} + 1\right)}r^d$$
      • 各辺の長さが $2r$ の超立方体の体積は $$(2r)^d$$
      • したがって、超球が超立方体に占める割合は $$V_d(r)/(2r)^d = \frac{\pi^{d/2}}{2^d\Gamma\left(\frac{d}{2} + 1\right)}\underset{d\rightarrow\infty}{\longrightarrow}0$$
  • クラスタ数をあらかじめ与える必要がある。
  • データが超球状に分布することを仮定している。（ユークリッド二乗距離で測るのならば）
  • 初期代表点の選び方に依存する。
    • KZZ法(Katsavounidis 1994)：最も離れた点を選ぶ。
    • K-means++(David Arthur 2007)：最近接の中心から遠いほど選ばれやすい。
      1. ランダムに1点選び、中心とする。
      2. 各点 $x$ に対し、最近接の中心との距離を $D(x)$ とし点 $x_i$ を確率 $\frac{D(x_i)^2}{\sum_i D(x_i)^2}$ で新たな中心として選ぶ。
      3. 2を繰り返して $k$ 個の中心を初期値として選ぶ。

$k$-means法の亜種¶

• kernel $k$-means：$k$-meansの評価関数をカーネルを使用して行う。非線形関数でデータを高次元に写像してクラスタを作成することが可能。
• spectral clustering：データからグラフを作成し、グラフの行列表現の固有値を解き、固有ベクトルのクラスタリングをすることで分類する手法。（重み付きkernel $k$-meansとnormalized cutのspectral clusteringが等価であることが示されている → Kernel k-means, Spectral Clustering and Normalized Cuts）

Q.1 Lloydのk-means実装

1. $K$ 個のクラスタにデータを分類するプログラムの実装
   • 入力：$K$ 個のクラスタにデータを分類する実数データ
   • 出力：データの各クラスタへの分類結果
2. 乱数データを生成し、「データの性質（次元・データ数・クラスタ数）」と「計算時間・誤差二乗平均・ループ回数」の関係性の測定
   • 超球内に一様分布したデータをクラスタリングする。
   • ランダムに設定した複数の中心の、超球体により構成されたデータで計測
3. k-means++の実装

解答

1¶

2¶

※ 超球体内部に一様分布する乱数は、ギブスサンプリングにより生成する。

なお、これから先の実行時間については、以下のプログラムファイルを用いた。

※ jupyter notebook上でプログラムを実装すると面倒なので、プログラムの実行結果をDataFrame で表示する。

3¶

$k$-means法の高速化¶

• $k$-means clusteringにおいて、最も計算時間を取るのは各点と代表点との距離計算（$O(kN)$）
• 繰り返しステップが進むと代表点の移動距離は減り、各点と最も近い代表点は変化しなくなる。
• この不要な距離計算を行わないことで計算の高速化を行う → 三角不等式(Triangle Inequality)
  • 「中心 $c$ の移動距離」 + 「移動前の中心 $c$ との距離」 >= 「移動後の中心 $c^{\prime}$ との距離」より、$\operatorname{dis}\left(\mathbf{c}_p^{\prime},\mathbf{x}\right) < \operatorname{dis}\left(\mathbf{c}_q^{\prime},\mathbf{x}\right)$ の判断が $c^{\prime}$ との距離」**より、$\operatorname{dis}\left(\mathbf{c}_p,\mathbf{x}\right) < \operatorname{dis}\left(\mathbf{c}_q,\mathbf{x}\right)$ を用いることでできる。
  • $50$ 次元以上で高速（Elkanの命題）：Charles Elkan. Using the triangle inequality to accelerate k-means. ICML, pages 147–153. AAAI Press, 2003.
  • $50$ 次元以下で高速（Hamerlyの命題）：Hamerly, G. Making k-means even faster. SDM (Symposium on Data Mining) 130-140, 2010.

Hamerlyの命題¶

　点 $x_i$ の

• 最近接の重心への距離の上界の一つを $u_i\geq\operatorname{d}\left(x_i, c(\operatorname{cls}(i))\right)$
• 2番目に近い重心への距離の下界の一つを $l_i\leq\min_{j\neq \operatorname{cls}(i)}\operatorname{d}\left(x_i, c(j)\right)$
• 所属するクラスターを $\operatorname{cls}(i)$
• $\operatorname{cls}(i)$ の重心を $c\left(\operatorname{cls}(i)\right)$
• $c\left(\operatorname{cls}(i)\right)$ と最も近い他のクラスタの重心との距離を $s\left(\operatorname{cls}(i)\right) = \min_{j\neq \operatorname{cls}(i)}\left\{\operatorname{d}\left(c\left(\operatorname{cls}(i)\right), c\left(j\right)\right)\right\}$

とする。この時、

1. $u_i\leq l_i$
2. $u_i\leq s\left(\operatorname{cls}(i)\right)/2$

のどちらかが成立するならば、$x_i$ に最も近いクラスターは重心移動後も $\operatorname{cls}(i)$ のまま。

証明¶

1. $u_i\leq l_i$ が成立する時： $$ \begin{aligned} \operatorname{d}\left(x_i, c\left(\operatorname{cls}(i)\right)\right)&\leq u_i &\left(\because \text{ definition.}\right)\\ u_i&\leq l_i & \left(\because \text{ definition.}\right)\\ l_i&\leq\min_{j\neq \operatorname{cls}(i)}\operatorname{d}\left(x_i, c(j)\right) &\left(\because \text{ hypothesis.}\right) \end{aligned} $$
2. $u_i\leq s\left(\operatorname{cls}(i)\right)/2$ が成立する時： $$ \begin{aligned} 2u_i&\leq s\left(\operatorname{cls}(i)\right) &\left(\because \text{ hypothesis.}\right)\\ s\left(\operatorname{cls}(i)\right)&\leq\operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right) + \min_{j\neq \operatorname{cls}(i)}\operatorname{d}\left(x_i,c(j)\right) & \left(\because \text{ Triangle Inequality.}\right)\\ \operatorname{d}\left(x_i, c\left(\operatorname{cls}(i)\right)\right) &\leq u_i & \left(\because \text{ definition.}\right)\\ \therefore 2u_i&\leq u_i + \min_{j\neq \operatorname{cls}(i)}\operatorname{d}\left(x_i,c(j)\right) \end{aligned} $$

より、どちらの場合でも以下が成立する。

Q.2 高速化手法の実装

1. Hamerlyの高速化手法の実装
   • 入力：$K$ 個のクラスタにデータを分類する実数データ
   • 出力：データの各クラスタへの分類結果
   • 初期代表点が同じ場合Lloydのアルゴリズムと等しいクラスタを出力することを確認。
2. 乱数データを生成し、「データの性質（次元・データ数・クラスタ数）」と「計算時間・誤差二乗平均・ループ回数」の関係性の測定
3. Elkanの高速化手法の実装

解答

1¶

2¶

3¶

Elkanの命題¶

　点 $x_i$ の

• 最近接の重心への距離の上界の一つを $u_i\geq\operatorname{d}\left(x_i, c(\operatorname{cls}(i))\right)$
• $j$ 番のクラスタの重心への距離の下界を $l_{ij}\leq\operatorname{d}\left(x_i, c(j)\right)$
• 所属するクラスターを $\operatorname{cls}(i)$
• $\operatorname{cls}(i)$ の重心を $c\left(\operatorname{cls}(i)\right)$
• $A$ の、一つ前のイテレーション時の値を $A^{\prime}$

とする。この時、

1. $\frac{1}{2}\operatorname{d}\left(c\left(\operatorname{cls}(i)\right),c(j)\right)\geq u_i$ が成立する時、$\operatorname{d}\left(x_i,c(j)\right)$ を計算する必要がない。
2. $l_{ij} = \max\left\{0, l_{ij}^{\prime}-\operatorname{d}\left(c(j), c^{\prime}(j)\right)\right\}$ と下界を計算できる。

※ Hamerlyの手法 では、全てのクラスタ重心への距離の下界の一つ $l_i$ を使用するのに対し、ここではそれぞれのクラスタ重心に対して距離の下界 $l_{ij}$ を保持している

証明¶

1. $\frac{1}{2}\operatorname{d}\left(c\left(\operatorname{cls}(i)\right),c(j)\right)\geq u_i$ が成立する時 $\left(j{\neq\operatorname{cls}(i)}\right)$： $$ \begin{aligned} \frac{1}{2}\operatorname{d}\left(c\left(\operatorname{cls}(i)\right),c(j)\right)\geq u_i&\geq\operatorname{d}\left(x_i, c(\operatorname{cls}(i))\right) &\left(\because\text{ definition.}\right)\\ \underset{-\frac{1}{2}\operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right)}{\Longleftrightarrow}\frac{1}{2}\operatorname{d}\left(c\left(\operatorname{cls}(i)\right),c(j)\right) - \frac{1}{2}\operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right)&\geq\frac{1}{2}\operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right)\\ \operatorname{d}\left(x_i,c(j)\right)&\geq\operatorname{d}\left(c\left(\operatorname{cls}(i)\right),c(j)\right) - \operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right)&\left(\because\text{ Triangle Inequality.}\right)\\ \therefore\operatorname{d}\left(x_i,c(j)\right) &\geq\operatorname{d}\left(x_i,c\left(\operatorname{cls}(i)\right)\right) \end{aligned} $$
2. この下界は、$l_{ij}^{\prime}$ が良い下界で、$j$ 番のクラスタの重心の移動距離が小さいほど良い推定となる。(trivial) $$ \begin{cases} \begin{aligned} \operatorname{d}\left(x_i,c(j)\right) &\geq \operatorname{d}\left(x_i,c^{\prime}(j)\right) -\operatorname{d}\left(c(j), c^{\prime}(j)\right) & \left(\because\text{ Triangle Inequality.}\right)\\ \operatorname{d}\left(x_i,c^{\prime}(j)\right)&\geq l_{ij}^{\prime} & \left(\because\text{ definition.}\right)\\ \end{aligned} \end{cases}\\ \begin{aligned} \operatorname{d}\left(x_i,c(j)\right) &\geq\max\left\{0, \operatorname{d}\left(x_i,c^{\prime}(j)\right) -\operatorname{d}\left(c(j), c^{\prime}(j)\right)\right\}\\ &\geq\max\left\{0,l_{ij}^{\prime} - \operatorname{d}\left(c(j), c^{\prime}(j)\right)\right\} \end{aligned} $$

1細胞RNA-seq(scRNA-seq)¶

• 微量なサンプルからcDNAライブラリーを作成する技術の発展により、1細胞単位での遺伝子発現量データを得ることが可能となった。（がん・メタゲノム・細胞系譜・神経生物学など。）
• 各細胞が数千~数万種類の遺伝子発現量のデータを持つscRNA-seq dataを解析するための手法として以下の手法がある。
  • 2次元または3次元に次元削減することによる可視化
  • 次元削減されたデータをクラスタリングすることによる細胞集団や細胞系譜の同定
• 例として、sc-RNAseqにより線虫(c.elegans)の原腸形成から最終分化細胞に至る発達段階の胚のトランスクリプトームを分析し、「sc-RNAseqの結果から細胞が分化する過程の細胞系譜木の再構成」をおこなったものがある（"A lineage-resolved molecular atlas of C.elegans embryogenesis at single-cell resolution"）

workflow¶

Single-cell RNA sequencing workflow	Cell isolation

Library Preparation¶

Full-length RNA-seq from single cells using Smart-seq2	Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets

次元削減法¶

scRNA-seqの遺伝子発現マトリックスでは各Sampleが高次元の情報を持つためのクラスタリングを行う前に次元削減が必要

• PCA:
  • 「データの分散最大化」 or 「次元削減による誤差の最小化」から主成分を求める。
  • 線形変換
• t-SNE:
  • データ点間の類似度を確率で表し低次元に埋め込む
  • 非線形変換
• UMAP:
  • 多様体学習を使用した非線型変換

※ scRNA-seqでは、t-SNE, UMAPが可視化のためによく用いられる。

t-SNE¶

• t-distribution: 低次元へ変換したデータ間の類似度を表すために自由度1のt-分布を使用する。
• Stochastic Neighbor Embedding: 高次元空間でのデータ間のユークリッド距離類似度を確率で表現し低次元に変換する。

「データ点間の類似度を表す確率分布」と「低次元に変換した点間の類似度を表す確率分布」を求め、その間のKullback-Leibler divergenceが最小化するように勾配法によって変換後の点の座標を求める手法。

確率分布¶

• 高次元のデータ点 $x_i$ に対するデータ点 $x_j$ の類似度を対象化された結合確率(joint distribution)で表す。 $$p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}$$
• データ点 $x_j$ は点 $x_i$ を中心とする標準偏差 $\sigma_i$ のガウス分布に比例して選び、条件付確率 $p_{j|i}$ は以下のように求める。 $$p_{j|i} = \frac{\exp\left(-\left\|x_i-x_j\right\|/2\sigma_i^2\right)}{\sum_{k\neq i}\exp\left(-\left\|x_i-x_k\right\|/2\sigma_i^2\right)}$$
• 低次元に変換した点 $y_i$ と点 $y_j$ の類似度を自由度$1$のt分布の同時確率で表す。 $$q_{ij} = \frac{\left(1 + \left\|y_i-y_j\right\|^2\right)^{-1}}{\sum_{k\neq l} \left(1 + \left\|y_k-y_l\right\|^2\right)^{-1}}$$

の$p_{ij}$ と $q_{ij}$ のミスマッチを最小化する。

Perplexity¶

全ての点のガウス分布に対して同じ標準偏差 $\sigma$ を使用するのではなく、「密度の高い領域の点は小さい標準偏差」を「密度の低い領域では大きい標準偏差」を用いる。

ここで、上記の条件を満たすために、データ点 $x_i$ 周辺のデータ密度を表現するPerplexityを以下のように定義する。

ここで、

• $\text{Perp}(P_i)$ に対して $H(P_i)$ は単調増加
• $H(P_i)$ に対して標準偏差 $\sigma_i$ は単調増加

より、ガウス分布の標準偏差 $\sigma_i$ が大きいほど、距離が遠い点でも条件付き確率 $p_{j|i}$ が大きな値(negligibleでない)を取ることが多くなる。

Pros and Cons.¶

• 利点：
  • Stochastic Neighbor Embeddingからの改善点であり、t-分布の裾が重いため低次元の空間において距離の遠いデータを遠くに配置できる。
  • データ間の近さに対称性がある
• 欠点：
  • 2,3次元へのマッピング以外では計算コストが高い。
  • 局所構造が次元の呪いの影響を受けやすい。
  • 最適化アルゴリズムの収束性が保証されていない。
  • perplexityの選択により出力に大きな影響がある。
  • マッピング後のクラスタ間の距離関係は元の距離を必ずしも反映していない。

DBSCAN¶

• 空間の中で密集している点をクラスタとしてまとめ、疎な領域の点を外れ値とするクラスタリング手法。
• 半径 $\varepsilon$、最小点数 $\text{minPts}$ が与えられた時に、点を $3$ 種類（「コア点」・「到達可能」・「外れ値」）に分類してクラスタリングを行う。
  • コア点：点 $P$ を含め半径 $\varepsilon$ 以内に $\text{minPts}$ 個以上の点が存在する点。
  • 直接到達可能：コア点 $P$ から $\varepsilon$ 以内にある点 $q$
  • 到達可能：コア点 $P_i$ から直接到達可能なコア点 $P_{i+1}$ のパス $P_1P_2\cdots P_nq$ が存在するとき、点 $q$ はコア点 $P_1$ から到達可能。
  • 外れ値：どの点からも到達可能でない点

アルゴリズム¶

def DBSCAN(data, epsilon, minPts):
    """
    @params data   : 点集合
    @params epsilon: 半径`epsilon`以内の点を直接到達可能と定義する。
    @params minPts : 半径`epsilon`以内に`minPts`個以上の点が存在する場合、コア点とする。
    """
    Clusters = []
    while(data):
        # 未探索の点Pを選ぶ。
        p = point_lists.pop()
        if p.searched:
            continue
        # 点pから半径ε以内にある点の個数。
        Ps_NeighborPts = p.reacheable(epsilon)
        # minPts個以上の点が存在する場合
        if minPts <= len(Ps_NeighborPts):
            # 点Pを新たなクラスタCに追加する。
            cluster = [p]
            for q in Ps_NeighborPts:
                # 未探索の場合、探索済みにする。
                if not q.searched:
                    q.searched = True
                # 点qから半径ε以内にある点の個数。
                Qs_NeighborPts = q.reacheable(epsilon)
                # minPts個以上の点が存在する場合
                if minPts <= len(Qs_NeighborPts):
                    # 点qから直接到達可能な点の集合を追加する。
                    Ps_NeighborPts += Qs_NeighborPts
                if q not in cluster:
                    cluster.append(q)
            Clusters.append(cluster)
        # minPts個未満の点しか存在しない場合
        else:
            # 探索済みとして終了。
            p.searched=True
    return Clusters

Pros and Cons.¶

• 利点：
  • k-meansと異なり、あらかじめクラスタ数を与える必要がない。
  • クラスタが超球状であることを仮定していない。
• 欠点：
  • クラスタ間で密度に大きな差がある場合正しくクラスタリングされない。
  • パラメータの選択により出力が影響を受ける。
  • 高次元では次元の呪いにより半径 $\varepsilon$ の選択が困難。

Q.3 t-SNE+DBSCAN

1. 点の集合をDBSCANを用いてクラスタに分類するプログラムの実装
   • 入力：データ集合・半径 $\varepsilon$・最小の近接数 minPts
   • 出力：クラスタへの分類結果
2. 線虫の細胞系譜の論文 "A lineage-resolved molecular atlas of C. elegans embryogenesis at single-cell resolution" のgene expression matrix dataのクラスタリング: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126954
   • GSE126954_cell_annotation.csv は、各細胞の情報（cell type, embryo time, lineage など。）
   • GSE126954_gene_annotation.csv は、各遺伝子の情報（wormbase での ID と short name.）
   • GSE126954_gene_by_cell_count_matrix.txt.gz は、遺伝子発現マトリックスのデータであり、疎な行列のため、圧縮した書き方をしている。
3. その他の細胞種のデータを使用してt-SNEを実行し、DBSCANでクラスタリングを行う。
   1. GSE126954_gene_by_cell_count_matrix.txtから特定のcell typeのデータを取り出す
   2. PCAによる次元削減
   3. t-SNEまたはUMAPによる次元削減
   4. DBSCANによるクラスタリング
   5. 各細胞種での遺伝子発現量のヒートマップを作成する

解答

GSE126954_gene_by_cell_count_matrix.txt のデータをクラスタリングする。

• 細胞総数は$89701$個
• 遺伝子総数は$20222$個
• celltypeは$37$種類
• embryo timeのbinは$12$種類

特定の細胞種に絞って解析を行う。

Decomposition¶

PCA¶

tSNE¶

UMAP¶

Results¶

Clustering¶

tSNE¶

UMAP¶

Results¶