Platform for Drug Discovery


TopHat2, CuffLinks2 and CummeRbund + GE (PE, with detection of de novo splicing variants)


Introduction



    Input formatFASTQ (Illumina, IonTorrent, 454)
    Library layoutFor paired-end read only
    SpeciesHuman(hg19), Rat(rn4), Mouse(mm9), Drosophila(dm3)
    Execution timeAbout a few days (20M paired-end reads[50bp+50bp])


    Results
    • GOアノテーション付き発現量テーブル
    • クラスタリング図 ...etc.
    Image
    クラスタリング結果

Table of contents




Inputs


  • ファイルはペアエンドのFASTQファイルを2~10サンプルの範囲で入力可能です。
  • FASTQについては、このページをご覧ください。
  • アップロードするファイル名にスペース「 」や、カッコ「(」、クオーテーション「'」などといった文字や日本語は使えません。(アップロードできても、解析時にエラーになります。) アルファベット+数字+アンダーバー「_」で名前を付けてください。
  • ペアエンドの場合、フォワード側は「xxx_1.fastq」、「xxx_1.fq」、「xxx_1.txt」、「xxx_1_sequence.txt」、「xxx_R1_yyy.fastq」などとなっている必要があり、リバース側は「xxx_2.fastq」、・・・となっている必要があります。
    (例:my_sample_1.fastq, my_sample_2.fastq)
  • 各FASTQファイルがそれぞれ50GB以下であれば、経験的には当サーバで解析可能です。ただし、ファイルサイズによって、解析時間が延びるため、手早く結果を見るには、それぞれ数GB程度の大きさにダウンサンプリングするとよいかもしれません。
  • 1リードの長さが50 bp以下の場合、オプションに「--no-coverage-search」と付けないと、計算時間が非常に長くなることがあります。

Outputs


    実行例1
    Image


    1.解析レポート(Genome Explorer)
    Image
    TopHatによるマッピング結果をGenome Explorerで可視化
    解析レポート(CummeRbund Density Plot)
    Image
    各サンプルの発現量の分布
    解析レポート(CummeRbund Dendrogram)
    Image
    サンプル間の類似性
    解析レポート(CummeRbund Scatter Plot)
    Image
    サンプル間の総当たり散布図
    解析レポート(CummeRbund Volcano Plot)
    Image
    サンプル間の総当たりVolcano Plot
    解析レポート(CummeRbund Cluster Plot)
    Image
    クラスタリング結果
    解析レポート(GOseqレポート)
    Image
    GOアノテーション後のHTMLレポート
    解析レポート(REVIGO)
    Image
    GOseq→REVIGOでGO Termの可視化
    解析レポート(発現量テーブル)
    Image
    発現量一覧(tsvファイル(Excelにて閲覧可能))
    2.TopHat-CuffLinks付随情報(スプライスジャンクション)
    Image
    TopHatのSplice junction部位リスト(BED形式)
    TopHat-CuffLinks付随情報(Insert)
    Image
    TopHatの挿入塩基部位リスト(BED形式)
    TopHat-CuffLinks付随情報(Deletion)
    Image
    TopHatの欠失塩基部位リスト(BED形式)



How to run this pipeline


    1.まずはMASERにログインし、Projectの一覧を表示させます。
    Image
    そして、「Create New Project」をクリックし、新しいプロジェクトを作成します。

    2.プロジェクトの名前を付け(ここでは仮に「human RNA-seq」とします)、OKを押します。
    Image

    3.作成したプロジェクトをクリックし、ページを開きます。
    Image

    4.プロジェクトを作成したら、次は入力ファイル(ここではFASTQファイル)をアップロードします。「Upload Data Here」をクリックしてください。
    Image

    5.ファイルの転送プロトコルを選択します。通常はHTTPSのほうをクリックすれば大丈夫です。
    Image

    6.Dataのラベル名(ここではSRX082565M1.2)、Dataの種類(ここではfastq (paired-end))、そしてアップロードするファイル(ここではC:\SRX082565M1.2_1.fastq)を選択します。
    Image
    paired-endのfastqはforwardとreverseの2つのファイルがセットですので、「Add file」ボタンを押して、項目を追加します。

    7.参照ボタンを押して、もう一つのfastqファイルを選択してください。
    Image
    アップロードするファイル名には次のような制限があります。
    スペース「 」や、カッコ「(」、クオーテーション「'」などといった文字や日本語は使えません。アルファベット+数字+アンダーバー「_」で名前を付けてください。
    ・ペアエンドの場合、フォワード側は「xxx_1.fastq」、「xxx_1.fq」、「xxx_1.txt」、「xxx_1_sequence.txt」、「xxx_R1_yyy.fastq」などとなっている必要があり、リバース側は「xxx_2.fastq」、・・・となっている必要があります。
    (例:my_sample_1.fastq, my_sample_2.fastq)
    詳細はこちらをご覧ください。

    また、FASTQファイルのクオリティにはいくつか種類があります。自分のデータがどの種類なのか分からない場合には、FASTQCなどのツールを先に使用して下さい。FASTQCの「Basic Statistics」項目で「Sanger / Illumina 1.8」などと表示されれば特に気にすることはありませんが、もしもIllumina 1.3 - 1.7の場合には、後述のtophatのオプション設定時に「--phred64-quals」を指定する必要があります。
    それと、基本的にはペアの向きがhead-to-head(→←の向き)でnon-directionalであることを想定しています。もしもメイトペアの場合や、directionalである場合には、適切なオプションを設定してください。

    8.ほかにも発現量を比較したいサンプルがある場合、同様にしてFASTQファイルをアップロードします。
    今回は、3つのデータをアップロードしました。
    データの下に表示されている「Select」ボタンを3つとも順にクリックしてください。
    クリックした順番に後ほど結果が並ぶので、見たい並び順となるように選択してください。
    Image
    もしも用意されたアノテーションとは別のものを使用したい場合(EnsGeneからGENCODEに変えたいなど)には、最初にGTFファイルもアップロードしておき、同時に選択しておいてください。
    また共生細菌など、宿主のゲノムとは別に追加したいリファレンス配列があれば、マルチFASTA形式のファイルを1つだけリファレンスに追加することが可能です。その場合には、FASTAファイルも同時に選択しておいてください。

    9.別枠に選択されたファイルの一覧が表示されているかと思います。
    Image
    「Analysis」ボタンをクリックします。

    10.「RNA-seq」を選択し、本パイプライン名の横の「Analysis」をクリックしてください。
    Image

    11.下にスクロールすると、入力ファイルの一覧が表示されます。
    Image
    「Set option and run」をクリックして下さい。

    12.オプションで変更可能な項目を説明します。
    必ず確認しなければいけないのは、使用するリファレンスの名前(ヒトであればhg19, マウスであればmm10など)です。また、アノテーションの種類を「refGene」、「ensGene」から選択できますが、多くの場合「refGene」のほうが無難です。
    先ほど選択したファイルの順番に従って、サンプル名を「input1 sample name」から入力してください。テクニカルリプリケートがある場合には、テクニカルリプリケート内では同じ名前にしてください。また、名前はアルファベットから始まる、アルファベット+数字+「_(アンダーバー)」で名前をつけてください。
    Image
    他のオプションも一通り確認して問題なければ、ページの最下部にある「Run」ボタンを押してください。

Result explanation


    1.結果や進行状況を見るには、プロジェクトページを開いて、「Show Module Flow」をクリックしてください。
    Image
    解析が終了したモジュールや、解析途中のモジュールが出現します。

    2.結果については、主に下記の赤丸1のレポートにまとめられています。また、下記の赤枠2のデータにはTopHat、CuffLinksが出力するいくつかのデータ(ジャンクション情報や、未マップのリードなど)があります。
    Image
    上の図で1番のアイコンをクリックし、解析結果のレポートを表示させます。

    3.まず一番上に、TopHatによってゲノムにマッピングされた結果がGenome Explorerで表示されています。別のウィンドウで開く場合には、「In a new window」をクリックしてください。Genome Explorerの使い方についてはこちらを参照してください。
    Image

    4.下にスクロールしますと、サンプルごとに遺伝子の発現量の分布がヒストグラムとbox plotで表示されます。
    Image
    Image

    5.その次には遺伝子の発現パターンからサンプル間の類似性をクラスタリングした結果が表示されます。テクニカルリプリケート間が最も近い距離にあるか、などを確認できます。
    Image

    6.その次は、二群間で遺伝子発現量について総当たりの散布図、Volcano Plotが表示されます。
    Image
    サンプル間の総当たり散布図

    Image
    サンプル間の総当たりVolcano Plot


    7.次は、いずれかの二群間比較で有意差のあった遺伝子を集めて、発現パターンで遺伝子をクラスタリング(k-means)した結果が表示されます。
    Image

    8.それぞれのクラスターごとに特徴的なGOをGOseqによって抽出した結果へのリンクが表示されます。
    Image
    たとえば、3番のクラスターに特徴的なGOのリストを表示させるには、「cluster: 3」をクリックしてください。すると、GOseqによって遺伝子長のバイアスなどを考慮した上で、有意に出現頻度が増加しているGO、減少しているGOのリストが表示されます。GOseqはバックグラウンドとして、5.のステップで選択した生物種が持っているGOを使用します。そのため、生物種と遺伝子アノテーションの種類によってはGOseqのデータベースに存在しないため、GOseqの結果が存在しません。今のところ、ヒト、マウス、ショウジョウバエでRefGeneアノテーションを選択した場合には、結果を確認しています。ほかの生物種、アノテーションへの対応をご希望の場合には、メールをいただけますと対応を検討いたします。
    Image

    9.各クラスターに特徴的なGOを、よりわかりやすくするために、GOseqの結果をREVIGOに入力します。12.で「p-value < 0.001 Revigo」のボタンをクリックします。
    Image
     REVIGOのウェブサイトが開かれ、クラスター3でp-value 0.001未満のGOリストが入力されている状態であると思います。画面下の方にある「Start Revigo」をクリックしてください。
     この際に、そのクラスターに含まれるコンティグ数が少ないなどの理由で、有意に出現頻度が上昇しているGOが0であることもあります。その際には、1つ戻って、「full Revigo」のほうをクリックすると、すべてのリストが入力されます。
     次の図はREVIGOのTreeMapの図です。
    Image

    10.解析結果レポートの上半分は遺伝子単位で発現量を見たときの結果で、下半分はトランスクリプト単位で発現量を見たときの結果が表示されています。新規のスプライスバリアントを見る場合にはトランスクリプト単位で見てください。ただし、トランスクリプト単位では擬陽性が多くなる傾向があるため、確実性を優先する場合には遺伝子単位で見た方がよいです。

    11.発現量などのテーブルは、レポートページの中央に遺伝子単位の場合、一番下にトランスクリプト単位の結果へのリンクが表示されています。
    Image
     リンクをクリックして、タブ区切りのテキストファイルをダウンロードし、エクセルなどで開いてください。
     ファイルを開くと、一番左の列に「tracking_id」という列があります。これはCuffLinksが予測した遺伝子もしくはトランスクリプトのIDです。そして、これらの遺伝子(もしくはトランスクリプト)が既知の遺伝子である場合には、「gene_short_name」の欄に遺伝子名が表示されます。この遺伝子名に紐付いたGOのリストはファイルの一番右端に記載があります(下記参照)。また、遺伝子(もしくはトランスクリプト)の染色体上の位置や、長さも表示されます。気になる発現パターンの遺伝子については、この染色体上の位置をGenome Explorerで見てみることもできます。
     次にCuffDiffが計算した、各サンプルの発現量が記述されています。FPKMとは1 kbpあたり、100万リードあたりの補正されたリード数で、原理的には異なるサンプル間の比較、異なる遺伝子間の比較ができるように補正された値です(ただし、実際には各種のバイアスにより、FPKMの値をそのまま比較することは難しいケースが多いです)。
    Image
     そして、CuffDiffの発現量のあとには、CummeRbundによってクラスタリングされたクラスターの番号が表示されています。このクラスターの数字は、11.の図と対応しています。その次からは、CuffDiffによる二群間比較解析の総当たり結果が表示されています。
    Image
     その次には、トランスクリプトの配列が表示されます(トランスクリプトのファイルに限定されます)。また、gene idに紐付いたGOのリストも表示されます。
    Image

    12.エクセルの機能で「cluster」番号でフィルタリングすることで、興味のある発現パターンの遺伝子リストを得ることができます。

Workflow


    TopHatによるマッピングを行った後、CuffLinksにて新規の遺伝子予測を行い、既存の遺伝子の発現量とともに発現量を出力する。
    さらにCuffDiffによる二群間比較を行い、その結果をCummeRbundによってグラフ化する。
    また、GOseqによるGO解析を行い、各クラスターで有意に出現頻度が上昇しているGOを抽出し、REViGOにより可視化する。

Comments


  • hg19, mm9, dm3以外でもTophat-Cufflinksまでは、比較的動作する可能性がありますが、GOseqで対応する情報が存在しないため、エラーとなる可能性が高いです。

Use case


  • 各50 bp・8000万ペアエンドリード(ERX011194, ERX11195, ERX11197)3組を入力とし、マッピング、発現量比較、GOアノテーションを行った。


Related information