Upgrade to Pro
— share decks privately, control downloads, hide ads and more …
Speaker Deck
Features
Speaker Deck
PRO
Sign in
Sign up for free
Search
Search
oku-slide-20240802
Search
Makito Oku
July 30, 2024
Research
0
67
oku-slide-20240802
幾つかのデータ解析手法の紹介
奥 牧人
2024/08/02
富山大学ムーンショット研究数理セミナー
Makito Oku
July 30, 2024
Tweet
Share
More Decks by Makito Oku
See All by Makito Oku
oku-slide-20231129
okumakito
0
99
oku-slide-20230827
okumakito
0
110
oku-slide-20230213
okumakito
0
200
oku-slide-20221212
okumakito
0
75
oku-slide-20221129
okumakito
0
140
oku-slide-20221115
okumakito
0
290
oku-slide-20220820
okumakito
0
260
oku-slide-stat1-1
okumakito
0
270
oku-slide-stat1-2
okumakito
0
300
Other Decks in Research
See All in Research
Practical The One Person Framework
asonas
1
1.4k
湯村研究室の紹介2024 / yumulab2024
yumulab
0
240
TransformerによるBEV Perception
hf149
1
380
大規模言語モデルのバイアス
yukinobaba
PRO
4
670
[CV勉強会@関東 CVPR2024] Visual Layout Composer: Image-Vector Dual Diffusion Model for Design Layout Generation / kantocv 61th CVPR 2024
shunk031
1
410
LLM based AI Agents Overview -What, Why, How-
masatoto
2
590
Physics of Language Models: Part 3.1, Knowledge Storage and Extraction
sosk
1
930
MetricSifter:クラウドアプリケーションにおける故障箇所特定の効率化のための多変量時系列データの特徴量削減 / FIT 2024
yuukit
2
110
marukotenant01/tenant-20240916
marketing2024
0
410
RCEへの近道
kawakatz
1
800
MIRU2024_招待講演_RALF_in_CVPR2024
udonda
1
330
ミニ四駆AI用制御装置の事例紹介
aks3g
0
160
Featured
See All Featured
We Have a Design System, Now What?
morganepeng
50
7.2k
The Power of CSS Pseudo Elements
geoffreycrofte
72
5.3k
Fireside Chat
paigeccino
32
3k
JavaScript: Past, Present, and Future - NDC Porto 2020
reverentgeek
47
5k
Unsuck your backbone
ammeep
668
57k
Creating an realtime collaboration tool: Agile Flush - .NET Oxford
marcduiker
25
1.8k
Why You Should Never Use an ORM
jnunemaker
PRO
53
9k
How GitHub (no longer) Works
holman
311
140k
RailsConf & Balkan Ruby 2019: The Past, Present, and Future of Rails at GitHub
eileencodes
131
33k
What's new in Ruby 2.0
geeforr
342
31k
The Illustrated Children's Guide to Kubernetes
chrisshort
48
48k
Navigating Team Friction
lara
183
14k
Transcript
幾つかのデータ解析手法の紹介 奥 牧人 2024/08/02 富山大学ムーンショット研究数理セミナー 発表後に一部修正 1 / 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 2 /
60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 3 /
60
背景 富山大学ムーンショット研究のメンバーの多くは実験系であり、 データ解析を担当出来る 数理系の者が慢性的に不足 している。 特に、習得難易度の高いデータ解析手法の場合、手法毎に扱える 者が一人しかいない状況も珍しくない。 また、比較的簡単な解析でさえも一部の者に集中しがちである。 しかし、数理系の者を増やすのは難しい。 4
/ 60
目的 一人しか扱えない手法について、万一に備えて、せめて概要だけ でも 複数名で理解しておく ことを目指す。 比較的簡単な解析が出来る者 を増やす。 → まずはお互いに知っている手法を紹介 5
/ 60
その他の問題 (医薬系) 医薬系の場合 授業でプログラミングや機械学習をほとんど習わない。 稀にデータ解析の出来る学生や研究員が出現するが、卒業や 異動の際に知識やノウハウが誰にも継承されず失われる。 数理の課題について相談できる者が少ない。 データ解析主体では論文になりにくい。 対策 マニュアルや相談体制を整備し、研究が軌道に乗る最小人数
(critical mass) を越えることを目指す。 基本的に 実験系の者が自分でデータ解析 できるようにし、 数理の専門家は助言や指導に徹する状況を目指す。 6 / 60
その他の問題 (理工系) 理工系の場合 生命科学分野の実データが入手しづらい。 異分野のデータの場合、結果の解釈がしづらい。 標準的な手法を使っただけでは研究と認めて貰えない。 対策 公開データや人工データを使う。 良き共同研究相手 を見つける。ダメな相手は切る。
独自性の高い手法 を使いこなせるようにする。 7 / 60
以前まとめた資料 遺伝子発現量データ解析の基礎, 和漢研年報, 44:2-11 (2018). http://hdl.handle.net/10110/00019294 前処理、次元削減、発現変動遺伝子、ベン図、階層的クラス タリング、エンリッチメント解析、ネットワーク図を説明 上記解説文を元にした講義資料 https://speakerdeck.com/okumakito/oku-slide-20210721
動的ネットワークバイオマーカー解析の説明資料 https://speakerdeck.com/okumakito/oku-slide-20220820 8 / 60
何について学ぶべきか データ解析手法について学ぶとき、必ずしも 動作原理 や背後に ある 数学の概念 を理解する必要はない。 むしろ、実際のデータ解析で使用する際の 注意点 や、手法毎の
利点と欠点を踏まえた 使い分け を覚える方が大事。 自動車で例えると、エンジンの仕組みを詳細に理解していなくて も運転は出来るが、交通ルールを分かっていないと事故を起こし てしまう。 9 / 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 10 /
60
統計学 医薬系では統計学をよく使う。医薬系の者の方が理工系の者より むしろ統計学に詳しい印象がある。 特に t検定 を多用する。 色々と批判されているものの、P値 を重視する者は多い。 私はそこまで詳しくないので、本発表では深入りしない。 11
/ 60
機械学習の区分のおさらい 大区分 中区分 例 教師あり学習 分類 線形判別分析、ロジスティック 回帰、ランダムフォレスト、 LightGBM 教師あり学習
回帰 単回帰、重回帰、Lasso, ガウス 過程回帰 教師なし学習 クラスタリング K平均法、階層的クラスタリング 教師なし学習 次元削減 PCA, t-SNE, UMAP 本発表では 教師なし学習のみ 扱う。 12 / 60
機械学習の区分の説明 教師あり学習は正解データあり、未知データの予測に使う。 教師なし学習は正解データなし、既存データの解釈に使う。 教師あり学習ではデータを訓練データとテストデータに分ける。 教師なし学習では全データをそのまま使う。 分類の出力は質的変数、回帰の出力は量的変数 クラスタリングではデータを複数のグループに分ける。 次元削減ではデータ点を (通常は) 2次元平面上に並べる。
ロジスティック回帰は、回帰と付いているが分類器 t-SNEやUMAPで通常グループ毎に色が違うのは、単に後から 追加でクラスタリングをかけているだけ 13 / 60
紹介する一般的な手法 主成分分析 t-SNE UMAP 倍率変化 FDR ベン図 階層的クラスタリング K平均法 エンリッチメント解析
14 / 60
主成分分析 項目 説明 名前 主成分分析、Principal Component Analysis (PCA) 区分 教師なし学習
> 次元削減 処理内容 多次元のデータを分散の意味で最もよく見える方向から 眺めた時の射影を計算する 関連用語 寄与率、ローディング、スクリープロット 利点 簡単、色々なことが分かる、解釈しやすい 欠点 scRNA-seqなど多数のグループがある場合に向かない 注意点 3次元は静止画だと意味がない (最適でない2次元への 射影になってしまっているため)、各軸毎に反転可 15 / 60
主成分分析、補足 私はPCAをデータ全体の傾向を把握する目的で使っている。 外れ値がないか、条件別に分かれているか、時間順に並んで いるかなどを確認 あくまで傾向把握目的なので、PC1, PC2のみ使用 重みの絶対値の大きい変数にも特に着目しない。重要変数の 選択はDEGでやる。 RNA-seqの場合、低発現遺伝子は予め除外する。 中心化
(平均を0にする) は必須だが、Zスコア化 (平均を0、標準 偏差を1にする) した方が良いかどうかはデータによって違う。 16 / 60
t-SNE 項目 説明 名前 t-SNE 区分 教師なし学習 > 次元削減 処理内容
多次元のデータ点を上手に2次元上に配置する 関連用語 - 利点 多数のグループをきちんと分離出来る 欠点 本当は同じグループの点が分離する場合がある、利用率 がUMAPに押され気味 注意点 sklearnの場合、学習率をデータ点数/12にし、初期配 置をPCAで決めるようにする。また、必要に応じて事前 にPCAで10次元程度に落としてからt-SNEをかける。 17 / 60
t-SNE、補足 t-SNEや後述のUMAPでは 各軸に優劣がない。寄与率 という概念 も存在しない。 簡単に言うと、 元の高次元空間で近い点は2次元平面上でもなるべく近く 元の高次元空間で遠い点は2次元平面上でもなるべく遠く なるように、各点の2次元平面上の位置をうまく調節する手法 つまり、PCAとは根本的に考え方が違う。
18 / 60
UMAP 項目 説明 名前 UMAP 区分 教師なし学習 > 次元削減 処理内容
多次元のデータ点を上手に2次元上に配置する 関連用語 - 利点 多数のグループをきちんと分離出来る 欠点 本当は同じグループの点が分離する場合がある、パラメ ータ調節が必要 注意点 パラメータ (min_dist, spread) を適宜調節する 19 / 60
倍率変化 項目 説明 名前 倍率変化、Fold change (FC) 区分 発現変動遺伝子解析 処理内容
2群間の比較において、平均値の比が (対数を取る前の 値で) 2倍より大きい遺伝子を選ぶ 関連用語 - 利点 簡単、再現性が高い、サンプル数に依存しない 欠点 P値至上主義者に理解して貰えない 注意点 閾値は2倍から変えても構わない 20 / 60
倍率変化、補足 複数時点ある場合、私は以下のようにしている。 各時点で平均を取る。 それらの最大値と最小値の差が (log2を取った値で) 1より 大きい遺伝子を選ぶ。任意の2時点間のDEGの和集合である。 取れた遺伝子を階層的クラスタリングにかける。 クラスタ毎に遺伝子一覧を調べたりエンリッチメント解析に かけたりする。
こうする理由は、多くの場合にDEGは複数時点で大きく重複する ので、時点別で見るよりクラスタ別で見た方が解釈しやすいから である。 21 / 60
FDR 項目 説明 名前 FDR、False Discovery Rate 区分 発現変動遺伝子解析、多重検定補正 処理内容
N個のP値に対応するN個のQ値を計算し、Q値が指定し たFDRの値以下のものを有意と判定する。FDRとは出力 の中に含まれる偽陽性の割合のこと。 関連用語 Benjamini-Hochberg法 利点 Nが大きい場合にBonferroni補正より検出力が高い 欠点 取れる遺伝子数がサンプル数に依存する 注意点 P値がほぼ同じでもQ値が大きく違ったり、P値が違うの にQ値が全く同じだったりする。P値とQ値は別物。 22 / 60
ベン図 項目 説明 名前 ベン図、Venn diagram 区分 発現変動遺伝子解析 処理内容 複数の集合間の重なり具合を計算して図示する
関連用語 和集合、積集合、差集合、ジャカード指数 利点 3集合以下なら分かりやすい 欠点 5集合以上だと領域数が多過ぎて意味不明 注意点 4集合の場合は正円だと描けないので楕円等を使う 23 / 60
階層的クラスタリング 項目 説明 名前 階層的クラスタリング 区分 教師なし学習 > クラスタリング 処理内容
似ているデータ点を順にまとめていき、閾値で切って 複数のクラスタに分ける 関連用語 樹形図、平均連結法、ウォード法、シルエットスコア 利点 ランダム性がない、生命科学の分野で一般的 欠点 汎用的な閾値自動決定法が無い 注意点 RNA-seqでは相関類似度を使うか遺伝子毎にZスコア化 する、単連結法はChaining現象が発生するため非推奨 24 / 60
階層的クラスタリング、補足 設定項目 類似度 or 非類似度: ユークリッド距離、相関類似度など 連結法: ウォード法、平均連結法など クラスタ分割基準: クラスタ数指定、非類似度の閾値指定など
私はRNA-seqデータ解析の場合では相関類似度、平均連結法、 非類似度の閾値指定を使っている。閾値は以下の式で計算 意味は、相関係数をフィッシャー変換した際の3σ点 1 − tanh ( 3 √N − 3 ) 25 / 60
K平均法 項目 説明 名前 K平均法、K-means 区分 教師なし学習 > クラスタリング 処理内容
クラスタ数をKとし、各クラスタの中心位置 (最初はラ ンダム) とそれに属するデータ点を交互に更新する 関連用語 - 利点 簡単 欠点 クラスタ数を決め打ちする必要がある、クラスタの初期 配置により結果が変わる 注意点 RNA-seqでは相関類似度を使うか遺伝子毎にZスコア化 する 26 / 60
エンリッチメント解析 項目 説明 名前 エンリッチメント解析 区分 その他 処理内容 入力した遺伝子一覧の中に、どのような注釈を持つ遺伝 子が多いかをデータベースに問い合わせ、出力する
関連用語 GO (Gene Ontology)、フィッシャーの正確検定 利点 個々の遺伝子について知らなくても結果の解釈が出来る 欠点 データベースによって結果が違う、重複数が極めて少な いものも出力される、詳しいことは分からない 注意点 P値だけ見てはダメ、個数も確認する 27 / 60
エンリッチメント解析、補足 フィッシャーの正確検定とは、ベン図で2つの集合の重なりが 有意に大きいかどうかを判定するもの。 エンリッチメント解析の場合の2つの集合とは、入力した遺伝子 群 (DEGなど) と特定の注釈を持つ遺伝子群のこと。 私は普段DAVIDを使うので、以下DAVIDについて補足する。 https://david.ncifcrf.gov/ P値はDAVIDの独自補正がかかっている。私はDAVIDから結果を
ダウンロードしてPythonでP値とQ値を再計算している。その 際、遺伝子の個数が少ない注釈も除外している。 GOで名前に"regulation of transcription" と付くものは解釈が 難しいので私は除外している。 28 / 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 29 /
60
同期性揺らぎ遺伝子解析 (DNB解析) 項目 説明 名前 同期性揺らぎ遺伝子解析、DNB解析 区分 その他 処理内容 揺らぎが大きく、かつ、互いに強く相関した遺伝子群を
選ぶ 関連用語 中央絶対偏差、スピアマンの相関係数、階層的クラスタ リング 利点 揺らぎに注目する点が独自的 欠点 言葉が独り歩きしていて過度に期待される、結果が合っ ているかどうかよく分からない、手法が確立していない 注意点 後述 30 / 60
主な流れ 主な流れ 1. 揺らぎの大きな遺伝子を選ぶ 2. その中から相関強度の高いクラスタを選ぶ 必要に応じて前処理や後処理をする。 本質的に 変数選択 の手法、発現変動遺伝子と同様
31 / 60
考えるべきこと 考えるべきことが沢山ある。 揺らぎの指標: 標準偏差、中央絶対偏差 揺らぎ判定法: 倍率変化、割合、仮説検定 倍率変化や割合の閾値 仮説検定の種類: F検定、Levene ル
ビ ー ン 検定、Brown-Forsythe検定 対照群: 使用する、使用しない 複数時点の扱い: 時点毎、中心化して結合、対照群のみプール 相関の指標: ピアソンの相関係数、スピアマンの相関係数 連結法: 平均連結法、完全連結法 クラスタ分割基準: 個数、閾値 32 / 60
考えるべきこと、続き クラスタ選択: 最大クラスタのみ選択、2番目以降も選択 最大クラスタに対する大きさ 最大クラスタ数 最小クラスタサイズ 事後クラスタリング: する、しない 時点毎の中心化: する、しない
相関係数の指標: ピアソン、スピアマン 最大クラスタ数 最小クラスタサイズ 遺伝子全体の揺らぎスコア: 選択条件に使う、使わない 33 / 60
RNA-seqに対する現在の方針 外れ値に強い 中央絶対偏差 と スピアマンの相関係数 は必須 結果の再現性を高めるため 対照群なし を選択 時点毎に計算するか、中心化して結合するかは、結果があまり
変わらないことが多いのでどちらでも良い。 対照群なし、かつ、時点毎の場合、異なる時点間で取れる遺伝子 が大きく重複するので、事後クラスタリング が必要 結果の妥当性を以下で確認 ヒートマップ 遺伝子リスト エンリッチメント解析 34 / 60
何が難しいのか? 手法が確立しておらず、選択肢の組み合わせの数が膨大なため、 データ毎に 数理の知識 と データ解析の経験 に基づき組み合わせ を変え、その妥当性を高度な ドメイン知識 に基づき多角的かつ
慎重に判断する必要がある。 私自身もマイクロアレイとRNA-seqのドメイン知識しか持ってい ないので、他の種類のデータには全く手が出せない。 今回の目的はお互いに知っている手法を教え合うことだが、正直に言 うと DNB解析の習得は勧められない。やるとしても、再実験が容易 で結果の再現性を確認しやすいものが良いと思う。 35 / 60
よくある間違い 1サンプルのみで高発現し、残りのサンプルでは全然揺らいでい なかった → 外れ値 に引きずられただけ サンプルを番号順にした場合、あるサンプルまでは一様に低く、 それ以降が一様に高かった → バッチ効果
を拾っただけ snoRNA (Snorから始まる) など短い転写産物ばかり → S/N比 の悪い遺伝子が取れただけ Olfrなど平均発現量の低い遺伝子ばかり → 同上 ミトコンドリア遺伝子 (mtから始まる) やリボソーム蛋白質 (Rp から始まる) ばかり → 単に 高発現 遺伝子を拾っただけ ほぼ全ての条件で同じ遺伝子が出てきた → 単に コントロール群 で偶然値が揃っていた 遺伝子を拾っただけ 36 / 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 37 /
60
エネルギー地形解析 項目 説明 名前 エネルギー地形解析、Energy Landscape Analysis (ELA) 区分 その他
処理内容 データを二値化し、イジングモデルに当てはめ、エネル ギーの谷や、谷と谷の間のエネルギー障壁を計算する 関連用語 イジングモデル、ベースングラフ、非連結性グラフ 利点 独自性が高い、解析ツールの使用自体は簡単 欠点 約7個の変数を選ぶ必要がある、結果の解釈が難しい 注意点 後述 38 / 60
エネルギー地形解析、補足 江崎先生がMatlabで作ったツールのPython移植版: https://github.com/okumakito/elapy エネルギーは仮想的なもので、確率と関連した値 単独のエネルギーの値に意味はなく、相対的な差が大事 頂点番号は、0と1から成る2進数を10進数に直したもの 例、1100111 = 64 +
32 + 4 + 2 + 1 = 103 ベースングラフでは、頂点が多いほど谷が広く、最小エネルギー が小さいほど谷が深いことを意味する。 非連結性グラフは、状態間 (正確にはエネルギー極小パターン間) を遷移する際のエネルギー障壁の高さを図示したもの。 p ∝ exp(−E) 2 39 / 60
エネルギー地形解析、補足続き 混乱を避けるため、ベースングラフの連結成分を 状態、個々の 頂点を パターン と呼び分ける。 ベースングラフの矢印は、1変数のみが異なる近傍パターンの中 でエネルギーが最小のものに向かって引いたもの、つまり、最急 降下方向を表す。実際は逆方向に遷移することもある。 非連結性グラフを元にした等高線図の作成も出来る。
40 / 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 41 /
60
お断り ここから先はPythonのコードを載せる。 RNA-seqに興味の無い方にも、データフレームの使い方の参考に なるとありがたい。 データフレームとは、各行と各列の名前付きの2次元データ 42 / 60
主な前処理 FASTQ → BAM → カウントデータ (私は出来ないので外注) ファイルの結合 不要な行や列の除外 ID変換
列名の変更 正規化 対数化 外れ値の除外 低発現遺伝子の除外 その他、対象データ特有の問題に対する処理 43 / 60
ファイルの結合 データファイルはなるべく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
ファイル読み込みの補足 データがコンマ区切りでなくタブ区切りの場合 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
不要な行や列の除外 ここは千差万別。比較的よく使うもの ある列 (仮に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
不要な行や列の除外、続き ある列の値がある文字列から始まる行のみ抽出、他 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
ID変換 Ensembl ア ン サ ン ブ ル IDなどの分かりにくいIDを遺伝子記号に変換 例、ENSMUSG00000024401
→ Tnf 複数の元IDが同じ遺伝子記号に対応する場合は最大値を採用 対応する遺伝子記号の無い転写産物のデータは削除 やり方は色々ある。例えば、遺伝子記号の列をindexにして groupbyするなど df = df.set_index('gene_sym').groupby(axis=0,level=0).max() 48 / 60
列名の変更 必要に応じて自分が分かりやすいように各列の名前を変更する。 手動で元ファイルを変更するのでなく、スクリプトを使う。 ミスを防ぐため、後で結果を再現できるようにするため データが数列程度なら直接書いても可 df.columns = ['ABC','DEF','GHI'] # 3列の列名を更新
長い場合は、for文などで修正名リストを作ってから更新 データの最初だけ見て番号順になっていると過信しない。 正規表現を使わずとも、大抵のことはreplaceとsplitで可能 49 / 60
正規化 マイクロアレイでは分位数正規化 (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
対数化 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
外れ値の除外 PCAとヒートマップを目視で確認し、外れ値のサンプルを除外 PCAでは同じ条件の他の点から遠く離れた位置に来る。 ヒートマップでは縦線が出る。 52 / 60
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
ヒートマップの使い方 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
ヒートマップの使い方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
低発現遺伝子の除外 低発現遺伝子の除外は必要 発現量が小さいほどS/N比が悪い。 0と非ゼロの差が対数化で増幅される。 平均値や中央値が閾値以上の遺伝子のみ抽出 df = df[df.mean(axis=1)>=5] # 平均が5以上の遺伝子を抽出
多臓器データの場合、臓器毎に計算する。 消化酵素など一部の臓器のみで発現する遺伝子がいるため 56 / 60
前処理が妥当かどうかの確認 私は以下の4つで前処理が妥当かどうかを確認している。 箱ひげ図 PCA DEGのヒートマップ DEGのエンリッチメント解析 大事なのは 違和感に気付く こと。 57
/ 60
Outline はじめに 一般的な手法 同期性揺らぎ遺伝子解析 (DNB解析) エネルギー地形解析 RNA-seqデータ前処理 まとめ 58 /
60
まとめ、その1 私の知っている複数のデータ解析手法について紹介した。 一般的な手法として、PCA、階層的クラスタリング、倍率変化に 基づくDEG取得、エンリッチメント解析などを説明した。このあ たりは実験系の方が自分で扱えるようになることが望ましい。 DNB解析が難しい理由として、手法が確立しておらず、選択肢の 組み合わせが膨大で、データ毎に結果が妥当かどうかを判断する 際に高度なドメイン知識が必要なためと説明した。 59 /
60
まとめ、その2 エネルギー地形解析を使うこと自体は簡単だが、解釈が難しいこ とを説明した。逆に言うと、慣れれば武器になる。 RNA-seqデータの前処理をPythonコードと共に説明した。色々 と書いたものの、結局のところ「これをやっておけば大丈夫」と いう保証は無いので、きちんと結果を確認して違和感を見落とさ ないことが何より大事。 60 / 60