Upgrade to Pro — share decks privately, control downloads, hide ads and more …

oku-slide-20240802

 oku-slide-20240802

幾つかのデータ解析手法の紹介
奥 牧人
2024/08/02
富山大学ムーンショット研究数理セミナー

Makito Oku

July 30, 2024
Tweet

More Decks by Makito Oku

Other Decks in Research

Transcript

  1. 何について学ぶべきか データ解析手法について学ぶとき、必ずしも 動作原理 や背後に ある 数学の概念 を理解する必要はない。 むしろ、実際のデータ解析で使用する際の 注意点 や、手法毎の

    利点と欠点を踏まえた 使い分け を覚える方が大事。 自動車で例えると、エンジンの仕組みを詳細に理解していなくて も運転は出来るが、交通ルールを分かっていないと事故を起こし てしまう。 9 / 60
  2. 機械学習の区分のおさらい 大区分 中区分 例 教師あり学習 分類 線形判別分析、ロジスティック 回帰、ランダムフォレスト、 LightGBM 教師あり学習

    回帰 単回帰、重回帰、Lasso, ガウス 過程回帰 教師なし学習 クラスタリング K平均法、階層的クラスタリング 教師なし学習 次元削減 PCA, t-SNE, UMAP 本発表では 教師なし学習のみ 扱う。 12 / 60
  3. 主成分分析 項目 説明 名前 主成分分析、Principal Component Analysis (PCA) 区分 教師なし学習

    > 次元削減 処理内容 多次元のデータを分散の意味で最もよく見える方向から 眺めた時の射影を計算する 関連用語 寄与率、ローディング、スクリープロット 利点 簡単、色々なことが分かる、解釈しやすい 欠点 scRNA-seqなど多数のグループがある場合に向かない 注意点 3次元は静止画だと意味がない (最適でない2次元への 射影になってしまっているため)、各軸毎に反転可 15 / 60
  4. t-SNE 項目 説明 名前 t-SNE 区分 教師なし学習 > 次元削減 処理内容

    多次元のデータ点を上手に2次元上に配置する 関連用語 - 利点 多数のグループをきちんと分離出来る 欠点 本当は同じグループの点が分離する場合がある、利用率 がUMAPに押され気味 注意点 sklearnの場合、学習率をデータ点数/12にし、初期配 置をPCAで決めるようにする。また、必要に応じて事前 にPCAで10次元程度に落としてからt-SNEをかける。 17 / 60
  5. UMAP 項目 説明 名前 UMAP 区分 教師なし学習 > 次元削減 処理内容

    多次元のデータ点を上手に2次元上に配置する 関連用語 - 利点 多数のグループをきちんと分離出来る 欠点 本当は同じグループの点が分離する場合がある、パラメ ータ調節が必要 注意点 パラメータ (min_dist, spread) を適宜調節する 19 / 60
  6. 倍率変化 項目 説明 名前 倍率変化、Fold change (FC) 区分 発現変動遺伝子解析 処理内容

    2群間の比較において、平均値の比が (対数を取る前の 値で) 2倍より大きい遺伝子を選ぶ 関連用語 - 利点 簡単、再現性が高い、サンプル数に依存しない 欠点 P値至上主義者に理解して貰えない 注意点 閾値は2倍から変えても構わない 20 / 60
  7. FDR 項目 説明 名前 FDR、False Discovery Rate 区分 発現変動遺伝子解析、多重検定補正 処理内容

    N個のP値に対応するN個のQ値を計算し、Q値が指定し たFDRの値以下のものを有意と判定する。FDRとは出力 の中に含まれる偽陽性の割合のこと。 関連用語 Benjamini-Hochberg法 利点 Nが大きい場合にBonferroni補正より検出力が高い 欠点 取れる遺伝子数がサンプル数に依存する 注意点 P値がほぼ同じでもQ値が大きく違ったり、P値が違うの にQ値が全く同じだったりする。P値とQ値は別物。 22 / 60
  8. ベン図 項目 説明 名前 ベン図、Venn diagram 区分 発現変動遺伝子解析 処理内容 複数の集合間の重なり具合を計算して図示する

    関連用語 和集合、積集合、差集合、ジャカード指数 利点 3集合以下なら分かりやすい 欠点 5集合以上だと領域数が多過ぎて意味不明 注意点 4集合の場合は正円だと描けないので楕円等を使う 23 / 60
  9. 階層的クラスタリング 項目 説明 名前 階層的クラスタリング 区分 教師なし学習 > クラスタリング 処理内容

    似ているデータ点を順にまとめていき、閾値で切って 複数のクラスタに分ける 関連用語 樹形図、平均連結法、ウォード法、シルエットスコア 利点 ランダム性がない、生命科学の分野で一般的 欠点 汎用的な閾値自動決定法が無い 注意点 RNA-seqでは相関類似度を使うか遺伝子毎にZスコア化 する、単連結法はChaining現象が発生するため非推奨 24 / 60
  10. 階層的クラスタリング、補足 設定項目 類似度 or 非類似度: ユークリッド距離、相関類似度など 連結法: ウォード法、平均連結法など クラスタ分割基準: クラスタ数指定、非類似度の閾値指定など

    私はRNA-seqデータ解析の場合では相関類似度、平均連結法、 非類似度の閾値指定を使っている。閾値は以下の式で計算 意味は、相関係数をフィッシャー変換した際の3σ点 1 − tanh ( 3 √N − 3 ) 25 / 60
  11. K平均法 項目 説明 名前 K平均法、K-means 区分 教師なし学習 > クラスタリング 処理内容

    クラスタ数をKとし、各クラスタの中心位置 (最初はラ ンダム) とそれに属するデータ点を交互に更新する 関連用語 - 利点 簡単 欠点 クラスタ数を決め打ちする必要がある、クラスタの初期 配置により結果が変わる 注意点 RNA-seqでは相関類似度を使うか遺伝子毎にZスコア化 する 26 / 60
  12. エンリッチメント解析 項目 説明 名前 エンリッチメント解析 区分 その他 処理内容 入力した遺伝子一覧の中に、どのような注釈を持つ遺伝 子が多いかをデータベースに問い合わせ、出力する

    関連用語 GO (Gene Ontology)、フィッシャーの正確検定 利点 個々の遺伝子について知らなくても結果の解釈が出来る 欠点 データベースによって結果が違う、重複数が極めて少な いものも出力される、詳しいことは分からない 注意点 P値だけ見てはダメ、個数も確認する 27 / 60
  13. 同期性揺らぎ遺伝子解析 (DNB解析) 項目 説明 名前 同期性揺らぎ遺伝子解析、DNB解析 区分 その他 処理内容 揺らぎが大きく、かつ、互いに強く相関した遺伝子群を

    選ぶ 関連用語 中央絶対偏差、スピアマンの相関係数、階層的クラスタ リング 利点 揺らぎに注目する点が独自的 欠点 言葉が独り歩きしていて過度に期待される、結果が合っ ているかどうかよく分からない、手法が確立していない 注意点 後述 30 / 60
  14. 考えるべきこと 考えるべきことが沢山ある。 揺らぎの指標: 標準偏差、中央絶対偏差 揺らぎ判定法: 倍率変化、割合、仮説検定 倍率変化や割合の閾値 仮説検定の種類: F検定、Levene ル

    ビ ー ン 検定、Brown-Forsythe検定 対照群: 使用する、使用しない 複数時点の扱い: 時点毎、中心化して結合、対照群のみプール 相関の指標: ピアソンの相関係数、スピアマンの相関係数 連結法: 平均連結法、完全連結法 クラスタ分割基準: 個数、閾値 32 / 60
  15. RNA-seqに対する現在の方針 外れ値に強い 中央絶対偏差 と スピアマンの相関係数 は必須 結果の再現性を高めるため 対照群なし を選択 時点毎に計算するか、中心化して結合するかは、結果があまり

    変わらないことが多いのでどちらでも良い。 対照群なし、かつ、時点毎の場合、異なる時点間で取れる遺伝子 が大きく重複するので、事後クラスタリング が必要 結果の妥当性を以下で確認 ヒートマップ 遺伝子リスト エンリッチメント解析 34 / 60
  16. 何が難しいのか? 手法が確立しておらず、選択肢の組み合わせの数が膨大なため、 データ毎に 数理の知識 と データ解析の経験 に基づき組み合わせ を変え、その妥当性を高度な ドメイン知識 に基づき多角的かつ

    慎重に判断する必要がある。 私自身もマイクロアレイとRNA-seqのドメイン知識しか持ってい ないので、他の種類のデータには全く手が出せない。 今回の目的はお互いに知っている手法を教え合うことだが、正直に言 うと DNB解析の習得は勧められない。やるとしても、再実験が容易 で結果の再現性を確認しやすいものが良いと思う。 35 / 60
  17. よくある間違い 1サンプルのみで高発現し、残りのサンプルでは全然揺らいでい なかった → 外れ値 に引きずられただけ サンプルを番号順にした場合、あるサンプルまでは一様に低く、 それ以降が一様に高かった → バッチ効果

    を拾っただけ snoRNA (Snorから始まる) など短い転写産物ばかり → S/N比 の悪い遺伝子が取れただけ Olfrなど平均発現量の低い遺伝子ばかり → 同上 ミトコンドリア遺伝子 (mtから始まる) やリボソーム蛋白質 (Rp から始まる) ばかり → 単に 高発現 遺伝子を拾っただけ ほぼ全ての条件で同じ遺伝子が出てきた → 単に コントロール群 で偶然値が揃っていた 遺伝子を拾っただけ 36 / 60
  18. エネルギー地形解析 項目 説明 名前 エネルギー地形解析、Energy Landscape Analysis (ELA) 区分 その他

    処理内容 データを二値化し、イジングモデルに当てはめ、エネル ギーの谷や、谷と谷の間のエネルギー障壁を計算する 関連用語 イジングモデル、ベースングラフ、非連結性グラフ 利点 独自性が高い、解析ツールの使用自体は簡単 欠点 約7個の変数を選ぶ必要がある、結果の解釈が難しい 注意点 後述 38 / 60
  19. エネルギー地形解析、補足 江崎先生がMatlabで作ったツールのPython移植版: https://github.com/okumakito/elapy エネルギーは仮想的なもので、確率と関連した値 単独のエネルギーの値に意味はなく、相対的な差が大事 頂点番号は、0と1から成る2進数を10進数に直したもの 例、1100111 = 64 +

    32 + 4 + 2 + 1 = 103 ベースングラフでは、頂点が多いほど谷が広く、最小エネルギー が小さいほど谷が深いことを意味する。 非連結性グラフは、状態間 (正確にはエネルギー極小パターン間) を遷移する際のエネルギー障壁の高さを図示したもの。 p ∝ exp(−E) 2 39 / 60
  20. 主な前処理 FASTQ → BAM → カウントデータ (私は出来ないので外注) ファイルの結合 不要な行や列の除外 ID変換

    列名の変更 正規化 対数化 外れ値の除外 低発現遺伝子の除外 その他、対象データ特有の問題に対する処理 43 / 60
  21. ファイルの結合 データファイルはなるべく1つにまとめる。 作業効率を上げるため、かつ、ミスを減らすため 私の場合、結合したい元ファイルを同じディレクトリに入れ、 以下のような流れで結合している。後述の処理を適宜足す。 dir_name = '../data/orig/' # 元データを置いたディレクトリ名

    file_list = os.listdir(dir_name) # ファイル名一覧を取得 file_list.sort() # 並び替え out_list = [] # 出力格納用リスト for file_name in file_list: file_path = os.path.join(dir_name, file_name) # フルパス df = pd.read_csv(file_path) # csvファイル読み込み ### ここに色々な処理を追加 ### out_list.append(df) # 出力格納用リストに追加 df = pd.concat(out_list, axis=1) # 全体を結合 df.to_csv('tmp.csv') # ファイルに保存、重要ファイル上書き防止のためtmp 44 / 60
  22. ファイル読み込みの補足 データがコンマ区切りでなくタブ区切りの場合 df = pd.read_table(file_path) # read_csvの sep='\t' と同じ 1列目をindexとして扱う場合

    (1行目は既定でheader扱い) df = pd.read_csv(file_path, index_col=0) 一部の列のみを読み込む場合 df = pd.read_csv(file_path, usecols=['ABC','DEF']) # 2列のみ ある文字から始まるコメント行を除く場合 df = pd.read_csv(file_path, comment='#') # '#" から始まる行を除外 他にもread_csv, read_tableのオプションは沢山ある。 45 / 60
  23. 不要な行や列の除外 ここは千差万別。比較的よく使うもの ある列 (仮にxyzとする) を除外 df = df.drop('xyz', axis=1) ある列が空欄でない行だけ抽出

    df = df[~df.xyz.isnull()] ある列の値がある値 (仮にabcとする) の行だけ抽出 df = df[df.xyz=='abc'] ある列の値があるリスト (仮にx_list) に含まれる行だけ抽出 df = df[df.xyz.isin(x_list)] 46 / 60
  24. 不要な行や列の除外、続き ある列の値がある文字列から始まる行のみ抽出、他 df = df[df.xyz.str.startswith('ABC')] # ABCで始まる df = df[df.xyz.str.endswith('ABC')]

    # ABCで終わる df = df[df.xyz.str.contains('ABC')] # ABCを含む (以下、条件部で利用) 文字列のある番号の文字のみ抽出 df.xyz.str[0] # 例、ABC, DEF > A, D 文字列をsplitした場合のある番号のブロックのみ抽出 df.xyz.str.split('-').str[1] # 例、ABC-123, DEF-456 > 123, 456 型変換 df.xyz.astype(int) # 整数に変換 df.xyz.astype(str) # 文字列に変換 47 / 60
  25. ID変換 Ensembl ア ン サ ン ブ ル IDなどの分かりにくいIDを遺伝子記号に変換 例、ENSMUSG00000024401

    → Tnf 複数の元IDが同じ遺伝子記号に対応する場合は最大値を採用 対応する遺伝子記号の無い転写産物のデータは削除 やり方は色々ある。例えば、遺伝子記号の列をindexにして groupbyするなど df = df.set_index('gene_sym').groupby(axis=0,level=0).max() 48 / 60
  26. 正規化 マイクロアレイでは分位数正規化 (quantile normalization) など も使われていたが、RNA-seqデータではほぼTPM TPMでは、まず行毎に正規化し、次に列毎に正規化する。 df= df.divide(len_sr, axis=0)

    # 各行を遺伝子の長さで割る df = df / df.sum() # 列毎に総和で割る df *= 10**6 # 1Mを掛ける ちなみに、RPKMでは行と列と正規化の順が逆 TPMで箱ひげ図がうまく揃わない場合、私は刈り込みTPMを使っ ている。2行目を以下のように変更 from scipy import stats trim_sum = stats.trim_mean(df, 0.05) * len(df) # 刈り込み総和 df = df / trim_sum # 列毎に刈り込み総和で割る 50 / 60
  27. 対数化 0は対数を取れない。 1などを予め足す場合もある。 私は、0は0のままとし、0以外の最小値を0にシフトしている。 df[df==0] = None # 0を一旦Noneに置換 df

    = df.apply(np.log2) # 対数化 df -= df.min().min() # 0以外の最小値を引く df = df.fillna(0) # Noneを0に戻す 0以外の場合を式で書くと . ただし、どうしても0と0以外の小さい値の差が不自然に大きく なる問題が生じる。当然、PCAやDEGに悪影響が出る。 結局、対数を取るなら後述する低発現遺伝子の除外が必要 log(x/xmin) 51 / 60
  28. PCAの使い方 色分け表示するにはseabornのscatterplotが便利 座標データをデータフレーム化し、各サンプルのグループ名を 3列目にして、hueでその列名を指定する。 from sklearn.decomposition import PCA model =

    PCA(n_components=2) # PCAオブジェクトを用意 score = model.fit_transform(df.T) # スコア行列の計算 pos_df = pd.DataFrame(score, columns=list('xy')) # データフレーム化 pos_df['group'] = group_list # 各サンプルのグループ名を追加 perc1, perc2 = 100 * model.explained_variance_ratio_ # 寄与率 fig, ax = plt.subplots() # 図を用意 sns.scatterplot(data=pos_df, x='x', y='y', hue='group', ax=ax) # 散布図の描画、色はグループ別 ax.set_xlabel(f'PC1 ({perc1:.1f} %)') # x軸のラベル、寄与率込み ax.set_ylabel(f'PC2 ({perc2:.1f} %)') # y軸のラベル、寄与率込み ax.legend(bbox_to_anchor=(1,1)) # 凡例の位置を枠外に移動 fig.show() # 図の表示 53 / 60
  29. ヒートマップの使い方 RNA-seqデータをヒートマップにする時は遺伝子毎のZスコア化 (標準化) が必須 clustermapの場合 (簡単だが微調整できない) sns.clustermap(df, z_score = 0,

    # 各行 (遺伝子) をZスコア化 cmap = 'RdBu_r', # 高発現が赤、低発現が青 vmax = 3, # カラースケールの最大値 vmin = -3, # カラースケールの最小値 row_cluster = True, # 各行をクラスタリングする col_cluster = False, # 各列をクラスタリングしない metric = 'correlation', # 相関係数を類似度とする method = 'average' # 平均連結法を使用 ) 微調整したい場合はsns.heatmapを使う。 54 / 60
  30. ヒートマップの使い方2 sns.heatmapを使う場合 from scipy.cluster.hierarchy import linkage, leaves_list # 行毎にZスコア化 df

    = ((df.T - df.T.mean()) / df.T.std()).T # 行方向をクラスタリング、相関類似度と平均連結法を使用 Z = linkage(df, metric='correlation', method='average') # 行方向を並べ替え df = df.iloc[leaves_list(Z)] # 図の描画 fig, ax = plt.subplots() # 図を用意 sns.heatmap(df, cmap='RdBu_r', ax=ax, vmax=3, vmin=-3) # ヒートマップの描画 fig.show() # 図の表示 55 / 60