Platform for Drug Discovery


Trinity, Bowtie, eXpress, DEGseq, clustering and GO analysis (PE)


Introduction


    このパイプラインは、ターゲットのゲノム配列が未知であり、新規にトランスクリプトーム配列を決定するところから解析する場合に用います。Trinityによるde novoトランスクリプトームアセンブルを行った後、bowtieによってリードをコンティグにマッピングし、eXpressによって発現量(FPKM)を算出します。さらにDEGseqによる二群間比較を行い、いずれかの二群間で有意に発現量が変動しているコンティグを抽出し、k-means法によるクラスタリングを行います。また、Uniprotに対するblastx、NRに対するblastnを行い、各コンティグの遺伝子アノテーションを付与します。最後に各クラスターで有意に出現頻度が上昇したGOの一覧を抽出し、Revigoによって可視化します。


    Input formatFASTQ (Illumina only)
    Library layoutFor paired-end read only
    Execution timeAbout a few days (20M paired-end reads[50bp+50bp])


    Results
    • Uniprot, NRアノテーション付き発現量テーブル
    • クラスタリング図 ...etc.
      Image
      クラスタリングを行った結果(HTMLレポート)

Table of contents


Inputs


  • ファイルはペアエンドのFASTQファイルを1~10サンプルまで入力可能です。発現量の比較解析を行わない場合には1サンプルでも動きますが、比較解析以降の結果に意味はありません。
  • FASTQファイルの合計がおよそ100GB以下であれば、経験的には当サーバで解析可能ですが、アセンブルに必要なメモリーを見積もることは難しいため、事前に判定などは行っておりません。
  • スペース「 」や、カッコ「(」、クオーテーション「'」などといった文字や日本語は使えません。アルファベット+数字+アンダーバー「_」で名前を付けてください。
  • ペアエンドの場合、フォワード側は「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については、このページをご覧ください。

Outputs


    実行例1
    Image

    1. 発現量テーブル
    Image
    発現量一覧と、Blastxによるアノテーションリスト(tsvファイル(Excelにて閲覧可能))
    2. クラスタリング図
    Image
    クラスタリングを行った結果(HTMLレポート)
    3. DEGseqによる統計解析 (MA plotなど)
    Image
    アセンブルされたコンティグの発現量を2群間で比較
    4. クラスターごとの特徴的なGOの抽出
    Image
    REVIGOによるGO解析
    5. アセンブル結果
    Image
    FASTAファイル
    6. マッピング結果
    Image
    BAMファイル (samtools viewで可視化)

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)
    詳細はこちらをご覧ください。
    さらにペアエンドのリード名は、リード名中の最初のスペース前までの文字を、フォワード、リバース両方とも同じ名前にして下さい。しばしばフォワード側のリード名には「/1」、リバース側には「/2」という文字が付与されていますが、eXpressが対応していないためペアエンドとしての発現量比較ができない可能性があります。

    また、FASTQファイルのクオリティにはいくつか種類があります。自分のデータがどの種類なのか分からない場合には、FASTQCなどのツールを先に使用して下さい。FASTQCの「Basic Statistics」項目で「Sanger / Illumina 1.8」などと表示されれば特に気にすることはありませんが、もしもIllumina 1.3 - 1.7の場合には、後述のbowtieのオプション設定時に「--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」を選択し、本パイプライン名(Trinity, Bowtie, eXpress, DEGseq, clustering and GO analysis (PE))の横の「Analysis」をクリックしてください。
    Image

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

    12.先ほど確認したファイルの順番に従って、サンプル名を「input1 sample name」から入力してください。
    Image

    13.二群間比較を行う際のp-valueの閾値と、k-meansクラスタリングによって分けられるクラスターの数を設定してください。
    Image
    クラスターの数ですが、サンプル数が3であれば20程度、サンプル数が5であれば60程度が目安です。サンプル数がそれ以上に増えてきますと、k-meansで少数グループに分類することが難しくなるため、パターンマッチングなど別の解析を後で行う必要があるかもしれません。

    他のオプションも一通り確認して問題なければ、ページの最下部にある「Run」ボタンを押してください。

Result explanation


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

    2.まずはアセンブル結果を見てみましょう。
    Image
    上の図で5番のアイコンをクリックし、ダウンロードページを開きます。

    3.次のように3つのファイルがあります。
    Image
    ・Trinity.fasta ・・・ アセンブルされたコンティグファイル
    ・Trinity.fasta.fai ・・・ コンティグファイルのインデックス(末尾に示すIGVなどで見る場合に必要)
    ・summary ・・・ アセンブル結果の統計値、およびマップ率
    summary横の「Https」をクリックし、データを開きます。

    4.アセンブルされたコンティグのコンティグ数、N50 (bp)(アセンブルにおいてよく用いられる平均長のような指標)、最大コンティグ長(bp)などが表示されます。
    Image
    下の方には、アセンブルされたコンティグに対して、bowtieでマッピングしたときのマップ率がinput1、input2と順番に表示されます。

    5.次に発現量のヒストグラムと、二群間比較の総当たり結果を見てみましょう。Outputの3番をクリックして、HTMLを開きます。
    Image
    横軸に発現量、縦軸にはカーネル密度推定による頻度値が表示されます。
    その下には、DEGseqのMATRによって二群間比較が行われた結果のMA plot図の一覧が表示されます。

    6.目的のサンプルで有意に発現量が増加・減少しているコンティグを調べましょう。
     Outputの1番を開き、データをダウンロードしてください。1番の「expression table」はタブ区切りのテキストデータです。エクセルなどで開いてください。横に長いテーブルですが、項目を一つずつ説明しますと、まずはコンティグ名が先頭にあります。
     コンティグ名は、「comp1000_c0_seq1」などとなっております。大まかには、「comp1000_c0」までが同じで、最後のseq1, seq2などが異なるコンティグはスプライシングバリアントの関係に当たります。コンティグ名の詳細はTrinityのマニュアルをご参照ください。
     その次に、コンティグの長さが表示されています。
     さらにその次からは、サンプルごとにeXpressによって算出された発現量がそれぞれ表示されます。eXpressはコンティグ長やリードのバイアスなどを考慮してFPKMを計算しています。解析する際によく使用する値は「fpkm」の項目で、補正された発現量となります。詳細はeXpressのマニュアルをご覧ください。
    Image
     右にスクロールしていきますと、コンティグの配列があります。エクセルで開いた場合の注意事項として、2万bpを超えるコンティグの場合には、配列の途中で改行されてしまいます。
     その次には、DEGseqによるサンプル間の二群間比較の総当たり解析の結果が表示されます。
    Image
     また右にスクロールしていきますと、クラスタリングを行った結果のクラスター番号が表示されています。番号が振られているコンティグは、二群間比較の総当たり解析の際に、どこかの二群間で有意に変動があったコンティグです。
     その次には、Uniprotに対してblastxを行った遺伝子アノテーションの情報が表示されます。UniprotのIDからはGOが豊富に紐づけられているため、その情報も併せて表示されています。
    Image
     さらに右にスクロールしますと、最後にNCBIのNRに対するblastnの結果が表示されます。NRのgiには生物種の情報が豊富に紐づけられているため、その情報も併せて表示されています。「superkingdom」の項目はtaxonomyのトップレベルカテゴリー(external link)のことです。
    Image

    7.ざっくりと発現パターンを分類して、目的の発現パターンに含まれる遺伝子を知ることもできます。(例:1,3日目では発現量が低いが、7,10日目では発現量が高い遺伝子、など)
     Outputの2番を開いてください。
    Image
     この例では、左から順番にサンプル1,2,3として、サンプル3で特異的に発現量が上昇しているクラスターは4番クラスターとなります。先ほどのエクセルファイル中の「cluster」項目に「4」と記入されているコンティグが、このクラスターに含まれるコンティグです。
     エクセルの機能で「cluster」番号でフィルタリングすることで、興味のある発現パターンの遺伝子リストを得ることができます。

    8.また各クラスターごとの大雑把な傾向を把握するため、各クラスターごとに特徴的なGOを可視化することができます。
     Outputの4番を開きます。
    Image
     その中の「cluster: 4」をクリックします。
    Image
     このリストは各クラスター vs バックグラウンド(二群間比較の総当たり解析でどの二群間でも有意に発現量が変化しないコンティグ)のGO比較をFET検定によって行い、各クラスターで出現頻度が上昇、減少しているGOのリストが表示されています。この図から、このサンプル3では「immune response」というキーワードの出現頻度の増加が目立つようにも見えますが、よりわかりやすくするために、REVIGOという外部サイトで可視化を行います。
     一つ戻って、クラスター4の「p-value < 0.001 Revigo」をクリックします。
    Image
     REVIGOのウェブサイトが開かれ、クラスター4でp-value 0.001未満のGOリストが入力されている状態であると思います。画面下の方にある「Start Revigo」をクリックしてください。
     この際に、そのクラスターに含まれるコンティグ数が少ないなどの理由で、有意に出現頻度が上昇しているGOが0であることもあります。その際には、1つ戻って、「full Revigo」のほうをクリックすると、すべてのリストが入力されます。
    Image
     REVIGOの結果を見ますと、クラスター4では「defense response」、「immune system process」に関わる遺伝子が多く、免疫系の防御反応が起きていると推測されます。

Comments

  • メタゲノムのようなサンプルの場合、de Bruijnグラフが複雑になりすぎ、リード数が5000万以下でも、Trinityの途中で計算が終わらなくなる場合があります。
  • 出力されるtsvファイルは、通常10万行以上の長さがあるため、Excel 2007以降でないと開けない可能性が高いです。
  • 少し前のIlluminaから出力されるFASTQは、クオリティスコアの表記が異なるため、"--phred64-quals"のオプション指定を行ったほうが良い場合があります。判断するには、FastQCによるチェックを行い、「Encoding:Sanger/Illumina 1.9」となっていなければ、指定したほうが良い可能性があります。
  • このパイプラインは、Illumina用にパラメータが設定されており、Solid, 454, IonTorrentのデータを処理することは基本的にはできません。特にbowtieでマッピングする箇所が、454、IonTorrentのデータには不向きです。もしも無理に使う場合は、bowtieのオプションを「-a -v 3」から「-k 10 -v 10」などとすると多少は動くかもしれません。
  • アセンブルされたコンティグのマッピングされた様子を確認するには、IGVなどのビューアで確認します。
以下にIGVでMASER上のファイルを直接読み込む方法を紹介します。

1.IGVのサイトを開き、IGVを起動する。
IGVのHP:http://www.broadinstitute.org/igv/home
Image
IGVのJava Web StartのURL:http://www.broadinstitute.org/igv/projects/current/igv.jnlp

2.IGVを開いたら、Genomes -> Load Genome from URLを選択します。
Image

3.MASERに戻り、データのダウンロードページからFASTAファイルのURLをコピーします。
Image

4.IGVにFASTAファイルのURLを貼り付けます。
Image
ユーザ名とパスワードを聞かれると思いますが、MASERにログインするときに使用するユーザ名とパスワードを入力してください。
IGVがリファレンスを読み込むのに数十秒ほどかかります。

5.次にマッピングデータを読み込みます。File -> Load from URLを選択します。
Image

6.先ほどと同様にMASERに戻り、BAMファイルのURLをコピーします。コピーしたURLをIGVに貼り付けます。
Image

7.IGVでデータを見る準備ができました。さらにBAMファイルを追加したり、GTFファイルを読み込ませることもできます。また、BAMファイルのトラック上で右クリックし「view as pairs」を選択するとペアエンドに対応した状態で表示させることなどもできます。詳しくはIGVのマニュアルをご覧ください。
Image

Use case

  • Example1: 各50 bp・8000万ペアエンドリード(ERX011194, ERX11195, ERX11197)3組を入力とし、アセンブル、アノテーションを行った。

Related information