機能ゲノム学（第6回）

Size: px

Start display at page:

Download "機能ゲノム学（第6回）"

なおみありはら
6 years ago
Views:

1 RNA-Seq データ解析リテラシー東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp 1

2009 年ごろの私次世代シーケンサー (NGS) 解析についての認識単に短い塩基配列が沢山あるだけでしょ得られる配列データって multi-fasta 形式のもので単にそれをリファレンス配列にマッピングしてカウントするだけでしょそれ以降の解析はマイクロアレイと同じなんじゃないのー私について

2 2009 年ごろの私次世代シーケンサー (NGS) 解析についての認識単に短い塩基配列が沢山あるだけでしょ得られる配列データって multi-fasta 形式のもので単にそれをリファレンス配列にマッピングしてカウントするだけでしょそれ以降の解析はマイクロアレイと同じなんじゃないのー私についてマイクロアレイを中心としたデータ解析手法の開発主に遺伝子発現行列の数値データのみを取り扱ってきた配列解析系のスキルはほぼゼロで用語がまるでわかっていないアグリバイオインフォマティクス教育研究プログラムの活動の一環で smallrna の NGS(Illumina) 解析をやりはじめた自分の研究テーマとして主体的にやり始めたのは 2011 年 ~ 2

3 取り扱う RNA-Seq データの基本形式データ取得完了! なんじゃこの変な記号は! 何をどうすれば... FASTQ 形式? 3

4 Contents RNA-Seq データ取得 ( マップする側 ) 基本形式 (FASTQ 形式 ) 公共データベースから取得する場合クオリティのカットオフマッピングに使うリファレンス配列 ( マップされる側 ) ゲノム配列 (RefSeq のような ) トランスクリプトーム配列リード数のカウントやデータの正規化 (RPKM) 分布の話 ( ポアソン分布負の二項分布 ) Rを使って二群間比較 ( 発現変動遺伝子検出の流れ ) なぜ倍率変化 ( 何倍発現が変化したかでの議論 ) がだめなのか ( 自分のデータ解析の際に路頭に迷わなくてすむよう ) 標準的な RNA-Seq データ解析を一通り眺める 4

5 参考 URL 5

6 FASTQ 形式 ( と FASTA 形式 ) FASTA 形式 > ではじまる一行の description 行と配列情報からなる形式 NGS の read 長は短いので実質的に一つのリードを二行で表現 >SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT FASTQ 形式一行目ではじまる一行の description 行二行目 : 配列情報三行目 : + からはじまる一行 ( の description 行 ) 四行目 : GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT +!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 6

7 塩基配列のクオリティ情報といえば Phred スコア Phred というベースコールプログラムから得られる Quality Value(QV 値 ) のことなぜ FASTQ 形式では Phred スコアそのものでクオリティ情報を表現しないの? 7

8 理由 :( 容量 ) 節約のため Cock et al., Nucleic Acids Res., 38: , 2010 FASTQ 形式中のクオリティ情報部分 Phred スコア (QUAL 形式 ) Phred スコアが X の場合 ASCII (X+33) に対応する文字コードを割り当てる 8

9 公共 DB からデータを取得する場合 ENA Sequence Read Archive (ERA; 欧 ) FASTQ 形式でダウンロード可能 NCBI Sequence Read Archive (SRA; 米 ) (sra 形式と )sra-lite 形式でダウンロード可能 DDBJ Sequence Read Archive (DRA; 日 ) FASTQ 形式と sra-lite 形式でダウンロード可能 9

sra.lite 形式 FASTQ 形式 *.lite.sra ファイルがあるフォルダにおいて Linux システム上で以下のコマンドを実行 ( 例 : SRR002324.

10 sra.lite 形式 FASTQ 形式 *.lite.sra ファイルがあるフォルダにおいて Linux システム上で以下のコマンドを実行 ( 例 : SRR lite.sra ファイル ) fastq-dump -A SRR D SRR lite.sra 13 分程度かかる 10

11 Q & A Q: なぜ sra.lite 形式で配布するんですか? A: ファイルサイズを大幅に圧縮できるからです SRR lite.sra ファイル : 約 0.9GB SRR fastq ファイル : 約 3.8GB Q:Linux が使えないとだめってことですよね?! A:( 今のところ ) そうですねしかもそれ以外の様々な局面で Linux 環境での作業が必要 NGS 解析は Linux 上で行うのが基本 11

12 様々な選択肢があります自前で大容量メモリ計算サーバ (Linux) を購入し必要なソフトのインストールからスタート特徴 : 難易度は高いが思い通りの解析が可能 Linux サーバをもつバイオインフォ系の人にお願いする特徴 : 気軽に頼める知り合いがいればいいがその人次第 DDBJ Read Annotation Pipeline を利用特徴 : 一番お手軽な選択肢だがサポートしているプログラムのみ 12

13 13

14 入力と出力のイメージ入力データ 1: 解析したい RNA-Seq データ ( マップする側 ) 2: マッピングに使うリファレンス配列 ( マップされる側 ) マッピングプログラム (bowtie, bwa など ) 許容するミスマッチ数複数個所にマップされるものの取り扱いなど多数のオプションが存在出力結果 (SAM/BAM 形式 BED 形式など ) リファレンス配列中のどこにマップされたかという座標情報を返すのが基本形例 1:4 番染色体のの領域 ( にマップされた ) 例 2:Ensembl Gene ID XXX( の yyyy-zzzz の領域にマップされた ) 14

15 Q & A Q: クオリティスコアでのフィルタリングは? A( 一般論 ): 研究者の哲学次第 A( 私の思想 ): スコアが極端に低いものは FASTQ ファイルの段階ですでに N になっているマップされる確率が低い RNA-Seq の場合は特に気にする必要はないのでは (R で ) 塩基配列解析のウェブページ上でも Phred スコアの任意の閾値でフィルタリングするやり方を紹介しています 15

16 Q & A Q: マッピングに使うリファレンス配列は? A: 好きなものを使ってくださいゲノム配列でもトランスクリプトーム配列でも結構です Q: どこから取得できるんですか? A: UCSC Sequence and Annotation Downloads などから取得できます ( アノテーション情報も ) 以下はほんの一例ヒト全ゲノム配列の場合 ftp://genome-ftp.cse.ucsc.edu/goldenpath/hg19/bigzips/hg19.2bit ヒトトランスクリプトーム配列 (RefSeq mrna) の場合 ftp://genome-ftp.cse.ucsc.edu/goldenpath/hg19/bigzips/refmrna.fa.gz ヒトアノテーション情報の場合 16

gov/refseq/h_sapiens/ mrna_prot/human.rna.fna.

17 こんな感じヒトトランスクリプトーム配列ヒトゲノム配列 1-22 番染色体 +X+Y 約 6200 万行のファイル約 3GB のサイズ 46,XXX mrna sequences ftp://ftp.ncbi.nih.gov/refseq/h_sapiens/ mrna_prot/human.rna.fna.gz アノテーション情報ファイル (refflat 形式 ) symbol 染色体名転写開始位置 name CDS start Exon 数 strand 転写終結位置 CDS end 17

18 Q & A Q: ドラフトゲノム配列しかないんですけど A: マッピングの際のオプションで許容するミスマッチ数を増やすなどの措置をとることである程度の解析は可能だと思います Q: 手持ちの RNA-Seq データしかないんですけど A:2010 年頃から提供されはじめた de novo transcriptome assembly 用のプログラム (Trinity や Trans-ABySS など ; もちろん Linux 用です ) を利用すればトランスクリプトームの配列セット ( RefSeq のようなイメージ ) を得ることができます入力 :RNA-Seq データ出力 : コンティグ ( 転写物配列 ) >read1 GGGGTTCAAAGCAGTATCGATCAAATAGTA >read2 GTTCAAAGCAGTATCGATCAAATAGTAAAT >read3 ACGATGCAGCCTTAACGATGGTCCACAATT >read4 >contig1(transcript1) GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAA CTCACAGTTTGGAGCTTATCAGTCAA >contig2(transcript2) ACGATGCAGCCTTAACGATGGTCCACAATTATCGGGAATCA >contig3(transcript3) 18

19 マッピングの基本的なイメージ基本的なマッピングプログラム (bowtie など ) を用いた場合リファレンス配列 : ゲノム count T1 サンプルの RNA-Seq データ mapping 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 リファレンス配列 : トランスクリプトーム count 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 ゲノム配列へのマッピングの場合複数のエクソンにまたがるリード (spliced reads) はマップされないのでは?! 19

20 対策 ( リードが短かったころ ;<50bp) 既知の splice junction 周辺配列を予め組み込んだリファレンスゲノム配列側にマッピング遺伝子 1 リファレンスゲノム配列への組み込み後のイメージ >chr1 GGGGTTCAAAGCAGTATCGATCAAATAGTA >chr2 GTTCAAAGCAGTATCGATCAAATAGTAAAT > 遺伝子 1 の Exon1 の end-20bp から exon2 の start+20bp ACGATGCAGCCTTAACGATGGTCCACAATT > 遺伝子 1 の Exon2 の end-20bp から exon3 の start+20bp ( 少なくとも ) 既知の exon 間をまたぐリードのマッピングは可能 20

21 対策 ( 一リード >75bp 程度の現在 ) 再帰的にマッピングする戦略 (recursive mapping strategy) ( 通常のマッピングプログラムでマップされなかったものに対して ) リードを短くしてマップされるかどうかを繰り返す >75bp 程度のマップされなかったリードの集団 mapping 遺伝子 1 マップされない遺伝子 1 マップされない遺伝子 1 マップされた ( 原理的に ) 未知のアイソフォームの発見も可能 ~ リード長などによっても戦略が異なる ~ 21

22 マッピング結果の出力ファイル形式 ( ゲノム配列の場合 ) どの染色体上のどの位置に ( どのリードが ) マッピングされたかあるいは ( トランスクリプトーム配列の場合 ) どの転写物配列上のどの位置に ( どのリードが ) マッピングされたかを表すファイル形式 ( フォーマット ) は複数あります : BED (Browser Extensible Data) format BEDtools (Quinlan et al., Bioinformatics, 26: , 2010) GFF (General Feature Format) format SAM (Sequence Alignment/Map) format SAMtools (Li et al., Bioinformatics, 25: , 2009) マッピング結果ファイルからどうやって転写物ごとのマップされたリード数をカウントするのか? 22

23 BED 形式 23

24 BED 形式あるトランスクリプトーム配列 (RefSeq) にマップした結果転写物 ID Start End 転写物 ID ごとの出現数 = マップされたリード数 24

25 ID ごとにマップされたリード数をカウント R 上で BED 形式ファイルを入力として下記スクリプトを実行 R を用いてコピペでマップされたリード数情報を得ることができます結果ファイル (output.txt) 25

26 データ解析の前に研究目的 ( と手持ちのデータ ) をおさらい一つのサンプル内でどの転写物 (or 遺伝子 ) の発現レベルが高いか低いかを調べたい場合 RPKM や FPKM などの転写物の長さを考慮して正規化されたデータで解析トータルのリード数を補正する必要はないがやってもよい遺伝子間の発現レベルの大小関係を調べたいだけなので解析データを定数倍したところで何ら影響を与えないからサンプル間比較 (sample A vs. B など ) で発現変動遺伝子 ( Differentially Expressed Genes; DEGs) を調べたい場合トータルのリード数を補正したデータで解析正確にはサンプル間で発現変動していない遺伝子 (non-degs) ができるだけ発現変動していないと判定されるように正規化したデータ既存の R パッケージを用いて解析を行う場合には ( 整数値のみからなる ) 生のリードカウントデータを入力とし内部的に上記正規化を行う研究目的によってやっていい正規化とやってはいけない ( と言われている ) 正規化がある 26

27 正規化 ( マップされたリード数発現レベル ) 基本的な考え Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらいサンプル間の総リード数の違いをいかに補正するか配列長由来の偏り ( 長いほど沢山 sequenceされる ) をいかに補正するか ( 長さの異なる複数の isoforms が存在する場合にその遺伝子の配列長をいかに定義するか ) RPKM (Mortazavi et al., Nat Methods, 2008; ERANGEの論文 ) Reads per kilobase of exon per million mapped reads NAC (Griffith et al., Nat Methods, 2010; ALEXA-seqの論文 ) Normalized average coverage FPKM (Trapnell et al., Nat Biotechnol., 2010; Cufflinksの論文 ) Fragments per kilobase of transcript per million mapped fragments FVKM (Lee et al., Nucleic Acids Res., 2010; NEUMAの論文 ) Fragments per virtual kilobase per million mapped reads 本質的に同じ Multiple isoforms 27

28 Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらいマイクロアレイデータの正規化各サンプルから測定されたシグナル強度の和は一定と仮定チップ上の遺伝子数が少ない場合は非現実的だが数千 ~ 数万種類の遺伝子が搭載されているので妥当 ( だろう ) グローバル正規化背景 : サンプル (or chip) ごとにシグナル強度の総和は異なる対策 : 総和が任意の値 ( 例では 100) になるような正規化係数を掛ける例 :sample1 の正規化係数 = 100 /

29 RNA-Seq データの正規化 ( の一部 ) 各サンプルから sequence された総リード数は一定と仮定 T1 Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらい遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 RPM 正規化 Reads Per Million mapped reads(rpm) 正規化後の総リード数が 100 万 (one million) になるように補正例 :T1 の正規化係数 = / 67 29

30 配列長の補正 Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらい Mortazavi et al., Nature Methods, 5: , 2008 配列長が長い遺伝子ほど沢山 sequence されるそれらの遺伝子上にマップされる生のリード数が増加傾向配列長が長い遺伝子ほど発現レベルが高い傾向になる発現レベルが同じで長さの異なる二つの mrnas AAAAAAA AAAAAAA 断片化して sequence マップされたリード数をカウント AAAAAAA AAAAAAA 30

31 配列長の補正 Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらい Mortazavi et al., Nature Methods, 5: , 2008 前提条件 : 配列長が既知補正の基本戦略 : 配列長で割る AAAAAAA AAAAAAA 1 / 配列長を掛ける場合塩基あたりの平均のリード数を計算しているのと等価 1000 / 配列長を掛ける場合その遺伝子の配列長が1000bpだったときのリード数と等価 Reads Per Kilobase of exon 31

RPKM RPM 正規化 ( マイクロアレイなどと同じところ ) Reads per million mapped reads サンプルごとにマップされた総リード ( 塩基配列 ) 数が異なる各遺伝子のマップされたリード数を総 read 数が100 万 (one million) だった場合に補正 raw counts:all reads= RPM : 1,000,000 A1BGの場合は

32 RPKM RPM 正規化 ( マイクロアレイなどと同じところ ) Reads per million mapped reads サンプルごとにマップされた総リード ( 塩基配列 ) 数が異なる各遺伝子のマップされたリード数を総 read 数が100 万 (one million) だった場合に補正 raw counts:all reads= RPM : 1,000,000 A1BGの場合は 744 : 5,087,097 = RPM : 1,000,000 1,000,000 1,000,000 RPM raw counts all reads 5,087,097 RPKM 正規化 (RNA-Seq 特有 ) Illumina webinar session 1(2011 年 9 月 8 日開催 ) のおさらい Mortazavi et al., Nature Methods, 5: , 2008 Reads per kilobase of exon per million mapped reads 遺伝子の配列長が長いほど配列決定 (sequence) される確率が上昇各遺伝子の配列長を 1000 塩基 (one kilobase) の長さだった場合に補正 RPKM 1,000,000 1,000 raw counts all reads gene length 1,000,000,000 raw counts gene length all reads RPM 32

33 Trinity 出力結果から FPKM 値を取得 Trinity (Grabherr et al., Nat Biotechnol., 29: , 2011) RNA-Seq データを入力としてトランスクリプトームのアセンブリを行うプログラム出力は multi-fasta 形式のファイル ( ファイル名 :Trinity.fasta) 転写物 ( コンティグ ) の塩基配列情報 description 行に配列長や発現レベルに相当する FPKM 値が含まれる R を使って簡単に FPKM 値の情報を取得することができます 33

34 利用可能な R パッケージたち DEGseq (Wang et al., Bioinformatics, 26: , 2010) ポワソン分布 (variance = mean) を仮定しているためばらつきを過少評価 edger (Robinson et al., Bioinformatics, 26: , 2010) 正規化法 :TMM 法負の二項分布 (variance > mean) を仮定 mean のみのパラメータを用いて現実のばらつきを表現 DESeq (Anders and Huber, Genome Biol., 11: R106, 2010) 正規化法 :RLE 法 (relative log expression) edger のモデルをさらに拡張 ( しているらしい ) bayseq (Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010) 正規化法 :RPM ( たぶん ) 配列の長さ情報を与えるオプションがあるデータセット中に占める DEG の割合 (P DEG ) を一意に返す NBPSeq (Di et al., SAGMB, 10:24, 2011) Oct 入力 : 生のリードカウントからなる遺伝子発現行列出力 : 遺伝子ごとの発現変動の度合い (p 値など ) 34

35 理想的な実験デザイン ( 二群間比較 ) サンプル A vs. B の比較 (Kidney vs. Liver; 腎臓 vs. 肝臓 ) 生のリードカウントのデータ ( 整数値 ) Biological replicates のデータ生物学的なばらつき ( 個体間の違い ) を考慮すべし A1: ある生物の腎臓 A2: 同じ生物種の別個体の腎臓 A3: 同じ生物種のさらに別個体の腎臓 B1: ある生物の肝臓 B2: 同じ生物種の別個体の肝臓 35

36 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) liver( 肝臓 ) Technical replicates のデータサンプル内の技術的なばらつき ( 例 : レーン間の違い ) の度合いを調べるためのデータでありこのようなデータで二群間比較し発現変動遺伝子がどの程度あるかといった数に関する議論は無意味解析例 : アリエナイ?! 数 (50% とか ) が発現変動遺伝子として検出される理由 :Biological variation > Technical variation 36

37 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) RPM 正規化 1,000,000 12,685 1,804,

分布の話例題 :Marioni et al., Genome Res., 18: 1509-1517, 2008 のデータ ( の一部 ) kidney( 腎臓 ) adjusted R-squared: 0.

38 分布の話例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) kidney( 腎臓 ) adjusted R-squared: y = x y = a + bx Technical replicates のデータは : ( 遺伝子の )VARIANCE はその MEAN で説明可能である VARIANCE MEAN ポアソン分布に従うポアソンモデルが適用可能 38

39 分布の話生物アイコン ( 例題 :Cumbie et al., PLoS ONE, 6: e25279, 2011 のデータ ( の一部 ) Arabidopsis( シロイヌナズナ ) adjusted R-squared: y = a + bx y = x Biological replicates のデータは : VARIANCE > MEAN 負の二項 (NB) 分布に従う NB モデルが適用可能 39

40 なぜ沢山の方法が存在しうるのか? DEGseq (Wang et al., Bioinformatics, 26: , 2010) ポワソン分布 (variance = mean) を仮定しているためばらつきを過少評価 edger (Robinson et al., Bioinformatics, 26: , 2010) 正規化法 :TMM 法負の二項分布 (variance > mean) を仮定 DESeq (Anders and Huber, Genome Biol., 11: R106, 2010) 正規化法 :RLE 法 (relative log expression) edgerのモデルをさらに拡張 ( しているらしい ) bayseq (Hardcastle and Kelly, BMC Bioinformatics, 11:422, 2010) 正規化法 :RPM ( たぶん ) 配列の長さ情報を与えるオプションがあるデータセット中に占めるDEGの割合 (P DEG ) を一意に返す NBPSeq (Di et al., SAGMB, 10:24, 2011) VAR (1 ) VAR (1 ) Ans. Variance と Mean の関係を表現する手段が沢山あるから VAR (1 1 ) VAR Oct

41 edger を使ってみる例題 :Marioni et al., Genome Res., 18: , 2008 のデータ kidney( 腎臓 ) liver( 肝臓 ) ファイル名 :SupplementaryTable2_changed.txt 内容 :A 群が最初の 5 列 B 群が残りの 5 列のデータ解析結果を hoge2.txt という名前でファイルに出力したい 41

42 edger を使ってみるファイル名 :SupplementaryTable2_changed.txt 内容 :A 群が最初の 5 列 B 群が残りの 5 列のデータ解析結果を hoge2.txt という名前でファイルに出力したい 42

43 edger を使ってみる R 上でスクリプトをコピペ! ( エラーメッセージが出ていなければ hoge2.txt というファイルができているはず ) 43

44 edger を使ってみる一番右側の数値が False Discovery Rate (FDR) この列 (O 列 ) で昇順にソートすれば任意の閾値を満たす遺伝子数がわかる 9,862 個が FDR < 0.01 を満たす 11,172 個が FDR < 0.05 を満たす 44

45 edger を使ってみる Top-ranked gene の生リードカウントを眺めても確かに発現変動 (Kidney << Liver) していることが分かる 45

46 edger を使ってみる M-A plot を描画 (FDR < 0.01 を満たすものを赤色で表示 ) hoge2.png 9877 個 ( 全遺伝子数のうち約 31% が FDR < 0.01 を満たす ) 46

47 edger を使ってみる M-A plot を描画 (2 倍以上発現変動しているものを赤色で表示 ) hoge2.png 個 ( 全遺伝子数のうち約 37% が 2 倍以上発現変動している ) このやり方はダメなんです 47

48 倍率変化がだめな理由をデモ例題 :Marioni et al., Genome Res., 18: , 2008 のデータ kidney( 腎臓 ) liver( 肝臓 ) 発現変動遺伝子がないデータで二群間比較をしてみる A 群 B 群 48

49 倍率変化がだめな理由をデモ例題 :Marioni et al., Genome Res., 18: , 2008 のデータ ( の一部 ) (A1, A2) vs. (A3, A4) の二群間比較結果 edger で FDR < 0.01 を満たすものは 0 個 (edger で )2 倍以上発現変動しているものは 3813 個低発現領域で log 比が大きくなる現象をうまくモデル化することが重要 49

50 まとめ RNA-Seq データ取得から標準的なデータ解析の流れを説明一般的なファイル形式 (FASTQ) について説明一通りの解析を自力で行うためにはLinux 系スキルが必要マッピング ( やアセンブル ) 以降は基本的にRで解析可能研究目的によってデータ解析時の入力データが異なるサンプル間比較 : 生のリードカウントデータサンプル内比較 : 長さ補正を行ったデータ (RPKMやFPKMなど) 分布を考えることは重要 (DEG 数を議論したい場合 ) technical replicatesやbiological replicates Rパッケージを用いれば発現変動遺伝子の検出から描画まで簡単二倍( 倍率変化 ) じゃだめなんですさん (Rで) 塩基配列解析のウェブページを用いてなるべく自力で解析 50

51 Top 400 Top 2000 低い全体的な発現レベル高い 51

52 adjusted R-squared: y = x y = a + bx adjusted R-squared: RPM 正規化データ RPKM 正規化データ 52

機能ゲノム学（第6回）

機能ゲノム学（第6回） R でトランスクリプトーム解析東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 ( かどたこうじ ) http://www.iu.a.u-tokyo.ac.jp/~kadota/ kadota@iu.a.u-tokyo.ac.jp 1 自己紹介 2002 年 3 月東京大学大学院農学生命科学研究科博士課程修了学位論文 : cdna マイクロアレイを用いた遺伝子発現解析手法の開発