Loading...
 
印刷
日本語

Cuffdiff出力にChIP-seqの結果を付加するフロー

目次


説明

  • ChIP-seqの結果で予測された転写因子等の結合部位がどの遺伝子付近に存在するかの情報を付加することは、ChIPpeakAnnoなど良く知られたツールが存在します。
  • 一方でRNA-seqの複数サンプルの解析結果と併せてこれらの情報を利用しようとすると、Rプログラミング等のプログラミングスキルがないと利用が困難な状況も続いてきました。
  • ここでは、ChIP-seqの結果に近傍遺伝子情報を付加するのに良く知られたChIPpeakAnnoというR/Bioconductorのパッケージの出力結果を遺伝子単位に再編集し、それを同じくRNA-seqの解析でよく利用されるCufflinks/Cuffdiffの出力結果を整形したものに付加するスクリプト群を紹介します。
  • Image
  • Image
  • ここでは、RNA-seqの結果とChIP-seqの結果を統合的に扱うために遺伝子アノテーション情報を両方で同一のものを利用するノウハウ、および、いくつかのピークコーラー(MACS, PeakSeq, SISSRs)の出力書式をこれらのインプットに利用するための整形手段についても記載します。

処理概要

  • 以下のステップからなります。
    • 前処理(MACSの例)
      #主要スクリプト名
      1bowtie等でマッピング
      2samtoolsでbamに変換
      3MACS2でピークコール
      4modifyBedForChIPpeakAnno.plで
      ChIPpeakAnno用整形
      5ChIPpeakAnnoで
      遺伝子アノテーション付加
    • 前処理(PeakSeqの例)
      #主要スクリプト名
      1bowtie等でマッピング
      2PeakSeq -preprocessで
      ピークコール前処理
      3PeakSeq -peak_selectで
      ピークコール
      4modifyBedForChIPpeakAnno.plで
      ChIPpeakAnno用整形
      5ChIPpeakAnnoで
      遺伝子アノテーション付加
    • 前処理(SISSRsの例)
      #主要スクリプト名
      1bowtie等でマッピング
      2BEDTools bamToBedで
      マッピング結果書式変換
      3SISSRsで
      ピークコール
      4簡易スクリプトで
      独自書式結果をBED変換
      5modifyBedForChIPpeakAnno.plで
      ChIPpeakAnno用整形
      6ChIPpeakAnnoで
      遺伝子アノテーション付加
    • 本処理
      #任意主要スクリプト名概要
      1-divideMultiCuffdiff.plで分割
      mergedCuffdiffToTsvForReport.plでマージ
      比較条件を横並びに結合
      2任意extractExtraInfFromGtf.plで情報抽出
      generalMergerOnmemory.plで結合
      アノテーション属性情報の付加
      3-convertChIPpeakAnnoForIntegration.plで
      ChIP-seq結合部位情報を近傍遺伝子単位に再構築
      ChIP-seq結合部位情報を近傍遺伝子単位に再構築
      4-AddChIPpeakAnnotation for cuffdiffで
      Cuffdiff結果への複数のChIP-seq結果を結合
      Cuffdiff結果へのChIP-seq結果結合
      5任意addUCSCandGELink.plゲノムブラウザへのダイレクトリンク付加
      6任意addEntrezGeneLink.pl遺伝子アノテーション情報の付加
      7-tsvFromCuffdiffOutToHtmlAfterMakeLink.plレポート化

各ステップの説明

比較条件を横並びに結合

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

ChIP-seq結合部位情報を近傍遺伝子単位に再構築

  • Image
    • ChIP-seqの解析結果で得られたピークのリスト(=転写因子等の予測結合領域)に近傍遺伝子情報が付加された書式は、RNA-seqの結果の遺伝子単位の情報と結合しにくい形式のため、書式変換して遺伝子単位の情報に再構成します。
    • オリジナルの書式では、ピーク一つあたり、複数の遺伝子が一定範囲内の近傍にある場合は、複数行になって出力されます。
    • 変換後の書式では、1遺伝子が1行になるように変換されます。この際、どのピークとの関連づかない"Gene B"は出力されていません。また、複数のピークが近傍に存在する"Gene A"では、マス内で2個分のピーク情報が出力されています。
    • 実際の出力書式は上記よりもう少々複雑になります。

Cuffdiff結果へのChIP-seq結果結合

  • Image
    • 遺伝子単位に再構築済みのChIP-seq解析結果を、RNA-seqの整形済みの結果に結合します。
    • 結合時には当然、ピークの無かった遺伝子は空になります。

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

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

レポート化

実処理

処理前のファイル

RNA-seqの結果であるcuffdiff出力ファイル

ChIP-seqの結果であるアノテーション付加済みbedファイル

比較条件を横並びに結合

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

ChIP-seq結合部位情報を近傍遺伝子単位に再構築

  • infile=`ls -LS1 ./input_1/*.* | egrep "(.bed|.txt)$" | head -n 1`;
    outfile=./output_1/`basename $infile | sed "s/\(.bed\|.txt\)//"`.peakAnno.txt;
    @CIP_HOME@/bin/processSeqIntegration/convertChIPpeakAnnoForIntegration.pl -i $infile -o $outfile --id_process #OPTION_1# --out_type #OPTION_2#;
    • "#OPTION_1#"は付加した遺伝子アノテーションのIDを結合前に加工するか否かを意味します。
      • #指定するオプション値説明
        1USE_ORIGINALIDの加工をしない。
        2USE_GENE_ID複合ID(gene_id___transcript_id___branching_number)からgene_idを取り出す。
        3USE_TRANS_ID複合ID(gene_id___transcript_id___branching_number)からtranscript_idを取り出す。
    • "#OPTION_2#"は付加するChIP-seq情報を簡易情報のみに留めるのか、詳細情報を含めるのかを指定します。
      • #指定するオプション値説明
        1DIGEST_ONLY簡易情報のみ含める
        2DETAIL_ONLY詳細情報のみ含める
        3BOTH両方(簡易情報と詳細情報)含める
    • "input_1"ディレクトリに遺伝子アノテーション付加済みの"*.bed"ファイルが存在する想定です。
    • "output_1"ディレクトリに処理結果の遺伝子単位で再構築されたChIP-seq結果ファイル(*.peakAnno.txt)が作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。

Cuffdiff結果へのChIP-seq結果結合

  • labels_org="#OPTION_1#,";
    for infile in `ls -LS1 ./input_1/*.* | egrep ".tsv$"`; 
    do
      labels=$labels_org;
      outname=`basename $infile`;
      echo "Now processing "$outname;
      @CIP_HOME@/bin/symlinkIfNotExist.sh $PWD/input_1/$outname ./output_1/$outname;
      for i in `seq 2 11`;
      do
        if [ -f "./input_${i}/.empty" ] ; then continue; fi;
        target=`ls -LS1 ./input_${i}/*.* | egrep "(.peakAnno.txt)$" | head -n 1`;
        label=`echo $labels | sed "s/^\([^,]\+\),.*$/\1/"`;
        label=`echo $label | sed "s/,//g;s/^ \+//;s/ \+$//"`;
        labels=`echo $labels | sed "s/^[^,]*,//"`;
        if [ "$label" == "" ] ; then
          label="Data"`expr "$i" - 1`;
        fi;
        mv ./output_1/$outname ./;
        header="id";
        if [ "#OPTION_2#" == "BOTH" -o "#OPTION_2#" == "DIGEST_ONLY" ] ; then header=$header","$label" peak digest"; fi;
        if [ "#OPTION_2#" == "BOTH" -o "#OPTION_2#" == "DETAIL_ONLY" ] ; then header=$header","$label" peak detail"; fi;
        echo -e $header;
        @CIP_HOME@/bin/processSeqIntegration/mergeWithMultiId.pl -i ./$outname -o ./output_1/$outname -t $target -k 3 --target_header "$header" --sort_unique_multi " " --multi_id_with_separator ",";
      done;
    done;
    • "#OPTION_1#"はChIP-seqの結果情報行のラベルをカンマ区切りで指定します。
      • 例:'MACS,PeakSeq,SISSRs'
      • 省略された場合は"Data1","Data2",..,"DataN"が利用されます。
    • "#OPTION_2#"は付加するChIP-seq情報を簡易情報のみに留めるのか、詳細情報を含めるのかを指定します。
      • #指定するオプション値説明
        1DIGEST_ONLY簡易情報のみ含める
        2DETAIL_ONLY詳細情報のみ含める
        3BOTH両方(簡易情報と詳細情報)含める
    • "input_1"ディレクトリに二つ前の処理で作成されたアノテーション属性情報の付加の結果の"*.tsv"ファイル群が存在する想定です。
    • "input_2"~"input_11"ディレクトリにそれぞれ、一つ前の処理で作成された遺伝子単位で再構築されたChIP-seq結果ファイル(*.peakAnno.txt)が一つずつある想定です。
    • "output_1"ディレクトリに処理結果の遺伝子単位で再構築されたChIP-seq結果ファイル(*.peakAnno.txt)が作られます。出力用のディレクトリは予め作成して置いて下さい。
    • "@CIP_HOME@"と記載のある箇所は適当にスクリプトを置いた場所に置き換えて下さい。

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

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

レポート化

前処理(ChIP-seq結果ファイルの準備)

処理前のファイル

  • ChIP-seqの解析では、bowtieやbowtie2はたまたBWAなどでマッピング結果のbamファイルを作成するところまでは共通で実行している前提です。
  • 別のマッピング結果書式を利用することもできますが、ゲノムブラウザ等での可視化との両立を考えるとindex付加済みの染色体並び順ソート済みのbamを利用することをお勧めします。

ピークコール by MACS2

  • infile1=`ls -LS1 ./input_1/*.* | egrep "(.bam)$" | head -n 1`; 
    infile2=`ls -LS1 ./input_2/*.* | egrep "(.bam)$" | head -n 1`; 
    outname=`basename $infile1 | sed "s/\.\(bam\)$//"`;
    if [ "#OPTION_2#" != "" -a "#OPTION_2#" != " " ] ; then 
    outname="#OPTION_2#"; 
    fi;
    outname=`echo $outname | sed "s/ /_/g"`; 
    outDir=./log;
    cd $outDir;
    export PATH=@CIP_HOME@/tool/#OPTION_1#/bin:@CIP_HOME@/tool/Python-2.7.3/bin:@CIP_HOME@/tool/numpy-1.6.2/bin:${PATH};
    export PYTHONPATH=@CIP_HOME@/tool/numpy-1.6.2/lib/python2.7/site-packages:@CIP_HOME@/tool/scipy-0.11.0/lib/python2.7/site-packages:@CIP_HOME@/tool/#OPTION_1#/lib/python2.7/site-packages;
    @CIP_HOME@/tool/#OPTION_1#/bin/macs2 callpeak #OPTION_3# -t ../${infile1} -c ../${infile2} --name=${outname} --format=#OPTION_4#;
    if [ -f ${outname}_model.r ] ; then
    R --vanilla < ${outname}_model.r;
    fi;
    cd ../;
    mv ./log/${outname}_peaks.bed output_1/;
    mv ./log/${outname}_summits.bed output_2/;
    
    mv ./log/${outname}_peaks.encodePeak output_4/;
    sed "s/\t-/\t -/" ./log/${outname}_peaks.xls > ./output_3/${outname}_peaks.xls;
    
    grep "^#" ./output_3/${outname}_peaks.xls > ./output_3/${outname}_summary.txt;
    grep -v "^#" ./output_3/${outname}_peaks.xls > ./output_3/${outname}_results.txt;
    
    if [ -f ./log/${outname}_model.pdf ] ; then
    mv ./log/${outname}_model.pdf output_3/;
    convert -density 600 -geometry 1000 ./output_3/${outname}_model.pdf ./output_3/${outname}_model.jpg;
    fi;
    
    cd ./output_3;
    ../myscript ${outname}_summary.txt ./${outname}_peaks.xls ${outname}_results.txt ./${outname}_model.pdf ${outname}_model-0.jpg ${outname}_model-1.jpg > ./index.html ;
    cd ../;
  • myscriptの内容
    #!/bin/sh
    echo "<HTML>";
    echo "<BODY>";
    echo "<title>MACS2 report</title>"; 
    echo "<h2>Analysis settings and information</h2>";
    echo "<object data=\"$1\" type="text/html" border="1" width="800" height="320"></object><br>"; 
    echo ""; 
    if [ -f $4 ] ; then
    echo "<h2>Paired-peaks model</h2>"; 
    echo "<a href=\"$4\" target="new">Download PDF</a><br><br>"; 
    echo "<img src=$5 width="300"><br>"; 
    echo "<img src=$6 width="300"><br>"; 
    echo ""; 
    fi;
    echo "<h2>Peak locations</h2>"; 
    echo "<a href=\"$2\" target=\"new\">Download Excel file</a><br><br>"; 
    echo "<object data=\"$3\" type=\"text/html\" border=\"1\" width=\"800\" height=\"250\"></object><br>"; 
    echo "</BODY>"; 
    echo "</HTML>";

ピークコール by PeakSeq

前処理

  • PeakSeqはbamファイルを入力書式として受け取らないので、独自の前処理用プログラムであるPeakSeq -preprocessで書式変換しておきます。
  • case用のファイルと、control用のファイルの両方を予め処理しておく必要があります。
  • infile=`ls -LS1 ./input_1/*.* | egrep "(.bam)$" | head -n 1`; 
    outname=`basename $infile | sed "s/\(.bam\)//"`;
    chromInf=@CIP_DATA@/genomes/#OPTION_3#/#OPTION_1#/chr_id_list.txt;
    if [ "$infile" == "" ] ; then
    touch ./output_1/.empty;
    else
    mkdir ./output_1/${outname};
    @CIP_HOME@/tool/#OPTION_2#/samtools view $infile | @CIP_HOME@/tool/#OPTION_1#/PeakSeq -preprocess SAM $chromInf stdin ./output_1/${outname};
    fi;

ピークコール

  • indir=`ls -LS1 ./input_1/*/* | head -n 1`;
    indir=`dirname $indir`;
    controlDir=`ls -LS1 ./input_2/*/* | head -n 1`;
    controlDir=`dirname $controlDir`;
    outname=`basename $indir`;
    if [ "#OPTION_1#" != "" -a "#OPTION_1#" != " " ] ; then 
    outname="#OPTION_1#"; 
    fi;
    outname=`echo $outname | sed "s/ /_/g"`;
    cd ./log;
    conf=./config.dat;
    echo "Experiment_id "$outname >> $conf;
    ln -s @CIP_DATA@/genomes/#OPTION_3#/#OPTION_2#/chr_id_list.txt ./;
    echo "chromosome_list_file chr_id_list.txt" >> $conf;
    ln -s @CIP_DATA@/genomes/#OPTION_3#/#OPTION_2#/Mapability_*.txt ./;
    mapability=`ls -1 ./Mapability_*.txt | head -n 1`;
    echo "Mappability_map_file $mapability" >> $conf;
    echo "ChIP_Seq_reads_data_dirs ../"$indir >> $conf;
    echo "Input_reads_data_dirs ../"$controlDir >> $conf;
    
    echo "Enrichment_fragment_length #OPTION_4# " >> $conf;
    echo "target_FDR #OPTION_5#" >> $conf;
    echo "N_Simulations #OPTION_6# " >> $conf;
    echo "Minimum_interpeak_distance #OPTION_7#" >> $conf;
    echo "max_Qvalue #OPTION_8#" >> $conf;
    echo "Simulation_seed #OPTION_9#" >> $conf;
    echo "Background_model #OPTION_10# " >> $conf;
    
    @CIP_HOME@/tool/#OPTION_2#/PeakSeq -peak_select ./config.dat;
    
    cd ../;
    
    sed -n '2,$p' ./log/${outname}.txt | awk '{print $1"\t"$2"\t"$3"\tPeakSeq_peak_"NR"\t"$4}' > ./output_1/${outname}.bed;
    @CIP_HOME@/bin/mvIfExists.sh ./log/${outname}_narrowPeak.txt ./output_2/;
    @CIP_HOME@/bin/mvIfExists.sh ./log/${outname}.txt ./output_3/;
    • スクリプトが長いですが、実行オプションをコンフィグファイルにしたり、後処理でPeakSeq独自書式からbedファイルを生成するなどを行っているのみです。

ピークコール by SISSRs

前処理

  • SISSRsはbamファイルを入力書式として受け取らないので、BEDtoolsのbamToBedで予めbed形式に変換しておきます。
  • case用のファイルと、control用のファイルの両方を予め処理しておく必要があります。
  • infile=`ls -LS1 ./input_1/*.* | egrep "(.bam)$" | head -n 1`; 
    outname=`basename $infile | sed "s/\(.bam\)//"`; 
    if [ "$infile" == "" ] ; then
    touch ./output_1/.empty;
    else
    @CIP_HOME@/tool/#OPTION_1#/bin/bamToBed -i $infile > ./output_1/${outname}.bed;
    fi;

ピークコール

  • infile=`ls -LS1 ./input_1/*.* | egrep "(.bed)$" | head -n 1`; 
    control=`ls -LS1 ./input_2/*.* | egrep "(.bed)$" | head -n 1`;
    if [ "$control" != "" ] ; then
    control="-b "$control;
    fi;
    outname=`basename $infile | sed "s/\(.bed\)//"`;
    if [ "#OPTION_4#" != "" -a "#OPTION_4#" != " " ] ; then 
    outname="#OPTION_4#";
    fi;
    genomeSize=`cat @CIP_DATA@/genomes/#OPTION_1#/chromInfo/totalGenomeSize`;
    if [ "#OPTION_3#" != "0" ] ; then 
    genomeSize="#OPTION_3#"; 
    fi;
    @CIP_HOME@/tool/#OPTION_2#/sissrs.pl -i $infile -o ./output_1/${outname}.sissrsout -s $genomeSize $control #OPTION_5#;

後処理

  • SISSRsのピークコールの結果は独自書式なので、bed形式に変換します。
  • infile=`ls -LS1 ./input_1/*.* | egrep "(.sissrsout)$" | head -n 1`;
    outname=`basename $infile | sed "s/\(.bam\)//"`; 
    sed -n "/^---/,/^===/p" $infile | grep -v "^---" | grep -v "^===" | awk '{print $1"\t"$2"\t"$3"\tSISSRs_peak_"NR"\t"$4}' > ./output_1/${outname}.bed;
    sed -n "1,/^---/p" $infile > ./output_2/statistics.txt;

ChIPpeakAnnoによる遺伝子アノテーション付加

前処理

  • bedファイルといっても、書式にはかなりのバリエーションがあり、そのままでは、ChIPpeakAnnoを使用するのに不適当な場合があります。
  • こちらの処理で不適当な箇所を修正します。
  • infile=`ls -LS1 ./input_1/*.* | egrep "(.bed)$" | head -n 1`; 
    outname=`basename $infile | sed "s/\(.bed\)//"`;
    @CIP_HOME@/bin/modifyBedForChIPpeakAnno.pl -i ${infile} -o output_1/${outname}_forChIPpeakAnno.bed;

アノテーション付加

  • bedfile=`ls -LS1 ./input_1/*.* | egrep "(.bed)$" | head -n 1`;
    outname=`basename $bedfile | sed "s/\(.bed\)//"`; 
    distancefile=${outname}_TSSdistance.tsv;
    gofile=${outname}_GO.tsv;
    fastafile=${outname}_sequence.fasta;
    @CIP_HOME@/tool/R-2.14.0/bin/R --vanilla --slave --args $bedfile #OPTION_3# $distancefile $gofile $fastafile #OPTION_1# #OPTION_2# < @CIP_HOME@/bin/runChIPpeakAnno.r 2&> ./log/chippeakanno_log;
    resultnum=`ls -LS1 | grep $distancefile |wc -l`;
    if [ "$resultnum" -ge "1" ] ;then mv $distancefile ./output_1/.;
    fi;
    resultnum=`ls -LS1 | grep $gofile | wc -l`;
    if [ "$resultnum" -ge "1" ] ;then mv $gofile ./output_2/.;
    fi;
    resultnum=`ls -LS1 | grep $fastafile | wc -l`;
    if [ "$resultnum" -ge "1" ] ;then mv $fastafile ./output_3/.;
    fi;


Last edited by nobfujii .
Page last modified on 火曜日 30 of June, 2015 04:52:07 JST.

List of attached files
ID 名前 desc uploaded サイズ ダウンロード数 操作
42 default symlinkIfNotExist.sh 2015-07-02 06:17 by nobfujii 898 b 314 表示 Download  
41 default mergeWithMultiId.pl 2015-07-02 06:17 by nobfujii 6.02 Kb 601 表示 Download  
40 default convertChIPpeakAnnoForIntegration.pl 2015-07-02 06:16 by nobfujii 4.50 Kb 900 表示 Download