Loading...
 
印刷
日本語

Cuffdiff出力を変換して絞込みに便利にするフロー

目次


説明

  • RNA-seq解析で最もよく使われるツールの一つがTopHatとCufflinksの組み合わせですが、Cufflinksプログラム群の中でサンプル間/群間比較を行い統計的有意な差のある転写物(や遺伝子、CDS、splicing junction)を検出してくれるサブプログラムがCuffdiffです。
  • Cuffdiffの出力するファイルはタブ区切りのテキストファイルなので、TSVファイルとしてMicrosoft Excel他表計算ソフトで開くことも、その応用で複数ファイルを切って繋げてということは、プログラムの書けない方にとっても不可能ではないとは思います。
  • 一方で、実際の解析現場では数サンプル(場合によっては十数サンプル)の比較を同時に行うことも稀ではないので、これらをいちいち手作業で切り貼りするのはヒューマンエラーの原因にもなり、管理が煩雑にもなります。
  • ここで記載するフローはせいぜい2~10サンプル程度を想定し横断的な絞込みができるようなファイルを作製する過程を自動化しています。
  • Cuffdiffの出力ファイルから出発して、いくつかのperlスクリプトで変換や情報の付加を行うことで、多サンプルでの複雑な条件(例えば、サンプル1とサンプル2では有意な発現上昇があり、サンプル1とサンプル3では有意かどうかに関わらず発現の上昇が見られ、染色体11番に由来の転写産物)をプログラムをかけなくとも可能にします。
  • 出力形式はブラウザで開くことで絞込みやソートを可能にしたHTML形式と、Excel等で開くことでExcelのオートフィルタや条件付書式の機能を使って好きなように絞込みや加工ができるCSV形式の2つを用意しています。

処理概要

  • 以下のステップからなります。
#任意主要スクリプト名概要
1-divideMultiCuffdiff.plで分割
mergedCuffdiffToTsvForReport.plでマージ
比較条件を横並びに結合
2任意extractExtraInfFromGtf.plで情報抽出
generalMergerOnmemory.plで結合
アノテーション属性情報の付加
3任意addUCSCandGELink.plゲノムブラウザへのダイレクトリンク付加
4任意addEntrezGeneLink.pl遺伝子アノテーション情報の付加
5-tsvFromCuffdiffOutToHtmlAfterMakeLink.plレポート化

  • 任意のステップは省略しても最終的な出力に付加情報がなくなるだけで問題ありません。

各ステップの説明

比較条件を横並びに結合

  • Cuffdiffに3サンプル以上のマッピング後のファイルを入力として渡すと(グループ化しない限りは)それらの総当り比較が行われ、それらの比較結果が上図のように縦並びで出力されます。
  • この状態では、複数の比較条件を組み合わせた絞込みに不都合なので、横並びに変換します。
  • この際、多数の列を全て横につなげると全体の見通しが悪くなることから、増減("log2(fold_change)"列の値が正か負か)、比較検定の妥当性("status"の列の値)、有意差の有無("significant"列の値)を勘案し、比較の結果を以下の6種の状態にカテゴリ分けしたダイジェスト(図中digest)列を付加することで、複数の比較条件の結果を少ない列数で一覧できるようにしています。
    • #表記概要定義
      1sigUp有意かつ発現増"significant"列の値="yes"かつ"status"列の値="OK"かつ”log2(fold_change)”の値が負
      2up有意でなく発現増"significant"列の値="no"かつ"status"列の値="OK"かつ”log2(fold_change)”の値が負
      3sigDown有意かつ発現減"significant"列の値="yes"かつ"status"列の値="OK"かつ”log2(fold_change)”の値が正
      4down有意でなく発現減"significant"列の値="no"かつ"status"列の値="OK"かつ”log2(fold_change)”の値が正
      5equal発現増減無し"status"列の値="OK"かつ”log2(fold_change)”の値が0
      6notEnough発現量判定不能"status"列の値が"OK"以外

アノテーション属性情報の付加

    • 特にアノテーション情報としてGENCODE由来のgtfを使用する場合、GENCODEアノテーション上で遺伝子のIDは"ENSG00000181449.2"のような番号、転写物のIDも同じく"ENST00000431565.2"のような番号で表されます。結果として出力はこれらIDに対する結果行が出力されます。
    • 多くの研究者は自身が着目する遺伝子を遺伝子名"SOX2"、転写物の"SOX2-201"のようなgeneSymbolの表記で記憶されていると思うので、これらの情報が付加していないと毎回IDを変換するという手間がかかります。
    • GENCODEのアノテーションgtfファイル中にはこれらgeneSymbolの情報が含まれているので、これらをこの処理で対応付けて付加してあげることで、毎回IDの変換を行う必要も無く、geneSymbolを用いて簡易検索したり、結果を絞り込むなどができるようにしています。
    • cuffdiffの出力中にはもともと3列目に"gene"列が存在し、ここに"SOX2"のような遺伝子のgeneSymbolの情報は含まれます。しかし、isoform_exp.diffといった、転写物単位の結果ファイルにも"SOX2-201"のような転写物のIDは含まれないので、これは後付で付加する必要があります。

ゲノムブラウザへのダイレクトリンク付加

  • RNA-seqの解析で特に複雑な転写制御を受けていたり、多数のsplicing variantが存在する座位の発現量をCuffdiffの結果のような統計値だけから判定すると思わぬ誤解をしてしまう場合がありえます。
  • そのような事態を軽減する意味でも、着目した遺伝子等の周辺のマッピング状況や予測転写物のゲノム上の領域をゲノムブラウザで目視確認するのは重要です。
  • ここでは、cuffdiffの出力ファイルにある座標情報から、各種ゲノムブラウザ(ここでは、UCSC genome browserと、当スクリプト開発チームで開発されたゲノムブラウザのGenome Explorerのみ対応)へのダイレクトリンク列を付加する処理を行います。
  • 今回紹介するスクリプトでは対応していませんが、IGVのようなPCにインストールするタイプのゲノムブラウザでも同様の処理を行うことは容易なのでスクリプトを置き換えて利用いただくのも良いかもしれません。

遺伝子アノテーション情報の付加

    • 上記のフローチャートでは、着目した各種アノテーション情報の参照先のポータルサイトとして利用できるEntrez geneへのダイレクトリンク列を付加する処理を示しています。
    • ダイレクトリンクはオプションで指定された列に格納された遺伝子名(この場合はLOCUS列に格納されたgeneSymbolを想定)と、同じくオプションで指定された生物種名をクエリに検索した結果と同一になります。

レポート化

    • タブ区切りテキストファイル(TSVファイル)でもExcel等で開いて、フィルタ機能で絞り込んだりはできますが、WEBブラウザ(Google chrome推奨)でフィルタやソートが可能な形式のレポートに整形を行います。
    • 絞込み等に使われるプログラムはgoogleのサイトから自動的にダウンロードされるため特別プログラムのインストールを行う必要はありませんが、インターネットに接続可能なパソコン等から接続する必要があります。
    • データ自体はインターネット上にアップロードされたりはしないので安心下さい。
    • 特にGENCODEアノテーション情報を利用したcuffdiff結果を整形すると、ファイルが非常に大きくなり、ブラウザでは開けない場合もあります。その場合は、同時に作成される'.csv'ファイルをExcel等で開いて絞込みやソートを行って下さい。
    • Internet Explorerで開くと、警告メッセージが頻発し、いつまでも開けない場合があるのでご留意下さい。
    • html番のレポートの利用方法は以下のページと共通する部分が多いのでそちらをご覧下さい。

実処理

処理前のファイル

  • $ ls -1 ./input_1/
    bias_params.info           isoforms.fpkm_tracking
    cds.count_tracking         isoforms.read_group_tracking
    cds.diff                   promoters.diff
    cds.fpkm_tracking          read_groups.info
    cds.read_group_tracking    run.info
    cds_exp.diff               splicing.diff
    gene_exp.diff              tss_group_exp.diff
    genes.count_tracking       tss_groups.count_tracking
    genes.fpkm_tracking        tss_groups.fpkm_tracking
    genes.read_group_tracking  tss_groups.read_group_tracking
    isoform_exp.diff           var_model.info
    isoforms.count_tracking
    • Cuffdiffプログラムの出力するファイル群は上記のように多数ですが、今回の処理対象は".diff"で終わる以下のファイルのみです。
  • $ ls -1 ./input_1/*.diff
    ./input_1/cds.diff
    ./input_1/cds_exp.diff
    ./input_1/gene_exp.diff
    ./input_1/isoform_exp.diff
    ./input_1/promoters.diff
    ./input_1/splicing.diff
    ./input_1/tss_group_exp.diff
  • $ head -n 5 ./input_1/gene_exp.diff
    test_id gene_id gene    locus   sample_1        sample_2        status  value_1 value_2 log2(fold_change)       test_stat       p_value q_value significant
    XLOC_000001     XLOC_000001     DDX11L1 chr1:11873-29370        q1      q2      NOTEST  0.00500908      0       -inf    0       1       1       no
    XLOC_000002     XLOC_000002     OR4F5   chr1:69090-70008        q1      q2      NOTEST  0       0       0       0       1       1       no
    XLOC_000003     XLOC_000003     LOC100132062,LOC100133331       chr1:323891-328581      q1      q2      NOTEST  0.0100167       0       -inf    0       1       1       no
    XLOC_000004     XLOC_000004     OR4F3   chr1:367658-368597      q1      q2      NOTEST  0       0       0       0       1       1       no
    • 説明中では主に、"gene_exp.diff"を処理したものを中心に説明します。処理自体は他のファイルも行われます。

比較条件を横並びに結合

  • for infile in `ls ./input_1/*.diff`;
    do
      numOfLine=`cat $infile | wc -l`;
      if [ "$numOfLine" -gt 1 ] ; then
      outname=`basename $infile | sed "s/.diff$//"`;
      @CIP_HOME@/bin/processCuffdiff/divideMultiCuffdiff.pl -i $infile -o ./${outname}.;
      numComp=`ls -1 ./${outname}.* | wc -l`;
      list=""; 
      for i in `seq 1 $numComp`;
      do
        list=$list" "./${outname}.$i;
      done;
      paste $list | @CIP_HOME@/bin/processCuffdiff/mergedCuffdiffToTsvForReport.pl #OPTION_1# -o ./output_1/${outname}.tsv;
      fi;
    done;
    • "input_1"ディレクトリに"*.diff"ファイルが置かれている想定です。
    • "output_1"ディレクトリに処理結果の"*.tsv"ファイルが作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。

アノテーション属性情報の付加

  • reffile=`ls -LS1 ./input_2/*.* | egrep "(.gtf)$" | head -n 1`;
    @CIP_HOME@/bin/processCuffdiff/extractExtraInfFromGtf.pl ./gene.inf.tsv ./trans.inf.tsv < $reffile;
    for infile in `ls -LS1 ./input_1/*.* | egrep ".tsv$"`;
    do
      outname=`basename $infile`;
      if [[ "$infile" =~ "gene_exp.tsv" ]] ; then
        awk "{if(NR==1){sub(\"gene_id\",\"test_id\")};print \$0;}" ./gene.inf.tsv | @CIP_HOME@/bin/generalMergerOnmemory.pl --keep_order_of=1 -k 1 $infile - | awk -F "\t" "\$2!=\"\"{print \$0}" > ./output_1/$outname;
      elif [[ "$infile" =~ "isoform_exp.tsv" ]] ; then
        awk "{if(NR==1){sub(\"transcript_id\",\"test_id\")};print \$0;}" ./trans.inf.tsv | @CIP_HOME@/bin/generalMergerOnmemory.pl --keep_order_of=1 -k 1 $infile - | awk -F "\t" "\$2!=\"\"{print \$0}" > ./output_1/$outname;
      else 
        @CIP_HOME@/bin/symlinkIfNotExist.sh $PWD/input_1/$outname ./output_1/$outname;
      fi;
    done;
    • "input_1"ディレクトリに一つ前の処理で作成された"*.tsv"ファイル群が存在する想定です。
    • "input_2"ディレクトリにgtfファイルが置かれている想定です。
    • "output_1"ディレクトリに処理結果の"*.tsv"ファイルが作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。

ゲノムブラウザへのダイレクトリンク付加

  • dataStr=`cat ./input_2/GEDatasetIdList.txt`;
    for infile in `ls -LS1 ./input_1/*.* | egrep "(.tsv)$"`;
    do
      outname=`basename $infile`;
      @CIP_HOME@/bin/mergedExp/addUCSCandGELink.pl -g #OPTION_1# -d $dataStr -i $infile -o ./output_1/$outname;
    done;
    • "#OPTION_1#"は'hg19'などゲノムバージョンの識別子を指定して下さい。
      • 生物種とゲノムIDの選択肢はこちら(external link)のページの"Genome ID"列などをご参考に下さい。
    • "input_1"ディレクトリに一つ前の処理で作成された"*.tsv"ファイル群が存在する想定です。
    • "output_1"ディレクトリに処理結果の"*.tsv"ファイルが作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。
    • 変数"$dataStr"は、GenomeExplorer利用者のみに関係します。データのID('DS00000001'など)をカンマ区切りで指定すると、そのデータトラックが直接開けます。
    • UCSC Genome Browserで直接データを開きたい場合は、予めデータをアップロードトラックを登録しておく必要があります。

遺伝子アノテーション情報の付加

  • for infile in `ls -LS1 ./input_1/*.* | egrep "(.tsv|.exprof)$"`;
    do 
      outname=`basename $infile`;
      @CIP_HOME@/bin/mergedExp/addEntrezGeneLink.pl #OPTION_2# #OPTION_3# -i $infile -o ./output_1/$outname;
    done;
    • "#OPTION_1#"は外部ユーザーが利用する機会は無いと思うので省略。
    • "#OPTION_2#"には以下の選択肢があります。
      • #指定する値説明
        1何も指定しない当該geneSymbolを持つ全ての生物種がリンク先のリストに含まれます。
        2'-s "Homo sapiens" ''-s'と生物種の指定を行う。空白が入る場合はクォート(")で囲む。
        指定された生物種の遺伝子情報のみがリンク先のリストに含まれます。
    • "#OPTION_3#"にはEntrez geneへのクエリに利用する遺伝子のID列の列番号を'-k'オプションで指定します。この場合はcuffdiffの出力のLOCUS列に相当する'-k 3'を指定します。
    • 実際の例
      ./bin/addEntrezGeneLink.pl -s 'Homo sapiens' -k 3 -i indir/gene_exp.tsv -o outdir/gene_exp.tsv;
    • "input_1"ディレクトリに一つ前の処理で作成された"*.tsv"ファイル群が存在する想定です。
    • "output_1"ディレクトリに処理結果の"*.tsv"ファイルが作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。

レポート化

  • mkdir ./output_1/data;
    for infile in `ls -LS1 ./input_1/*.* | egrep "(.tsv)$"`;
    do
      outname=`basename $infile | sed "s/\(.tsv\)//"`;
      @CIP_HOME@/bin/processCuffdiff/tsvFromCuffdiffOutToHtmlAfterMakeLink.pl -i $infile -o ./output_1/data/${outname}.html -c ./output_1/data/${outname}.csv -p ./output_1/data/${outname}.prop;
    done;
    cd ./output_1/; 
    @CIP_HOME@/command/wrapHTMLmaker.pl ./data/ ./report.html; 
    cd ../;
    • "input_1"ディレクトリに一つ前の処理で作成された"*.tsv"ファイル群が存在する想定です。
    • './output_1/data/’ディレクトリ以下に処理結果の'.html'ファイルと'.csv'ファイルが生成されます。
    • wrapHTMLmaker.plは単純にdataディレクトリ以下の'.html'と'.csv'ファイルへのリンクページを作成します。
    • '.prop'という拡張子のファイルはwrapHTMLmaker.plの制御用に用いられる属性情報です。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。


Last edited by nobfujii .
Page last modified on 月曜日 29 of June, 2015 22:45:12 JST.

List of attached files
ID 名前 desc uploaded サイズ ダウンロード数 操作
39 default wrapHTMLmaker.pl 2015-06-29 22:36 by nobfujii 6.61 Kb 494 表示 Download  
38 default tsvFromCuffdiffOutToHtmlAfterMakeLink.pl 2015-03-22 12:35 by nobfujii 33.27 Kb 1452 表示 Download  
37 default extractExtraInfFromGtf.pl 2015-03-22 12:35 by nobfujii 3.98 Kb 2844 表示 Download  
36 default mergedCuffdiffToTsvForReport.pl 2015-03-22 12:35 by nobfujii 8.06 Kb 1146 表示 Download  
35 default divideMultiCuffdiff.pl 2015-03-22 12:34 by nobfujii 2.07 Kb 516 表示 Download  
34 default addEntrezGeneLink.pl 2015-03-22 12:33 by nobfujii 4.31 Kb 490 表示 Download  
33 default addUCSCandGELink.pl 2015-03-22 12:33 by nobfujii 5.52 Kb 859 表示 Download  
32 default generalMergerOnmemory.pl 2015-03-22 12:32 by nobfujii 4.03 Kb 522 表示 Download