Platform for Drug Discovery


[日本語]  

cuffdiffOutToHtmlWithChIPpeakAnnotation(for multi comparison, add extra inf from GTF, with link for GE list, max 10 peaks)


Introduction


  • このパイプラインはRNA発現解析結果(Cuffdiff出力)に対して、ChIP-seqによるピークの有無(転写因子等の結合領域の有無)の情報を付加することで、両方の情報を利用した絞込み等を可能にします。
  • 多種類の転写因子に関心がある場合や、ピーク検出の妥当性を複数手法を比較して参照したい場合などに利用が期待できます。

  • Input1 formatcuffdiffOutput
    Input2 formatgtf
    Input3-12 formatChIPpeakAnno output
    Input13-32 formatge_db
    Library layoutAny library layout including single-end and paired-end
    SpeciesAny species with Genome sequence and with known gene annotation
    Execution timeAbout 20 mins.

    Results
    Image

Table of contents

Inputs


入力の補足情報(RNA-seq実行結果の準備)

  • RNA-seqの結果にChIP-seqの結果を統合するのに先立ち、RNA-seqの解析を実行しておく必要があります。
  • 以下のパイプラインなどで生成される"cuffdiffOutput"形式のデータの準備が必要です。入手方法は後述の 入力1: cuffdiffOutput (Input cuffdiff output files)に詳しい記載があります。
  • ここでは、フローが複雑になりすぎるのを防ぐためプロジェクトの部分的コピー機能を使って以下のデータのみ抽出しています。
    • 1. RNA-seq配列のゲノムへのマッピング結果のbamファイル(必須ではありません。)
    • 2. 既知アノテーション情報のgtfファイル
    • 3. RNA-seqから予測された新規遺伝子候補を含む遺伝子構造情報gtfファイル (必須ではありません)
    • 4. 上記1~3をゲノムビューワGenomeExplorer用にトラックとしてロード(読み込み)したもの
    • 5.Cuffdiffの出力結果
    • 6. 上記5を整形した簡易絞込み等可能なレポート。(必須ではありません。)
  • 結果として以下のような実行履歴ができます。
    • 下図の橙色の注釈で示されたデータがこのRNA-seqとChIP-seq統合パイプラインで入力が必須となるデータです。
  • 解析フローが複雑なので模式図で表現すると以下のような処理が行われていることになります。

入力の補足情報(ChIP-seq実行結果の準備)

  • RNA-seqの結果にアノテーション情報として付加するためのChIP-seqの解析結果を準備します。
  • この例では同一のChIP-seqデータから3種の異なるピークコーラーにより転写因子結合領域予測を行った結果を利用します。
  • それぞれの予測結合領域情報に対して、領域近傍の既知遺伝子アノテーション付加を行っています。
  • RNA-seqとの統合解析を行うには、このChIP-seq解析で使用する遺伝子アノテーション情報をRNA-seq解析と同一ファイルにするのが無難です。
    • なぜなら遺伝子情報はバージョンによってIDの変更や追加などが発生するのが一般的です。
    • 全く異なる由来の遺伝子情報の場合ID体系が異なり対応付けができない場合もありえます。
    • ここでは遺伝子アノテーション情報ファイルとして一般的なgtfファイルの例を記載しています。
      • gtf形式のアノテーション情報の入手方法はこちら
      • ChIP-seqのアノテーション付加にはあらかじめbed形式に変換する必要があるため、以下のパイプラインを利用しています。
      • パイプライン gtf->BED
    • アノテーション情報としてgtfファイルでなく、UCSCから提供されているgenePred(refFlat)形式のファイルを利用することも可能ですが、予め以下のような変換が必要です。
      • genePred(refFlat)形式のアノテーション情報の入手方法はこちら
      • RNA-seq用 パイプラインrefFlat<genePred>->gtf
      • ChIP-seq用 パイプラインrefFlat->BED12
  • この例ではChIP処理配列と、Control配列をそれぞれ以下のパイプラインでマッピング後、ユニークヒットなもののみ取り出しています。
  • その後以下の3種のパイプラインを用いて、ピークコール(結合領域予測)を行います。
  • 最後に以下のパイプラインを用いて、予測結合領域に遺伝子アノテーションを付加します。
  • 結果として以下のような実行履歴ができます。
    • 下図の橙色の注釈で示されたデータがこのRNA-seqとChIP-seq統合パイプラインで入力が必須となるデータです。
    • この例ではアノテーション付予測結合領域データが3つですが、1~10個の任意個数指定できます。また、この例のように同一ChIP-seqデータ由来の必要もありませんし、フローが繋がっている必要もありません。
  • 解析フローが複雑なので模式図で表現すると以下のような処理が行われていることになります。

RNA-seqプロジェクトとChIP-seqプロジェクトのマージ

入力1: cuffdiffOutput (Input cuffdiff output files)

入力2: gtf (Annotation GTF)

入力3~12: ChIPpeakAnno output (ChIP peak with gene annotation)

  • 上述の入力の補足情報(ChIP-seq実行結果の準備)において"[ピークコーラ何某]検出ピーク領域 アノテーション付加後"と表現されているデータです。
  • bed形式の領域情報に対して以下のパイプラインを実行した結果です。
  • 説明例ではChIP-seqの異なる3つのピークコーラーによるピーク領域(この場合は転写因子の予測結合領域)に遺伝子アノテーションを付加したデータを利用していますが、以下のような入力も指定が可能です。
    • 異なる解析条件(パラメータの変更等)で実行した予測結合領域情報
    • 異なる抗体を用いたChIP-seqの解析結果
    • 有意なメチル化率の変動を持つbisulfite-seq由来などのCpG領域情報
    • その他遺伝子近傍に付加される情報
  • 最大20個まで指定できるゲノムビューワ(GenomeExplorer)用のトラックデータです。
  • 指定は任意で指定しなくても実行結果は返されます。
  • 出力結果のレポート中のGenomeExplorerへのリンクを利用する際に、指定された順番のトラックデータが選択された状態で起動します。
  • 結果のリストに関係するRNA-seqのマッピングデータや、遺伝子アノテーション情報、ChIP-seqのマッピングデータ、検出されたピーク領域のデータに関するトラックを追加すると良いでしょう。
  • この例ではHow to run this pipelineにも記載された以下の順でトラックを選択してます。
    • #トラック説明
      1GenomeExplorer用 既知遺伝子トラック
      2GenomeExplorer用 予測遺伝子構造トラック
      3GenomeExplorer用 RNA-seqサンプル1トラック
      4GenomeExplorer用 RNA-seqサンプル2トラック
      5GenomeExplorer用 RNA-seqサンプル3トラック
      6GenomeExplorer用 ChIP(CASE)トラック
      7GenomeExplorer用 ChIP(Control)トラック
      8GenomeExplorer用 MACS検出ピークトラック
      9GenomeExplorer用 PeakSeq検出ピークトラック
      10GenomeExplorer用 SISSRs検出ピークトラック

Pipeline hisory


  • Image
  • 上の点線で囲まれた範囲がこのパイプラインの取り扱い範囲になります。

Outputs

(1) Cuffdiff with ChIPpeakAnno RNA-seq and ChIP-seq integrated report

  • プログラミング言語が使えない方でも、ChIPのピーク有無と、RNA-seqの発現解析結果を組み合わせて、絞込みを可能にしたインターフェイス付きのレポートを返します。
  • 絞り込んだ結果を幣所開発のGenomeExplorerで確認するなどしてより詳細な状況確認を行うことができます。
    • Image
    • ビューワ上で"MAP2K"遺伝子がサンプル3のみ発現が殆どみられないことや、近傍の"ABCA10"遺伝子が異なる発現パターンを見せていること、"MAP2K"遺伝子周辺に複数のChIPによる予測結合部位が見られることが確認できます。
    • Image
    • 予測結合領域を確認することで馴染みのないChIP-seqのピークコールプログラムの予測傾向なども把握することができるでしょう。
  • レポートの機能と動作上の注意について
    • ファイルサイズが大きい場合Internet Explorerでは警告メッセージがでて閲覧できない場合がありますが、google chromeなどのブラウザでの閲覧を推奨します。
    • ブラウザを通じてメモリを大きく使用するため、古いパソコンなどで閲覧する場合は、他のアプリケーションを閉じるなどの対応をとって下さい。
  • Image
    • 簡易フィルタ機能を含むhtml版と、Excel等のスプレッドシートで開くことのできるcsv版が提供されていますが、本質的にはデータの内容は同一です。
    • 以下のように利点と難点がありますので、状況に応じて使い分けて下さい。
      #ファイル名の例利点難点
      1html版(何某).html簡易的なフィルタリングが容易に行えます。検定対象のエントリー数が多いと表示でません。
      Internet Explorer非対応。
      データの加工は困難。
      2csv版(何某).csvExcel等のスプレッドシート用アプリケーションで開き複雑な絞り込みや、データの加工ができます。
      データサイズが増えても大きな問題ありません。
      開くまでに手間がかかります。
      3共通点-ダウンロードしてオフラインで利用できます。
      ファイルそのものにデータを含んでいます。
  • 列の説明
    #列名説明
    1tss_idCuffdiffが付番する一意な転写開始点ID
    2gene_idCuffdiffが付番する一意な遺伝子ID
    3gene既知遺伝子名(複数関連づく場合はカンマ区切り)
    4locus当該遺伝子座標
    5<RNA-seq A>/<RNA-seq B> Up/downRNA-seq群Bに対する群Aの増減および、検定結果を組み合わせた記号。詳しくは下表Up/Down参照。
    6<RNA-seq A>/<RNA-seq B> sample_1RNA-seq群Aのサンプルラベル
    7<RNA-seq A>/<RNA-seq B> sample_2RNA-seq群Bのサンプルラベル
    8<RNA-seq A>/<RNA-seq B> statusRNA-seq群Aと群Bの比較における検定の実施有無です。OK:検定が正常に実施、NOTEST:検定に必要なアライメントが充分得られない、LOWDATA:遺伝子構造が複雑すぎるか配列が薄すぎる、HIDATA:領域内に配列が多すぎる、FAIL(その他検定に失敗した場合)。
    9<RNA-seq A>/<RNA-seq B> value_1RNA-seq群AのFPKM値
    10<RNA-seq A>/<RNA-seq B> value_2RNA-seq群BのFPKM値
    11<RNA-seq A>/<RNA-seq B> log2(fold_change)RNA-seq群Aに対する群Bのlog2 fold change値。
    12<RNA-seq A>/<RNA-seq B> test_statRNA-seq群Aと群Bの比較における統計量です。
    13<RNA-seq A>/<RNA-seq B> p_valueRNA-seq群Aと群Bの比較における検定結果のp値です。
    14<RNA-seq A>/<RNA-seq B> q_valueRNA-seq群Aと群Bの比較におけるFDR補正後のq値です。
    15<RNA-seq A>/<RNA-seq B> significantRNA-seq群Aと群Bの比較における有意差の有無です。yes:有意差あり、no:有意差無し。
    16chrom染色体名。
    17strPosゲノム中の領域の開始座標。
    18endPosゲノム中の領域の終了座標。
    19GE linkGenome Explorerへのリンク。
    20USCS linkUCSC Genome Browserへのリンク。
    21Entrez gene linkEntrez Geneへのリンク。ID列と、オプションで指定した生物種をキーに問い合わせた結果のページにリンクしています。
    22ChIP A peak digestChIP Aにおけるピークの有無のダイジェスト。ピーク名と遺伝子からの相対位置関係を含む。(複数ピークある場合は複数行表示)
    23ChIP A peak detailChIP Aにおけるピークの有無の詳細情報。ChIPpeakAnnoの出力に含まれる値を":"区切りで連結したのものです。※
    24以降列名不定アノテーション情報のgtfファイルに含まれるgene_idとtranscript_id以外の情報が順番に格納されます。
    • Up/downの説明
      notEnoughstatus列が"OK"以外の場合です。
      sigUpstatus列が"OK"でsignificant列が"yes"でlog2(fold_change)列が正の値の場合。
      upstatus列が"OK"でsignificant列が"no"でlog2(fold_change)列が正の値の場合。
      sigDownstatus列が"OK"でsignificant列が"yes"でlog2(fold_change)列が負の値の場合。
      downstatus列が"OK"でsignificant列が"no"でlog2(fold_change)列が負の値の場合。
      equalstatus列が"OK"でlog2(fold_change)列が0の場合。
  • GenomeExplorerへのリンク
    • html版からのGEへの遷移
      csv版からのGEへの遷移
    • html版、csv版どちらも、右端近くの"GE link"列からGenome Explorerに遷移することができます。
    • html版ではwebページ上の"Link"ボタンを押せば別画面でサンプル選択済みのGenomeExplorerが起動します。
    • csv版ではLinkマスにあるURLをコピーしてブラウザで開くとサンプル選択済みのGenomeExplorerが起動します。
    • 青線と、緑線矢印にあるように、リンクのあった行が表示されていることが確認出来るかと思います。
    • Genome Explorerでの参照方法は、こちらをご参照ください。
    • 現在のバージョンでは、GenomeExplorerへの遷移の度に、初期表示に戻ってしまうため、複数の領域を逐次的に見たい場合は、GenomeExplorerのsearch機能で、場所だけ移動すると便利です。
  • 開閉式の簡易フィルタと列説明
      • html版では、画面の左上開閉式の簡易フィルタの指定箇所(column filter)と、列の説明文(legend)があり、マウスでクリックするたびに開閉します。
  • 簡易フィルタとソート
    • 簡易フィルタ
      ソート
    • 画面左上の"column filter"をクリックすると開かれる選択ボックスやスライダーを利用すると、結果列をフィルタリングすることができます。
    • フィルタの例
      • 遺伝子名を"MAP2"を含むものにフィルタ
      • Up/down列が有意かつ発現増加または有意かつ発現減少のものにフィルタ
      • RNA-seq1とRNA-seq2の比較のfoldChange値が3以上のものにフィルタ
      • 染色体13番の遺伝子のみにフィルタ
      • 3手法の全てでChIPピークが検出された遺伝子にフィルタ
    • html版の結果の任意のカラムラベルをクリックすると、その列で結果のソートが出来ます。
      • 図中の例では、FPKM値による降順(灰色の下三角が表示)と、log2(fold_change)値による昇順(灰色の上三角が表示)を記載しています。

How to run this pipeline

  • 解析実行の流れの詳細はこちらをご確認ください。
  • 実行するプロジェクト内に入力となるアノテーション付加対象bed領域情報ファイルと、アノテーションファイル(bed形式)を準備しておきます。
  • アップロードしたファイルのデータアイコンの下の"Select"ボタンを押した後、ポップアップするサブ画面中の"Analysis"ボタンを押します。
    • 上図のカッコ内数字の順にデータを選択していきます。
    • Image
    • データを選択していくと、上のようなデータ選択ダイアログが表示されるので、データ選択が終了したら、このダイアログ中の"Analysis"ボタンを押します。
  • 選択されたファイルのデータ型に応じた解析の選択肢がリストされるので、その左のメニューから"Tools"-"Integration"を選び、このパイプラインを選びます。
  • 選択された入力ファイルの確認と、プロジェクト(ワークスペース)に誤りが無いかを確認します。
    • 特にアノテーション付加後のピークデータ(ChIPpeakAnno output)と、GenomeExplorer用のトラック(ge_db)が意図した順番になっているか確認してください。意図しない順番になっていた場合は、ドラッグアンドドロップで置き換えができます。
  • レポート出力のカスタマイズなどの設定が必要がない場合はオプションを設定せずに"Run"することもできますが、以下のような場合は、"Set option and run"を押して、オプション選択画面に遷移します。
    • レポート中のラベルを制御したい場合。
    • ヒト以外の生物で実行する場合。
    • GenomeExplorerやEntrez geneへのリンクを有効にするため、生物種等の指定が必要な場合。
    • その他、解析のカスタマイズを行いたい場合。
  • オプション選択についての詳細は下のOptionsの項目をご覧ください。

Options


    Image

    • 設定可能な実行オプション値と説明は以下にあります。
    • 特にオプション自由入力欄に値を設定する場合は間違いやすいので、「解析条件の確認と実行ログの確認について」を参照して、トラブルシューティングする必要があるでしょう。
    • #オプション番号オプション名重要性説明
      11.1cuffdiff report labelsレポート列中のRNA-seqのサンプル名のペアを{}で囲った上でカンマ区切りで指定します。Cuffdiffの実行時にラベル指定している場合は設定する必要はありません。この例では'{RNA-seq1}{RNA-seq2},{RNA-seq1}{RNA-seq3},{RNA-seq2}{RNA-seq3}'を指定しています。省略時はCuffdiffのデフォルトラベルのq1,q2,q3,…が使用されるため、3サンプル総当りの今回は'{q1}{q2},{q1}{q3},{q2}{q3}'と同義です。
      21.2ID process type-アノテーション付けされたChIP-seqの情報の複合IDから、gene_idやtranscipt_idを抜き出す場合は指定します。分らない場合は初期値のgene_idを指定します。Cuffdiff中に付加されるアノテーションはgtfファイル中のgene_idが使用されますが、ChIP-seqのアノテーション付けでgene_idとtranscript_idを複合したIDを使用している場合は、その複合ID中からgene_idを抜き出す処理を実施した上でCuffdiff結果とChIP-seq中のアノテーション情報の対応付けが行われます。
      31.3convertChIPpeakAnnoForIntegration.pl --out_type-レポートに含めるChIPアノテーション行を各サンプルごと、DIGEST_ONLY:ダイジェスト列のみ表示、DETAIL_ONLY:詳細列のみ表示、BOTH:ダイジェスト列と詳細列の両方表示から選択します。
      41.4AddChIPpeakAnnotation labelsChIPピークの入力それぞれに対する列ラベルのプレフィックスをカンマ区切りで指定します。この例では"MACS,PeakSeq,SISSRs"としています。省略時には、'Data1,Data2,..,DataN'が使用されます。
      51.5GENOME_ID-ゲノムのバージョンで、レポート中のGenomeExplorerへのリンクに使用されます。GenomeExplorerを利用しない場合はこのオプションは何が設定されていても問題ありません。
      61.6Gene ID type-特に理由が無ければAUTOを指定します。レポート中のEntrez Gene IDのリンクを付加する際、Cuffdiff出力からgene_idの情報を取得し利用します。特別な加工を行ったgtfでない限りは変更の必要はありません。
      71.7Organism for Entrez Gene Link-レポート中のEntrez Gene IDのリンクを付加する際、どの生物種のIDかを指定するのに利用します。省略可ですが、既知遺伝子情報GTFファイル中のgene_idに利用されることの多い、geneSymbolは複数の生物種で共通なので、生物種を指定しないと、リンク先にすべての生物種の当該遺伝子の詳細がリストされてしまいます。ヒトであれば学名の'Homo sapiens'と指定しておくと、直接ヒトのエントリーのみがリンクから開けます。入力が省かれた場合は、リンクをクリックすると全生物のリストが表示されることになるでしょう。
      81.8tsvFromCuffdiffOutToHtmlAfterMakeLink suppress_compare_detail-レポートに発現解析の比較の詳細列を含めない場合にチェックをオンにします。
      91.9tsvFromCuffdiffOutToHtmlAfterMakeLink suppress_position-レポートに位置情報の列を含めない場合にチェックをオンにします。
      101.10tsvFromCuffdiffOutToHtmlAfterMakeLink suppress_links-レポートにリンク列を含めない場合にチェックをオンにします。
      111.11tsvFromCuffdiffOutToHtmlAfterMakeLink suppress_others-レポートにその他の列を含めない場合にチェックをオンにします。その他の列とは、入力2番目に指定されたgtfファイル中のgene_id,transcript_id以外の付加情報などです。
      121.12tsvFromCuffdiffOutToHtmlAfterMakeLink infinity_value-htmlレポート作成の際に、スライダーによる数値フィルタを有効にするためには、無限大の値が含まれていてはエラーになるため、無限大を特定の大きい数値に置き換えています。このオプションでは置き換える数値を指定します。デフォルトでは9999に置き換わります。

Comments

  • このパイプラインでは、予測結合領域を幣所開発のゲノムビューワ(GenomeExplorer)にインポートして可視化するために、予めビューワに登録された生物種のゲノムバージョンを選択する必要がありますが、この機能を無視し、使用する生物種に対応するアノテーション情報と、ピークコール結果を使用すれば、その他の生物でも利用できるはずです。ゲノムの選択肢があり、利用できない場合はご連絡下さい。

Use case

  • この説明ページで使用されているデータは、以下の処理を施しています。
    • RNA-seqデータは、ENCODE ProjectのRNA-seqデータ(SRA039973)で、多数あるサンプルの内、説明用に3サンプル抜き出したものを、パイプラインTopHat-Cufflinks-Cuffmerge-Cuffdiff for novel transcript detection and expression comparison pipeline (less than 10 N=1 paired-end fastq files)で処理したものです。
    • #SRX番号サンプル名含まれるレーン番号
      1SRX084666GSM765388: CshlLong_RnaSeq_MCF-7_cell_longPolyASRR315301,SRR315302
      2SRX082565GSM758559: CshlLong_RnaSeq_GM12878_cell_longPolyASRR307897,SRR307898
      3SRX084683GSM765405: CshlLong_RnaSeq_K562_cell_longPolyASRR315336,SRR315337
    • ChIP-seqデータはENCODEプロジェクトが公開しているGSM1010826細胞株での転写因子FOXA1のChIP-Seqデータ(SRR577861(external link))とそのコントロールデータ(SRR577663(external link))を入力として、MACS, PeakSeq, SISSRsのそれぞれのデフォルトオプションで検出されたピークに対してiGenome由来の遺伝子アノテーション付加を行っています。

Related Tools