Platform for Drug Discovery


[日本語]  

TopHat, HTSeq, DESeq multi samples differentially expressed gene detection analysis (less than 10 N=1 samples)


Introduction


    TopHat, HTSeq, DESeq multi samples differentially expressed gene detection analysis (less than 10 N=1 samples)

    Input formatfastq (paired-end), gtf
    Library layoutFor paired-end read only
    SpeciesUnspecified
    Execution timeAbout 3 days (3 samples: each sample contains 70 million reads)


    Results

    Image Image

Table of contents


Inputs


    入力1~10:fastq (paired-end)

  • ペアエンドのショートリードRNA-seq配列ファイルをfastq(ペアエンド)形式で指定します。
  • 比較するサンプルごとに2~10ペアを指定します。入力3~10は省略可能です。
  • データ書式の説明はこちらをご覧ください。
  • データセット型: fastq <paired-end>

    入力12:Null dataset (This input always must be empty)

  • 内部的に利用される空のデータで、何も指定できませんし、指定する必要がないものです。

Pipeline history


    Image

Outputs

Image
(1) Raw read level mapping status on Genome Explorer
Image
(2) Raw tag count per gene

Image
(3) Tag count matrix

Image
(4) Statistics about reads not assigned any gene

Image
(5) Statistical comparison results raw data

Image
(6) Figures about each comparison

Image
(7) Report of statistical comparison results

(1) Raw read level mapping status on Genome Explorer

  • 当プロジェクトで開発されたゲノムブラウザGenomeExplorerで生リードレベルのマッピング状況を確認できます。
  • 場合によってはミスマッチなどの多型情報も確認できます。詳しい使い方は以下をご覧ください。
  • 以下に全サンプル(この場合は3サンプル)と入力に使用した既知遺伝子情報を同時に閲覧する手順を記載します。
    • パイプラインの実行後にできた履歴画面のge_db(GenomeExplorer DataBase)型のアイコンの下に表示された"Select"ボタンを順にクリックします。
    • "Select"ボタンを押す度に、サブ画面に選択したデータが追加されていきます。
    • サブ画面で実行可能なメニューが表示されるので"Genome Explorer"と記載のあるボタンを押します。
    • 選択されたデータの確認 兼 データの追加・削除用リストが表示されますが、何もせず"Genome Explorer"ボタンを押します。
    • GenomeExplorerが起動しますが、現行の仕様では既知アノテーション情報が見づらい書式で初期表示されるので、表示形式を変更します。
    • 右上の設定ボタンを押し、設定画面の"Plot Style"タブを押します。
    • キーボードのShiftボタンかCtrlボタンを押しながら、既知遺伝子アノテーショントラックの順鎖(+)と逆鎖(-)を選択し、Plot Typeを"Connected Box"に変更します。
    • この操作で、既知遺伝子情報がExon-Intronが区別できる箱型で表示されるようになったので、遺伝子構造が見やすくなりました。
  • Image
    • 上の例は、表示形式変更後のGenome Explorer画面です。
    • デフォルトではRNA-seqパイプラインはdirectionalなRNA-seqを想定し、順鎖方向の発現と逆鎖方向の発現を区別して見られるよう、別トラックとして表示しています。
    • ここで表示されている箇所は、後述するレポートで有意な発現変動があった箇所を表示しています。
    • 以下のような点が見て取れるかと思います。
      • 中央のSTARD4遺伝子が、サンプル1(図中RNA-seq1)に比べてサンプル2(図中RNA-seq2),サンプル3(図中RNA-seq3)で発現量が増加している。
      • 画面表示上のラベルは切れている左側に位置するCAMK4遺伝子がサンプル2のみで高発現している。
      • 画面右のLOC100505678遺伝子がサンプル3で発現上昇している。
    • 尚、この箇所でビューワ上では目立たないLOC100505678は、下の群間比較時はp-valueがとても低いため、ビューワ上の目視と、比較検定は補完的な役割を果たしているとも言えるかと思います。
    • 画面中央にある黒い箱に緑の文字の記述のあるポップアップは、既知遺伝子をマウスでクリックした際に表示されるID情報と白字の領域座標情報です。

(2) Raw tag count per gene

  • Image Image
    ブラウザ(google chrome)で開いたデータ  Microsoft Excelで開いたデータ

  • HTSeq countプログラムが出力する遺伝子ごとに張り付いたタグ(リード)数をカウントした値です。
  • HTSeq countではgtfファイル中のgene_idでグループ化した領域に対してタグカウントをするので、gene_idを共通する複数のスプライシングバリアントの和集合とみなすことができます。
  • ファイルの末尾には、どの遺伝子にも割り当たらなかったリード数が計上されます。
    • 詳しくは NGS Sufer's Wiki(HTSeq count)やHTSeqのマニュアルをご覧ください。
  • テキストファイルなのでブラウザによってはそのまま参照できますが、少々見づらいかと思います。
  • ダウンロードしたテキスト形式のデータファイルをMicrosoft Excelで開いた図の方が見やすいかと思いますが、開く際には以下のような注意が必要です。

(3) Tag count matrix

  • Image Image
    ブラウザ(google chrome)で開いたデータMicrosoft Excelで開いたデータ

  • サンプルごとに分かれていたタグカウント情報を全サンプルでマージしマトリクス化したファイルです。
  • R/Bioconductorなどの各種群間比較プログラムに読み込ませるのに適した形式です。
  • サンプルごとのファイルに存在していた、どの遺伝子にも割り当たらなかったリードに関する行は除外しています。
  • テキストファイルなのでブラウザによってはそのまま参照できますが、少々見づらいかと思います。
  • ダウンロードしたテキスト形式のデータファイルをMicrosoft Excelで開いた図の方が見やすいかと思いますが、開く際には以下のような注意が必要です。

(4) Statistics about reads not assigned any gene

  • Image
    • 上の(3)で除外されたどの遺伝子にも割り当たらなかったリード数の全サンプルマージ結果ファイルです。
  • 上の画面の例はダウンロードしたテキスト形式のデータファイルをMicrosoft Excelで開いています。

(5) Statistical comparison results raw data

  • Image
    • 検定結果を全てマージしたファイルです。閲覧目的よりはより高次の解析パイプラインの入力用途で利用を想定しています。
  • 上の画面の例はダウンロードしたテキスト形式のデータファイルをMicrosoft Excelで開いています。開く際には以下のような注意が必要です。

(6) Figures about each comparison

  • RNA-seqの群間比較を行うDESeqが出力する比較に関係する図です。
  • オプションで指定された比較条件ごとに1列出力されます。比較条件の指定がない場合は、全サンプルをN=1とみなした総当り比較が行われます。
  • #図の名称説明
    1EstimatedDispersionばらつきの推定に関する散布図。縦軸:dispersion 横軸:mean of normalized counts
    2HistgramOfPValuep-valueのヒストグラム。縦軸:Frequency 横軸:P-value
    3MA-plot.DESeqVersionDESeq独自形式のMA-plot。縦軸:Log2 Fold Change 横軸:mean of normalized counts
    4MA-plot.common一般的なMA-plot。縦軸:Log2 Fold Change 横軸:Log2([mean of sample1]x[mean of sample2])/2
    • グラフの意味については詳しくはDESeqのマニュアルをご覧ください。

(7) Report of statistical comparison results

    • gtfファイルによって指定された既知遺伝子領域ごとのサンプル間検定結果と関連情報をマージしたレポートデータです。
    • Internet Explorerでは警告メッセージがでて閲覧できないため、google chromeなどのブラウザでの閲覧を推奨します。
    • ブラウザを通じてメモリを大きく使用するため、古いパソコンなどで閲覧する場合は、他のアプリケーションを閉じるなどの対応をとって下さい。
  • Image
    • 簡易フィルタ機能を含むhtml版と、Excel等のスプレッドシートで開くことのできるcsv版が提供されていますが、本質的にはデータの内容は同一です。
    • 以下のように利点と難点がありますので、状況に応じて使い分けて下さい。
      #ファイル名の例利点難点
      1html版(何某).perEntry.tested.html簡易的なフィルタリングが容易に行えます。検定対象のエントリー数が多いと表示でません。
      Internet Explorer非対応。
      データの加工は困難。
      2csv版(何某).perEntry.tested.csvExcel等のスプレッドシート用アプリケーションで開き複雑な絞り込みや、データの加工ができます。
      データサイズが増えても大きな問題ありません。
      開くまでに手間がかかります。
      3共通点-ダウンロードしてオフラインで利用できます。
      ファイルそのものにデータを含んでいます。
    • 列の説明
      #列名説明
      1IDエントリーの識別子。通常はgtfファイル中のgene_idで指定された文字列が利用されます。
      2(A群名)/ (群B名): UpDownサンプルBに対するサンプルAの増減および、検定結果を組み合わせた記号。詳しくは下表Up/Down参照。
      3(サンプル名) (tag count)サンプルの生タグカウント(補正前のオリジナル集計値)。
      4(A群名)/ (群B名): baseMeanA群とB群の比較におけるA群、B群全てサンプルの正規化されたカウントの平均。
      5(A群名)/ (群B名): baseMeanAA群とB群の比較におけるA群の正規化されたカウントの平均。
      6(A群名)/ (群B名): baseMeanBA群とB群の比較におけるB群の正規化されたカウントの平均。
      7(A群名)/ (群B名): foldChangeA群からB群への正規化されたカウントのFold Change値。baseMeanBをbaseMeanAで割ったものに等しいので、B群が2倍なら2、B群が1/2倍なら0.5になります。
      8(A群名)/ (群B名): log2FoldChangeA群からB群への正規化されたカウントのFold Change値の2底のlogをとったもの。B群が2倍ならlog2(2)=1、B群が1/2倍ならlog2(1/2)=-1になります。
      9(A群名)/ (群B名): pvalA群とB群のDESeqによる検定結果のp値。検定手法の詳細についてはDESeqのマニュアルをご覧ください。
      10(A群名)/ (群B名): padj検定結果のBH法による補正後のq値。詳しくはR関数のp.adjustをご覧ください。
      11(A群名)/ (群B名): significanceTRUE/FALSE/NA。検定結果の符号。NAは検定できなかった場合。
      12chrom染色体名。
      13strPosゲノム中の領域の開始座標。
      14endPosゲノム中の領域の終了座標。
      15strandゲノムの順鎖(+)・逆鎖(-)のどちらにマップされた情報の集計かを示します。区別しない場合は"*"。
      16GE linkGenome Explorerへのリンク。
      17UCSC linkUCSC Genome Browserへのリンク。
      18Entrez gene linkEntrez Geneへのリンク。ID列と、オプションで指定した生物種をキーに問い合わせた結果のページにリンクしています。
    • UpDownの説明
      equalbaseMeanが、A群とB群で一致する場合
      upbaseMeanが、サンプルBに対してサンプルAが大きい場合かつ検定結果が有意でない場合
      sigUpupかつ検定結果が有意の場合
      downbaseMeanが、サンプルBに対してサンプルAが小さい場合かつ検定結果が有意でない場合
      sigDowndownかつ検定結果が有意の場合
      noData両群のbaseMeanが0の場合。
  • GenomeExplorerへのリンク
    • html版からのGEへの遷移
      csv版からのGEへの遷移
    • html版、csv版どちらも、右端近くの"GE link"列からGenome Explorerに遷移することができます。
    • html版ではwebページ上の"Link"ボタンを押せば別画面でサンプル選択済みのGenomeExplorerが起動します。
    • csv版ではLinkマスにあるURLをコピーしてブラウザで開くとURL中のパラメータで指定されたサンプルが選択された状態でGenomeExplorerが起動します。
    • 橙線と、緑線矢印にあるように、リンクのあった行が表示されていることが確認出来るかと思います。
    • 上で記載したGenomeExplorerの参照方法のときのように全サンプルのアイコンを選ぶ必要はありません。
    • Genome Explorerでの参照方法は、こちらをご参照ください。
    • 現在のバージョンでは、GenomeExplorerへの遷移の度に、初期表示に戻ってしまうため、複数の領域を逐次的に見たい場合は、GenomeExplorerのsearch機能で、場所だけ移動すると便利です。
  • 開閉式の簡易フィルタと列説明
      • html版では、画面の左上開閉式の簡易フィルタの指定箇所(column filter)と、列の説明文(legend)があり、マウスでクリックするたびに開閉します。
  • 簡易フィルタとソート
    • 簡易フィルタ
      ソート
    • 画面左上の"column filter"をクリックすると開かれる選択ボックスやスライダーを利用すると、結果列をフィルタリングすることができます。
    • フィルタの例
      • ID(遺伝氏名など)に対する部分一致フィルタ。
      • UpDown列のsigUp,sigDownを選択して有意なもののみにフィルタ。あるいは発現増のみ/減のみフィルタ。
      • tag count列のスライダーで特定サンプルが10タグ以上のものにフィルタ。
      • log2FoldChangeが2以上(=4倍増)のものかつ、検定結果が有意なもののみフィルタ。
      • 染色体番号11番のみにフィルタ
    • html版の結果の任意のカラムラベルをクリックすると、その列で結果のソートが出来ます。
      • 図中の例では、タグ数による降順(灰色の下三角が表示)と、p値による昇順(灰色の上三角が表示)を記載しています。

How to run this pipeline

  • 解析実行の流れの詳細はこちらをご確認ください。
  • 実行するプロジェクト内に入力となるRNA-seq配列ファイル(2~10個)と、アノテーションファイルをアップロードしておきます。
  • アップロードしたファイルのデータアイコンの下の"Select"ボタンを押した後、"Analysis"ボタンを押します。
  • 選択されたファイルのデータ型に応じた解析の選択肢がリストされるので、その左のメニューから"RNA-seq"を選び、このパイプラインを選びます。
  • 選択された入力ファイルの確認と、プロジェクト(ワークスペース)に誤りが無いかを確認します。
    • サンプルの順番が違っていた場合などは、データの欄をドラッグアンドドロップで置き換えができます。
  • ヒトのデータであれば、そのままデフォルトオプションで"Run"ボタンを押し実行することもできますが、以下のような場合は、"Set option and run"を押して、オプション選択画面に遷移します。
    • ヒト以外の生物で実行する場合
    • レポートやFig中のサンプル名や比較群のラベルを指定したい場合(デフォルトだとData1,Data2..になってしまいます)
    • directional RNA-seqでない場合
    • Read1順方向、Read2逆方向でライブラリ調整していない場合
    • N=1の1サンプル1群以外の比較条件で比較したい場合
    • その他、解析のカスタマイズを行いたい場合。
  • オプション選択についての詳細は下のOptionsの項目をご覧ください。

Options


  • 設定可能な実行オプション値と説明は以下にあります。
  • 特に後半のオプションを変更される場合は、機能ごとに実行可能な単品パイプラインで動作確認を行った後に利用した方が、無駄な待ち時間を減らせる可能性が高いので、検討下さい。
  • #オプション番号オプション名重要性説明
    1パイプライン共通tophat library typeRNA-seqサンプル調整時の方向性の考慮の有無を指定します。厳密にはdirectional RAN-seqの場合は、サンプル調整方法に応じて、fr-firststrandかfr-secondstrandのどちらかを選択し、non-directionalの場合はfr-unstrandedを選択するべきですが、分らない場合はどちらのケースでもデフォルトのfr-unstrandedで大きな問題は生じないケースが殆どです。詳しくはTopHatのマニュアルを参照ください。
    2パイプライン共通tophat2 version-スプライシングサイトを考慮してゲノムにRNA-seqリードをマップするTopHat2のバージョンを指定します。
    3パイプライン共通bowtie version 1 or 2 for tophat-TopHatから呼び出されるbowtieのバージョンを指定します。
    4パイプライン共通samtools version-HTSeqから呼び出されるsamtoolsのバージョンを指定します。
    5パイプライン共通GENOME_IDゲノムのバージョンで、ヒト以外の生物種の配列を指定する場合は必須です。選択肢にない生物はご相談ください。
    6パイプライン共通tophat pair insert size -rペアエンドのRNA-seqのペア間の距離を指定します。省略可能です。Insert size selectionを行っている場合は指定したほうが精度が向上する可能性があります。
    7パイプライン共通TopHat2 other options-TopHatのその他のオプション指定欄です。マップ率が悪いなどチューニングしたい場合に利用します。
    8パイプライン共通Load directionalityゲノムビューワに配列の生データをロードする時の順鎖・逆鎖の区別をするかしないかです。directionalなRNA-seqの場合は、directionalを指定すると順方向と逆方向のゲノムの鎖へのそれぞれのマップ状況が区別されてビューワ上で視認できます。non-directionalなRNA-seqの場合は、順鎖と逆鎖を分ける意味はないので、non-directionalを指定します。
    9パイプライン共通Consider not proper pair reads as two single end resds-マッピングプログラムにも依存しますが、ペアエンドのリードのマッピングに、通常ありえない距離を隔ててマップするケースがあり、このようなケースを"not proper pair"として区別可能にして出力されます。ゲノムビューワ上で、このようなペアリードの間を線で繋いでしまうと、可視化の弊害になるため、デフォルトではこのような"not proper pair"のペアは異なる2つのみなしごリードとみなしてマップしています。デフォルトでオンになっているので、この機能をはずしたい場合はオフにします。
    10パイプライン共通Load directional second as reverse-ペアエンドのサンプル調整は通常リード1が順方向、リード2が逆方向に読まれるため、directionalなRNA-seqでは、リード1は順鎖、リード2が逆方向の場合は順鎖にマップする方が、オリジナルのRNAの方向を適切に表現できます。デフォルトでは、リード2の方向をリード1と逆にする処理がオンなので、あえてオフにしたい場合はチェックをはずします。
    11パイプライン共通HTSeq version-HTSeqのバージョンを指定します。HTSeq countは、発現情報に使用する遺伝子ごとのタグカウント情報を、マッピング結果を入力とし、算出するプログラムです。
    12パイプライン共通HTSeq count mode主に、遺伝子情報の境界領域にマップされたリードや、複数の遺伝子領域にまたがるリードの扱いをどうするかを指定するモードを指定します。詳しくはHTSeqのマニュアルをご参照ください。
    13パイプライン共通HTSeq count strandHTSeqでの発現量集計時の方向性の考慮の有無を指定します。non-directional RNA-seqの場合は"no"を指定します。
    14パイプライン共通HTSeq count other options-HTSeqのその他のオプション値を指定します。
    151.1mergeExpressionFiles column labelsカンマ区切りでRNA-seqのサンプルそれぞれのラベルに利用する文字列を指定します。3サンプルであれば、"Sample1,Sample2,Sample3"のように指定します。デフォルトでは"Data1","Data2",..,"DataN"が利用されます。
    161.2Common option for output file prefix発現情報を全サンプルでマージしたファイルのファイル名のプレフィックスを指定します。デフォルトでは、"merged"が使用されます。
    172.1R version-DESeqを動作する際のRのバージョンを指定します。
    182.2DESeqController control stringDESeqでのサンプル比較条件を指定します。省略時は全サンプルをデータ数N=1の異なる群とみなし、全群総当りで比較を行います。この文字列を適切に指定すればN=1以外の比較など、任意の組み合わせの比較だけ行うことも可能です。こちらの編集は失敗しやすいので、失敗した際の手戻りを減らすため、初期状態で実行した結果が得られた後、"DESecController"単体パイプラインで再実行したほうが待ち時間を減らせますのでご検討下さい。記載の方法はこちらをご参照ください。
    192.3DESeq method-群間比較を行うDESeqのmethodを指定します。biological replicateのN=1の場合は"blind"、N>1の場合は明示的な指定なしとして機能する"DEFAULT"かその他の指定を検討します。詳しくはDESeqのマニュアルをご参照ください。
    202.3DESeq sharing_mode-群間比較を行うDESeqのsharing modeを指定します。biological replicateのN=1の場合は"fit-only"、N>1の場合は明示的な指定なしとして機能する"DEFAULT"かその他の指定を検討します。詳しくはDESeqのマニュアルをご参照ください。
    212.5DESeq fitType-群間比較を行うDESeqのfitTypeを指定します。biological replicateのN=1の場合は"local"、N>1の場合は明示的な指定なしとして機能する"DEFAULT"かその他の指定を検討します。詳しくはDESeqのマニュアルをご参照ください。
    222.6DESeqController --fig_size-DESeqの比較検定時に出力される比較条件ごとのMA-plotや分布などの図のサイズを指定します。初期値は200ピクセルで、最大1000ピクセルまで指定できます。
    232.7DESeqController --depth_col_rm-列名に総マップリード数の情報を含む場合にそれを除去するかを指定するオプションですが、HTSeqで発現量をカウントした場合は含まれないためチェックする必要はありません。
    243.1Gene ID type-GTFファイルで入力された既知遺伝子情報に対する、ゲノムビューワ上でラベルの種類を指定します。デフォルト値の"AUTO"は"c:Use a composition ID made of gene symbol, accession number and a branching number delimited by ___"と同義で、GTFファイル中のgene_idとtranscript_id、通し番号の3つを"__"区切りで連結した文字を、ラベルに表示可能にします。gene_idとtranscript_idが同一の場合などは"g:Use gene symbol as ID"などを指定した方が表示がシンプルになります。以下の緑文字の表記に利用されます。
    253.2Load directionalityゲノムビューワに既知遺伝子アノテーション情報をロードする時の順鎖・逆鎖の区別をするかしないかです。通常、既知遺伝子情報は順鎖方向と逆鎖方向の区別があるので、"directional"を指定します。
    264.1Organism検定結果レポート中のEntrez Gene IDのリンクを付加する際、どの生物種のIDかを指定するのに利用します。省略化ですが、既知遺伝子情報GTFファイル中のgene_idに利用されることの多い、geneSymbolは複数の生物種で共通なので、生物種を指定しないと、リンク先にすべての生物種の当該遺伝子の詳細がリストされてしまいます。ヒトであれば学名の'Homo sapiens'と指定しておくと、直接ヒトのエントリーのみがリンクから開けます。入力が省かれた場合は、リンクをクリックすると全生物のリストが表示されることになるでしょう。
    274.2MergedDESeqOutTohtml.pl suppress_tag_count-レポートにタグカウントの列を含めない場合にチェックをオンにします。
    284.3MergedDESeqOutTohtml.pl suppress_rpkm-レポートにrpkm値の列を含めない場合にチェックをオンにします。このケースではrpkm値は入らないので、チェックしても意味がありません。
    294.4MergedDESeqOutTohtml.pl suppress_rpm-レポートにrpm値の列を含めない場合にチェックをオンにします。このケースではrpm値は入らないので、チェックしても意味がありません。
    304.5MergedDESeqOutTohtml.pl suppress_fpkm-レポートにfpkm値の列を含めない場合にチェックをオンにします。このケースではfpkm値は入らないので、チェックしても意味がありません。
    314.6MergedDESeqOutTohtml.pl suppress_compare_detail-レポートに比較の詳細情報の列を含めない場合にチェックをオンにします。比較の詳細情報とは、baseMean,baseMeanA,baseMeanB,foldChange,log2FoldChange,pval,padj,significance列です。
    324.7MergedDESeqOutTohtml.pl suppress_position-レポートにゲノム上の座標に関する列を含めない場合にチェックをオンにします。ゲノム上の座標に関する列とはchrom,strPos,endPos,strandです。
    334.8MergedDESeqOutTohtml.pl suppress_links-レポートにリンク列を含めない場合にチェックをオンにします。
    344.9MergedDESeqOutTohtml.pl suppress_others-レポートにその他の列を含めない場合にチェックをオンにします。初期状態ではその他の列はないので、チェックする意味はありません。
    354.10MergedDESeqOutTohtml.pl infinity_value-htmlレポート作成の際に、スライダーによる数値フィルタを有効にするためには、無限大の値が含まれていてはエラーになるため、無限大を特定の大きい数値に置き換えています。このオプションでは置き換える数値を指定します。デフォルトでは9999に置き換わります。

DESeqController control stringの設定方法

  • DESeqController control stringは、DESeqの比較検定を複数回任意の組み合わせで実行するために使用されます。
  • このcontrol stringはR/BioconductorパッケージのDESeqとは独立に私達が開発した機能なので、DESeqのマニュアルなどには、説明がありません。
  • この例で記載している、3サンプルをN=1のそれぞれ3群とみなした3群総当り比較を実施していますが、内部的には自動的に以下のcontrol stringが生成されを用いて実行されます。
    • {2:3:Data1:Data2:q0.1},{2:4:Data1:Data3:q0.1},{3:4:Data2:Data3:q0.1}
    • 上記の"{}"で囲まれた範囲が1比較条件を意味し、それを","区切りで指定しています。
    • ":"で区切られた文字の意味は以下になります。
      #要素の意味
      1A群の列番号
      2B群の列番号
      3A群の群ラベル
      4B群の群ラベル
      5検定の閾値
    • A/B群の列番号は、1列目がID列になるので、2から順番に指定可能です。この場合は3サンプルなので、2~4が指定できます。
    • A/B群の群ラベルは、レポートの列名やMA-plot等で使用されるラベルに使用されるそれぞれの群のラベルです。
    • 検定の閾値は0.01であれば、p値<0.01が有意(significant=TRUE)となり、p1e-05はp値<1e-05が有意、q0.05はFDR補正後のq値<0.05が有意となります。初期値はq0.1です。
    • このため、{2:3:Data1:Data2:q0.1}は、群Aのデータを2列目(サンプル1)、群Bのデータを3列目(サンプル2)とし、群Aのラベルを"Data1"、群Bのラベルを"Data2"とし、閾値はq値<0.1に設定されています。
    • biological replicateを設定して検定したい場合は、以下のように記述すれば2と3をA群、4と5をB群とみなして検定できます。
      {2,3:4,5:groupA:groupB:p0.05}
    • technical replicateで、内部的に複数サンプルをマージして使用したい場合は、以下のように記述します。(この場合は2,4,6がマージ、3,5,7がマージ)
      {{2,4,6}:{3,5,7}:g1:g2:q0.01}
    • もっと複雑な条件で検定を行いたい場合(methodやfitTypeを条件ごとに変えたい)場合は、実行履歴にあるRスクリプトを参考にご自身でDESeqを起動して実行することをお勧めします。
      • 上は実行履歴のDESeqController.plモジュールから自動生成され実行されたLinuxコマンドやRスクリプトの置き場所の開き方。
      • Image
      • 上は1_C000000062/log/controlStringLog.txtファイルの中身で、自動生成されたcontrol stringの内容。
      • Image
      • 上は1_C000000062/log/exec.Rファイルの中身で、control stringで指定された条件を元に、DESeqの複数回検定を実行するためのRスクリプト。
      • 実行ログのその他のファイルの意味についてはこちらをご覧下さい。

Comments



Use case

  • 上記で比較に使用した例は、ENCODE ProjectのRNA-seqデータ(SRA039973(external link))で、多数あるサンプルの内、説明用に3サンプル抜き出したものです。
    • #SRX番号サンプル名含まれるレーン番号
      1SRX084666GSM765388: CshlLong_RnaSeq_MCF-7_cell_longPolyASRR315301,SRR315302
      2SRX082565GSM758559: CshlLong_RnaSeq_GM12878_cell_longPolyASRR307897,SRR307898
      3SRX084683GSM765405: CshlLong_RnaSeq_K562_cell_longPolyASRR315336,SRR315337
  • 今回は説明の簡略化のために、1サンプル2レーン分あるものをあらかじめマージしたfastqから開始していますが、実際には片方が実験的その他の理由で失敗している場合もありえるので、マッピングをそれぞれ独立に実施し、検定時にTechnical replicateをマージして実行したほうが好ましいかもしれません。
  • ただし、このような長いパイプラインの後半の検定でマージが失敗すると、このマニュアルに記載のない単体パイプラインで途中から実行しなおすなど、多少難易度が上がりますので一概にどちらが良いかは言い難い側面もあります。
  • 慣れてきたら、ステップごとに実行する方が、実験的なフィードバックや、解析のチューニングもし易いので、ご検討ください。基本的には、履歴のモジュール名が、単体パイプライン名になっています。

Related information