講義資料PDF - アグリバイオインフォマティクス教育研究ユニット

USBメモリ中のhogeフォルダをデス
クトップにコピーしておいてください。
前回（5/19）のhogeフォルダが
デスクトップに残っているかも
しれないのでご注意ください。
機能ゲノム学
第3回
大学院農学生命科学研究科
アグリバイオインフォマティクス教育研究プログラム
門田幸二（かどたこうじ）
[email protected]
http://www.iu.a.u-tokyo.ac.jp/~kadota/
May 26 2015
1
講義予定

第1回（2015年5月12日）





遺伝子発現行列作成（データ正規化）
クラスタリング（データ変換や距離の定義など）
教科書の3.2節周辺
第3回（2015年5月26日）



原理、各種データベース、生データ取得
教科書の1.2節、2.2節周辺
第2回（2015年5月19日）


細胞中で発現している全転写物（トランスクリプトーム）
の解析技術は、マイクロアレイから次世代シーケンサ（
RNA-seq）に移行しつつあります。しかしRNA-seqデー
タ解析の多くは、マイクロアレイの知識を前提としていま
す。本科目では、マイクロアレイデータを主な例として、
各種トランスクリプトーム解析手法について解説します。
実験デザイン、発現変動解析（多重比較問題）、M-A plot
教科書の3.2節と4.2節周辺
第４回（2015年6月9日）

機能解析（Gene Ontology解析やパスウェイ解析）、分類など
May 26 2015
教科書
2
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
3
教科書p107-108
実験デザイン（§3.2.2）

「個体間のばらつき」を考慮した実
験デザインを組むべし、という話
Affymetrix GeneChip

Ge et al., Genomics, 86: 127-141, 2005



GSE2361、GPL96 (Affymetrix Human Genome U133A Array)、22,283 probesets
ヒト36サンプル：Heart (心臓）、Thymus (胸腺)、Spleen (脾臓)、Ovary (卵巣)、Kidney (腎
臓)、Skeletal Muscle (骨格筋）、Pancreas (膵臓）、Prostate (前立腺)、…
Nakai et al., Biosci Biotechnol Biochem., 72: 139-148, 2008 8匹のラットを使用
GSE7623、 GPL1355 (Affymetrix Rat Genome 230 2.0 Array)、31,099 probesets
 ラット24サンプル：Brown adipose tissue (褐色脂肪組織; BAT）8サンプル、White adipose
tissue (白色脂肪組織; WAT）8サンプル、 Liver (肝臓; LIV）8サンプル
 BAT 8サンプル：通常（BAT_fed） 4サンプル対 24時間絶食（BAT_fas） 4サンプル
 WAT 8サンプル：通常（WAT_fed） 4サンプル対 24時間絶食（WAT_fas） 4サンプル
 LIV 8サンプル：通常（LIV_fed） 4サンプル対 24時間絶食（LIV_fas） 4サンプル
 Kamei et al., PLoS One, 8: e65732, 2013 10匹のラットを使用
 GSE30533、 GPL1355 (Affymetrix Rat Genome 230 2.0 Array)、31,099 probesets
 ラット10サンプル：全てLiver (肝臓）サンプル
 iron-deficient diet （Iron_def） 5サンプル対 control diet （Control） 5サンプル

May 26 2015
4
2群間比較が主な目的であり、各群につき
5反復(five replicates)とっている。生物学
的なばらつき(biological variation)を考慮す
 Kamei et al., PLoS One, 8: e65732, 2013 べく、反復データは別々の個体からとって
いる(biological replicates)
 GSE30533、 GPL1355 (Affymetrix Rat Genome 230 2.0 Array)、31,099 probesets
 ラット10サンプル：全てLiver (肝臓）サンプル
 iron-deficient diet （Iron_def） 5サンプル対 control diet （Control） 5サンプル
教科書p107-108
実験デザイン（§3.2.2）
Fe
May 26 2015
Fe
Fe
Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
5
このやり方で得られる結論は
限定的!できるだけ多様な別個
体サンプルを沢山用いるべし!
教科書p107-108
実験デザイン（§3.2.2）

Kamei et al., PLoS One, 8: e65732, 2013
対比的な用語は技術的なばらつき（technical
GSE30533、 GPL1355
(Affymetrix Rat Genome 230 2.0 Array)、31,099 probesets
variation）であり、同一個体由来サンプルを分割
 ラット10サンプル：全てLiver (肝臓）サンプル
して得られた反復データ（technical replicates）
 iron-deficient diet （Iron_def） 5サンプル対 control diet （Control） 5サンプル

ラットB
ラットA
Fe
Fe
Fe
Fe
Fe
May 26 2015
Fe
Fe
Fe
Fe
Fe
6
2群間での発現変動遺伝子(DEG)
検出結果は多くなる傾向。多けれ
ばいいというものではない!
教科書p107-108
実験デザイン（§3.2.2）

Technical replicatesだと…
1. 自分は「鉄欠乏対通常」の違いを見ているつもりでも、個体間の他の違い（身長、体重な
ど）由来要因との区別がつかない
高身長対低身長、低体重対高体重、他の病気の有無、家系の違いなど
2. 得られる結果から導き出される結論は、そのラット間のみで成立する事象であり、ラットという
生物種全体に適用可能なわけではない
ラットB
ラットA
Fe
Fe
Fe
Fe
Fe
May 26 2015
Fe
Fe
Fe
Fe
Fe
7
教科書p107-108
実験デザイン（§3.2.2）

普遍的な結果を得たいのなら、できるだけ
多様な別個体サンプルを沢山用いるべし!
Expression Atlasも3 biological replicates
以上を基本としているようだ。
Biological replicatesでも多様性が不十分な場合はイマイチ…
家系A
Fe
Fe
家系B
Fe
Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Fe Fe Fe
Fe Fe
Fe
Technical repliatesよりはましだろうが、鉄
欠乏状態による結果なのか家系の違いに
よるものなのかの区別はつけられない。
May 26 2015
8
クラスタリングと発現変動解析
May 26 2015
クラスタリング結果を眺めれ
ば、発現変動遺伝子（DEG）
数に関するおおよその見当
がつきます。
→ クラスタリングって重要
9
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）




コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
10
対数変換後のデータを用いて2群間比較
データ解析もいろいろ
発現変動遺伝子同定
クラスタリング
遺伝子発現行列
機能解析
・Gene Ontology（GO）
・パスウェイ解析
分類（診断）
遺伝子ネットワーク推定
May 26 2015
11
比較したいグループ間で発現変動
している遺伝子または転写物を同
定することはデータ解析の基本
2群間比較



農産物の栽培条件の違い（通常対低温、通常対乾燥）
味の違い（おいしい対まずい）
サンプルの状態の違い（癌対正常）
G1群 G2群
高発現
G1群
N genes
低発現
May 26 2015
発現変動
度合いで
ソート
G2群
G1群
G2群
2群間で発現の異なる遺伝子
(Differentially Expressed Genes;
DEGs)を抽出
12
2群間比較

相関係数rの絶対値が大きいほど
発現変動の度合いが大きいと解釈
パターンマッチング法

理想的なパターンyとの類似度が高い順にランキング
1 n
 ( xi  x)( yi  y)
n  1 i 1
相関係数 r 
(1  r  1)
n
n
1
1
2
2
(
x

x
)
(
y

y
)
 i
 i
n  1 i 1
n  1 i 1
May 26 2015
13
パターンマッチング法
May 26 2015
テンプレートパターン情報を含
むファイルを読み込んでパタ
ーンマッチングを行ってみよう
14
パターンマッチング法
May 26 2015
入出力の全体像。出力ファイルは、
入力ファイルin_f1の右側に相関係数
計算結果が追加されたもの。
15
コードの解説
May 26 2015
入力ファイルの読み込みがう
まくできていることがわかる
16
コードの解説
May 26 2015
読み込み時にheader=TRUEや
row.names=1の記述がない点に注目!
17
コードの解説
May 26 2015
hogeの中身は入力ファイルと同じだが、
欲しいのはhogeオブジェクトの2列目部分
18
読み込み時にrow.names=1をつけてもよい
Tips
May 26 2015
19
apply関数は行ごとや列ごとに同じ関数
を繰り返して実行させたい場合に便利。
①dataオブジェクトの、②各行に対して
、③cor関数を適用せよ。その際、④テ
ンプレートyはdata.clとし、⑤相関係数
の種類はparamで指定したものとする。
コードの解説：apply
①
May 26 2015
② ③
④
⑤
20
Tips：cor, データの型
May 26 2015
apply関数を使わずに、遺伝子（つまり行
）ごとにcor関数を用いて相関係数を計
算する手順。as.numeric関数は、data.cl
オブジェクトとデータの型を揃える目的
で利用している。「互換性のない次元で
す」と言われているが、ここでは、相関係
数を計算するときの入力データの見栄
えが異なると文句を言われていると解釈
すればよい。この「見栄え」が「データの
型」と呼ばれるものに相当する。
21
Tips：as.vector, mode
May 26 2015
慣れてくるとas.numericを使えばいいと
すぐにわかるが、はじめのうちは
as.character, as.integer, as.vectorなど
既知のas…関数候補の中からそれっぽ
いものを試して型（見栄え）が揃うもの
を探すのが一般的かも…。modeという
データの型を表示する関数もある。
22
コードの解説：cbind
May 26 2015
cbind関数を用いて出力させたい順番で列
(column)方向で結合(bind)したものがtmp
オブジェクト。ちなみに行(row)方向で結合
させたい場合は、rbind関数を用いる。
23
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
24
教科書p108-116
t-test：仮想データ
このウェブページを用いたDEG検出手順
の一般的なグループごとの反復数指定方
法です。G1群は6反復、G2群は5反復。同
一群でまとまっている、という前提です。
①
②
May 26 2015
25
教科書p108-116
t-test：仮想データ
May 26 2015
このウェブページを用いたDEG検出手順
の一般的なグループごとの反復数指定方
法です。G1群は6反復、G2群は5反復。同
一群でまとまっている、という前提です。
gene1のような比較するグループ間（G1群
対 G2群）で明らかに発現の異なる遺伝子
（DEG）のp値は0に近い値となり、明らか
にDEGではないもの（non-DEG）のp値は1
に近い値になるという感覚は重要。
26
教科書p108-116
t-test：乱数
May 26 2015
N = 10,000遺伝子（行）からなる遺伝
子発現行列（各群3サンプル）を入力
として、遺伝子ごとにt-testを実行し、
p < 0.05を満たす遺伝子数を眺めるこ
とを通じて多重比較問題を実感する
27
教科書p108-116
t-test：乱数
May 26 2015
data.clオブジェクトは、テンプレートパ
ターンのようなもの。入力データに相
当するdataオブジェクトの1-3列がG1
群、4-6列がG2群由来サンプルだと
いうことを指し示すクラスラベル情報
28
教科書p108-116
確かに（標準）正規分布乱数になっている
t-test：乱数
May 26 2015
29
教科書p108-116
t-test：乱数
May 26 2015
特定の条件を満たす
列のみ取り出すやり方
30
教科書p108-116
t-test：乱数
May 26 2015
有意水準α = 0.05とすると、1行
目の遺伝子のp値は0.05未満で
はない。2行目、3行目、…。
31
教科書p108-116
Tips
May 26 2015
str関数を利用すれば、t.test関数に限ら
ず、出力結果からどのような情報を取り出
せるかを知ることができます。
32
教科書p108-116
t-test：乱数
May 26 2015
出力ファイル(hoge3.txt)は、入力ファイル
情報の右側にp.value, q.value, rankingとい
う列の情報が追加されたものです。
33
p.value列かranking列で昇順にソ
ートすると、発現変動順になる。
教科書p108-116
t-test：乱数
hoge3.txt
May 26 2015
34
教科書p108-116
EXCELでのソートの仕方
Tips
①
②
④
③
May 26 2015
35
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
36
教科書p108-116
多重比較問題のイントロ
DEGの存在しないnon-DEGの
みからなるデータなので妥当
①
p < 0.05を満たす
遺伝子数は492個
②
p < 0.10を満たす
遺伝子数は959個
May 26 2015
37
教科書p108-116 参考
Type-I error （false positive）
ランダムデータの場合


有意水準αでN回の検定（多重比較）を行うと、
（N×α）個のFalse Positiveが得られる。
10000個の遺伝子（N=10000）に対してp < 0.05を満
たすものを調べる（有意水準αを0.05に設定すること
と同義）と（N×α）個程度が本当は発現変動遺伝子
(Differentially Expressed Genes; DEGs)でないにも
かかわらず発現変動遺伝子と判断されてしまう。
May 26 2015
38
ここではq-valueとしているが
、adjusted p-valueと呼ばれる
ものと同じです。Rパッケージ
出力結果でもadjPやPadjや
 うれしくない結果：「実際に得られた発現変動遺伝子数 ≒
FDRなどの列がありますが、
（解析遺伝子数N×設定した有意水準α）個」
それに相当する議論です。
 このデータ中には「発現変動遺伝子 (DEG)はない」と判断する。
教科書p108-116 参考
p値だけである程度判断できるが…

うれしい結果：「実際に得られた発現変動遺伝子数 >> （解析
遺伝子数N×設定した有意水準α）個」


このデータ中には「真の発現変動遺伝子が存在する」ことが期
待される。
実際に利用されているRパッケージの多くは、(多重比較を考
慮した補正後のp-valueに相当する)q-valueの値を出力する
 (p値利用時の有意水準αに相当する) False Discovery
Rate (FDR)の閾値を満たす遺伝子数を頼りに発現変動
遺伝子の有無を判断する
May 26 2015
39
参考
Benjamini and Hochberg J. Roy. Stat. Soc. B, 57: 289-300, 1995.
多重比較問題：FDRって何?

p-value (false positive rate; FPR)



本当はDEGではないにもかかわらずDEGと判定してしまう確率
全遺伝子に占めるnon-DEGの割合（分母は遺伝子総数）
例：10,000個のnon-DEGからなる遺伝子をp-value < 0.05で検定すると、
10,000×0.05 = 500個程度のnon-DEGを間違ってDEGと判定することに相当




実際のDEG検出結果が900個だった場合：500個は偽物で400個は本物と判断
実際のDEG検出結果が510個だった場合：500個は偽物で10個は本物と判断
実際のDEG検出結果が500個以下の場合：全て偽物と判断
q-value （false discovery rate: FDR）



DEGと判定した中に含まれるnon-DEGの割合
DEG中に占めるnon-DEGの割合（分母はDEGと判定された数）
non-DEGの期待値を計算できれば、p値でも上位x個でもDEGと判定する手段は
なんでもよい。以下は10,000遺伝子の検定結果でのFDR計算例



May 26 2015
p < 0.001を満たすDEG数が100個の場合：FDR = 10,000×0.001/100 = 0.1
p < 0.01を満たすDEG数が400個の場合：FDR = 10,000×0.01/400 = 0.25
p < 0.05を満たすDEG数が926個の場合：FDR = 10,000×0.05/926 = 0.54
40
Benjamini and Hochberg J. Roy. Stat. Soc. B, 57: 289-300, 1995.
多重比較問題：FDRって何?

DEGかnon-DEGかを判定する閾値を決める問題



DEG数に関するよりよい結果
を得たい場合には、p-value
ではなくq-valueを利用しまし
ょう。（閾値を有意水準αでは
なくFDRで設定しましょう。）
有意水準5%というのがp-value < 0.05に相当
False discovery rate (FDR) 5%というのがq-value < 0.05に相当
発現変動ランキング結果は不変なので上位x個という決め打ちの場合
にはこの問題とは無関係
May 26 2015
41
教科書p108-116
多重比較問題
q < 0.05を満たす遺伝子数は0
個。DEGの存在しないnon-DEG
のみからなるデータなので妥当
①
p < 0.05を満たす
遺伝子数は492個
②
p < 0.10を満たす
遺伝子数は959個
May 26 2015
42
FDR
教科書p108-116 expected = 全遺伝子数×p-value =
expected = 10,000×0.0002304 = 2.304。
FDR = 偽物検出割合 =
expected/observed = 2.304 / 2 = 1.1518。
2は、p = 0.0002304以下をDEGとみなした
場合に2個がDEGということ。
基本的にこの2つは同じものという理
解でよい。より正確には、FDR列の情
報をもとに値の分布が滑らかになるほ
うに細工しているのがq.value列の数値
May 26 2015
43
FDR
May 26 2015
教科書p108-116 参考 p.valueが高いもののFDR値から順に
見ていって、FDR値が低くなるように
置換していったものがq.value列の値
44
教科書p115
参考
自力でq-value (FDR)計算
May 26 2015
・FDR = 偽物検出割合
・FDR = expected/observed
45
教科書p115-116 参考
自力でq-value (FDR)計算
May 26 2015
p値計算結果が手元にあれば（
つまりp.valueオブジェクトがあれ
ば）このコードを実行することに
よってFDRの概要がわかります
46
教科書p115-116 参考
自力でq-value (FDR)計算
May 26 2015
ここで指定しているのはp-value
の閾値（つまり有意水準）です
47
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
48
教科書p108-116
t-test：DEGを含むデータ
May 26 2015
10,000個中1,000個がG1群で高
発現の2群間比較用シミュレーシ
ョンデータ。確かにG1群で高発現
になっていることがわかります
49
教科書p108-116
t-test：DEGを含むデータ
①p < 0.05を満たす遺伝子数は
1,226個。期待値は500個なので、
(1,226 – 500)個程度が本物だと
判断する。②q < 0.20を満たす遺
伝子数は388個。FDR = 0.20なの
で、388*0.2 = 77.6個は偽物で残
りの80%は本物だと判断する。
①
②
May 26 2015
50
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
51
教科書p167- limmaという発現変動解析用Rパッケ
ージを用いてDEG検出を行います。
発現変動解析(limma)
May 26 2015
52
教科書p167-
発現変動解析(limma)

GSE7623データを用い、様々な2群
間比較を行い、クラスタリング結果
とDEG検出結果の関係をみてみよう
Nakai et al., Biosci Biotechnol Biochem., 72: 139-148, 2008


GSE7623、 GPL1355 (Affymetrix Rat Genome 230 2.0 Array)、31,099 probesets
ラット24サンプル：Brown adipose tissue (褐色脂肪組織; BAT）8サンプル、White
adipose tissue (白色脂肪組織; WAT）8サンプル、 Liver (肝臓; LIV）8サンプル



BAT 8サンプル：通常（BAT_fed） 4サンプル対 24時間絶食（BAT_fas） 4サンプル
WAT 8サンプル：通常（WAT_fed） 4サンプル対 24時間絶食（WAT_fas） 4サンプル
LIV 8サンプル：通常（LIV_fed） 4サンプル対 24時間絶食（LIV_fas） 4サンプル
①
②
rcode_clustering_png.txtの実行結果。
①肝臓と脂肪間で大きく二つのクラス
ターに分かれている。
②脂肪の中でも白色脂肪と褐色脂肪
に分かれている。
③褐色脂肪は空腹（24時間絶食）と
満腹（通常）できれいに分かれている。
③
May 26 2015
53
教科書p167-
発現変動解析(limma)
May 26 2015
解析1の予想：DEGなし。解析2~4
の予想：DEGあり。予想される
DEG数：解析2 < 解析3 < 解析4
54
教科書p167-
発現変動解析(limma)
May 26 2015
rcode_limma_basic.txt。解析1用。テ
ンプレートからの変更点および追加
点。Macのヒトは区切り文字に注意
55
rcode_limma_basic.txt
入力ファイル（data_mas_EN.txt）
読み込み後のdataオブジェクト
は24サンプルからなる
¥¥
¥¥
May 26 2015
56
rcode_limma_basic.txt
May 26 2015
posiで指定した列番号の
みからなるサブセットを抽
出できていることがわかる
57
教科書p167-
発現変動解析(limma)
data.clで指定している情報は
、グループラベル情報（どの
列がどの群由来かということ）
解析1の予想：DEGなし
May 26 2015
58
教科書p167- 通常（私）は、0.05～0.30あたりの
FDR閾値を調査→DEGがないと判断
発現変動解析(limma)
q < 0.70を満たす遺伝子数は330個。
解析1の予想：DEGなし
FDR = 0.70なので、330*0.7 = 231個は
偽物で残りの30% (つまり330*0.3 = 99
個)は本物だと判断することになる…
q < 0.65を満たす遺伝子数は167個。
FDR = 0.65なので、167*0.65 = 108.55
個は偽物で残りの35% (つまり167*0.35
= 58.45個)は本物だと判断する…
May 26 2015
59
rcode_limma_basic.txt
rcode_limma_all.txt（の一部）
May 26 2015
rcode_limma_basic.txtで動作確認を
してから、param_posiのように変更予
定箇所を上のほうに移動して、解析24用のコードを作成する（のが門田流）
60
教科書p167-
解析2用のコード
rcode_limma_all.txt（の一部）
発現変動解析(limma)
解析2の予想：DEGあり
May 26 2015
61
教科書p167- 通常（私）は、0.05～0.30あたりの
FDR閾値を調査→DEGがあると判断
発現変動解析(limma)
q < 0.1を満たす遺伝子数は38個。FDR
解析2の予想：DEGあり
May 26 2015
= 0.1なので、38*0.1 = 3.8個は偽物で
残りの90% (つまり38*0.9 = 34.2個)は
本物だと判断することになる…
62
教科書p167- 通常（私）は、0.05～0.30あたりの
FDR閾値を調査→DEGがあると判断
発現変動解析(limma)
解析3の予想：DEGあり
May 26 2015
63
教科書p167- 通常（私）は、0.05～0.30あたりの
FDR閾値を調査→DEGがあると判断
発現変動解析(limma)
解析4の予想：DEGあり
May 26 2015
64
教科書p167-
DEG検出結果まとめ
G2群
G2群
G2群
解析1の予想：DEGなし。解析2~4
の予想：DEGあり。予想される
DEG数：解析2 < 解析3 < 解析4
G2群
G1群
May 26 2015
65
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
66
参考
4.と5.と6.は実質的に同じです
Tips：関数の利用法
May 26 2015
67
参考
Tips：関数の利用法
May 26 2015
ネット接続環境であれば、ここで
提供している関数を利用可能
68
参考
Tips：関数の利用法
May 26 2015
ネット接続環境でなくても、一旦
R_functions.Rファイルを作業ディ
レクトリにダウンロードしておけば
Students_ttest関数を利用可能
69
参考 limma以外にも様々なパッケージがありま
多数の方法があります
May 26 2015
す。ウェブページでリストアップされていな
い方法(例：右下のFCROS)も多数あり。こ
こではlimmaと似た系列のIBMT法を紹介。
70
参考
多数の方法があります
IBMT法は、limmaパッケージ中の関数
を内部的に用いています。limmaを基
本としつつ、改良を加えた関数部分の
み提供しているという解釈でもよい。
IBMT法の実体であ
るIBMT関数の読み
込みを行っている
May 26 2015
71
参考
発現変動解析(IBMT)
IBMT法はlimma系なので、全体的な
傾向は変わりません。そして違いの
程度を把握してもらうのが目的です。
解析1の予想：DEGなし
解析2~4の予想：DEGあり
予想されるDEG数：解析2 < 解析3 < 解析4
May 26 2015
72
参考
発現変動解析(IBMT)
May 26 2015
rcode_ibmt_basic.txt。解析1用。テン
プレートからの変更点および追加点
。Macのヒトは区切り文字に注意
73
解析結果(IBMT)
May 26 2015
rcode_ibmt_all.txtの結果。同じDEG検出法でも入
力データ（前処理法：MAS5, RMA, RobLoxBioC）が
異なると結果も異なることが分かる。
74
課題：limmaとIBMTの解析結果を考察せよ
課題：解析結果(limma)
May 26 2015
75
Contents


実験デザイン（教科書の§3.2.2）
2群間比較：発現変動遺伝子（DEG）検出

パターンマッチング法（相関係数の利用）


多重比較問題とFalse Discovery Rate (FDR)



正規分布乱数由来のDEGが存在しないデータでStudent’s t-test
10% DEGが存在する正規乱数でデータ（10,000個中1,000個がDEG）でStudent’s t-test
発現変動解析用Rパッケージの利用（§4.2.1, p167-）





コードの中身をおさらい、apply関数の基本的な利用法など
limmaパッケージ (Smyth GK, SAGMB, 2004)
関数の利用法
IBMT法 (Sartor et al., BMC Bioinformatics, 2006)
課題
描画（M-A plot）


May 26 2015
作成法
同一群内のばらつき（前処理法間の違い）
76
M-A plot
May 26 2015
limma解析結果（解析4のFDRが0.05を満た
す2,892 probesets）のM-A plotを描画しよう
77
rcode_limma_MAplot_basic.txtは、
②の例題7をテンプレートにしてい
ます。最終形は右下の図ですが、
順を追って説明していきます。
M-A plot
①
②
May 26 2015
78
M-A plot
May 26 2015
M-A plotの説明。横軸のAは平均発現
レベル、縦軸のMはlog2(G2/G1)に相当
79
mean_G1とmean_G2は、単にグルー
プごとの平均値を算出しているだけ
横軸のAは平均発現レベル、
縦軸のMはlog2(G2/G1)に相当
May 26 2015
80
①(mean_G2 – mean_G1)の計算結果を
縦軸のMとして計算できるのは、発現
レベルが対数変換後のデータだから
①
横軸のAは平均発現レベル、
縦軸のMはlog2(G2/G1)に相当
May 26 2015
81
plotの基本形
May 26 2015
82
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
May 26 2015
黒丸の塗りつぶしにすべく、
pch=20オプションを追加
83
プロットの大きさをデフォルトの10%
にすべく、cex=.1オプションを追加
May 26 2015
84
グリッド線を追加
May 26 2015
85
指定したFDR閾値を満た
すDEGの位置情報を取得
May 26 2015
86
objベクトルがTRUEの場所を
magenta色で描画。cexとpchオプ
ションの値を同じにすることで色
だけを変更していることに相当
May 26 2015
87
参考
1
May 26 2015
上位x個のみ色を変えることも簡単にで
きます。rcode_limma_MAplot_basic.txtの
下のほうのコードです。これは、Top10
のみマゼンタ色にするやり方です。
88
参考
2
May 26 2015
上位x個のみ色を変えて大きくすることも
できます。rcode_limma_MAplot_basic.txt
の下のほうのコードです。cexのところの
値を大きくして、通常の1.5倍の大きさに
しています。
89
参考
May 26 2015
表示範囲を自在に変更可能です。
rcode_limma_MAplot_basic.txtの下
のほうのコードです。
90
param_posiをc(5,6,7,8)に変更すれ
ばBAT_fas内（同一群内）のばらつ
きの程度を調べることに相当。
rcode_limma_MAplot_basic2.txtです
May 26 2015
91
rcode_limma_MAplot_basic2.txt
の下の方のコードです。FoldchangeによるDEG検出の危険
性がよくわかります。
May 26 2015
92
同一群内のばらつきを概観
IBMT法実行結果。同一群内の
ばらつきの程度は、前処理法内
では概ね同じだが前処理法間で
は大きく異なる。→前処理法の
選択やDEG検出法との相性あり
MAS5データ
RMAデータ
May 26 2015
93

Download Report