DROMPA3: その11 複製解析（出芽酵母）
久しぶりのDROMPAのエントリです。今回は出芽酵母（S. cerevisiae）の複製解析を行います。複製開始後のDNAを複製前のDNAで割り算することで、ゲノム上のどこでどの程度複製フォークが進んでいるかを可視化することができます。 データは以下の論文のものを用います。
ここで用いるサンプルは、SOLiDで生成されたcolor spaceデータになります。 color spaceデータはbowtieでマッピング可能です（専用のindexファイルが必要）。
インストール
DROMPA3のインストール方法はこの記事を参照してください。
データダウンロード
今回は、必要なreference data類をGoogle Driveにアップしてみましたので、以下のコマンドでダウンロードしてください。
wget -O genome_table "https://drive.google.com/uc?export=download&id=11xhQpJt-kYJKc38HWXUfwoim2utSby8Q" wget -O ARS-oriDB.txt "https://drive.google.com/uc?export=download&id=12A8EigD3J_DB_BaMyDRAd51LlOBg9F9u" wget -O SGD_features.tab "https://drive.google.com/uc?export=download&id=13OPXenZyxXmF8hkQ1my3GlW6QrH7Y1Xr"
bowtieのcsインデックスはbowtie-buildで作成可能ですが、ここでは作成済のものをダウンロードします。
wget "https://drive.google.com/uc?export=download&id=138QOuZUjOMyXeFVs6Uv6SbpyKixt9N-G" -O S_cerevisiae_cs.zip unzip S_cerevisiae_cs.zip
csfastqファイルのダウンロード
以下ではシェル変数とfor文を用いていきます。よくわからないという方は DROMPA3: その5 シェル変数を使う - Palmsonntagmorgen を参考にしてください。
今回用いるデータをSRAからダウンロードします。
以下では fasterq-dump を使っていますが、これは fastq-dump の後継で、マルチコアで高速にfastqファイルをダウンロードするものです。fastq-dump を使ってダウンロードしてもかまいません。
for num in $(seq 398609 398624) do fasterq-dump SRR$num done
手元の環境ではダウンロード完了まで４時間ほどかかりました。
私の場合は通常、ファイルサイズを削減するためにダウンロードしてから pigz で圧縮するのですが、今回は実行時間節約のため圧縮せずにいきます。（pigz は sudo apt install pigz でインストール可能です。）
ファイルを確認しましょう。
$ ls -lh SRR* -rw-r--r-- 1 rnakato rnakato 7.3G 11月 6 11:40 SRR398609.fastq -rw-r--r-- 1 rnakato rnakato 6.7G 11月 6 11:49 SRR398610.fastq -rw-r--r-- 1 rnakato rnakato 7.4G 11月 6 12:06 SRR398611.fastq -rw-r--r-- 1 rnakato rnakato 5.5G 11月 6 12:17 SRR398612.fastq -rw-r--r-- 1 rnakato rnakato 8.2G 11月 6 12:32 SRR398613.fastq -rw-r--r-- 1 rnakato rnakato 6.0G 11月 6 12:42 SRR398614.fastq -rw-r--r-- 1 rnakato rnakato 8.2G 11月 6 12:55 SRR398615.fastq -rw-r--r-- 1 rnakato rnakato 6.8G 11月 6 13:07 SRR398616.fastq -rw-r--r-- 1 rnakato rnakato 6.3G 11月 6 13:19 SRR398617.fastq -rw-r--r-- 1 rnakato rnakato 7.2G 11月 6 13:32 SRR398618.fastq -rw-r--r-- 1 rnakato rnakato 6.8G 11月 6 13:43 SRR398619.fastq -rw-r--r-- 1 rnakato rnakato 6.0G 11月 6 13:54 SRR398620.fastq -rw-r--r-- 1 rnakato rnakato 7.0G 11月 6 14:19 SRR398621.fastq -rw-r--r-- 1 rnakato rnakato 6.4G 11月 6 14:35 SRR398622.fastq -rw-r--r-- 1 rnakato rnakato 7.1G 11月 6 14:54 SRR398623.fastq -rw-r--r-- 1 rnakato rnakato 7.1G 11月 6 15:14 SRR398624.fastq
複製前・複製後のペアが8ペア、計16サンプルです。
fastqファイルの中身を見てみましょう。 head コマンドで最初の数行を表示します。
$ head SRR398609.fastq @SRR398609.1 2010_018_bcSample1_1_19_677 length=50 T0002.2000..000..0000.00000000..0..02...0.........0 +SRR398609.1 2010_018_bcSample1_1_19_677 length=50 !''*+!%'*'!!)&)!!*(-)!.-)'/*%,!!.!!)&!!!)!!!!!!!!!8 @SRR398609.2 2010_018_bcSample1_1_20_852 length=50 T0100.1010..000..0200.00000020..0..00...0.........1 +SRR398609.2 2010_018_bcSample1_1_20_852 length=50 !%)).!%(1%!!(&'!!4%))!))*-)%%)!!1!!**!!!*!!!!!!!!!& @SRR398609.3 2010_018_bcSample1_1_20_996 length=50 T3102.0001..102..2220.20001230..2..23...3.........2
fastqファイルは 4行で１つのリードを表します。
@ で始まる行がリード名、次の T で始まる行が リード配列、4行目は読まれた各塩基のクオリティ（信頼度）を表します。
（このデータでは、3行目は意味を持ちません。）
リード配列が 0,1,2 の数字で表されていますね。 これが color space データ で、【２つの塩基の組み合わせ】を数字で表現したものです。 SOLiD sequencer で読まれた配列はこのような形式で出力されます。 color space形式のfastqファイルはcsfastq形式と呼ばれます。
マッピング
以下のfor文は、各csfastqファイルに対して、bowtieでマッピングした後そのままsamtoolsで sorted bamに変換しています。詳細は SAMtoolsとリダイレクト - Palmsonntagmorgen を参照してください。
mkdir -p bam for num in $(seq 398609 398624) do bowtie -S -C S_cerevisiae-cs/S_cerevisiae-cs \ SRR$num.fastq -n2 -k1 --chunkmbs 2048 -p12 \ | samtools sort \ > bam/SRR$num-n2-k1.sort.bam done
最初にbam という名前のディレクトリを作成し、マップファイルをその中に出力するようにしています。
\ は１行のコマンドを複数行に分割して書く時に使う記号です。
-S は出力ファイルをSAM形式にするオプション、 -C はcsfastqデータである時に指定します。
S_cerevisiae-cs/S_cerevisiae-cs はさきほどダウンロード・解凍したcolor space用のbowtieインデックスファイルです。
-n2 はリード配列のミスマッチを2まで許容し、 -k1 はリードがゲノムの複数箇所にマッチした時にベストマッチのみ出力するオプションです。
--chunkmbs 2048 はbowtieに許容するメモリ上限、 -p 12 はCPUを12コア使うことを意味します。
bedGraphファイルの変換
bamファイルをbedGraphに変換します。genome_table は先ほどダウンロードしたものを用います。
-of 3 が出力をbedGraphにするオプションです。
for num in $(seq 398609 398624) do parse2wig -f BAM -i bam/SRR$num-n2-k1.sort.bam \ -o SRR$num-n2-k1 -n GR -of 3 -binsize 100 -gt genome_table done
このコマンドの結果、parse2wigdir の中にbedGraphファイルが生成されます。
pdf可視化（1サンプル）
それでは可視化してみましょう。 今回は「複製後/複製前」の割り算 (ratio) を可視化したいので、 PC_ENRICH コマンドを使います。 サンプル数が多いので、まず1サンプルペアのみ可視化します。 アノテーションファイルは先ほどダウンロードしたものを指定してください。
dir=parse2wigdir drompa_draw PC_ENRICH \ -gene SGD_features.tab -gftype 3 \ -ars ARS-oriDB.txt \ -gt genome_table \ -i $dir/SRR398612-n2-k1-scer-raw-mpbl-GR,$dir/SRR398611-n2-k1-scer-raw-mpbl-GR,YST1019_Gal \ -p pdf/yeast1 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \ -if 3
-i で入力サンプルを指定します。<ChIP>,<Control>,<Label> の順にコンマで区切って指定します。 今回の場合だと、 <複製後>,<複製前>,<株名> となります。 入力サンプルには染色体までのprefixを指定します。
-p 出力ファイル名のprefix
-sm 1000 スムージング幅 (ここでは1kbp)
-ls 200 1段あたりの length
-rmchr 染色体別のファイルを作成しない
-nosig ピーク抽出を行わない
-lpp 2 で1ページあたりの段数を2に
-bn 3 でy軸の罫線の数を３に
-ystep 14 で罫線１あたりの長さ（pixel）を14に（デフォルトは20）
-scale_ratio 1 で罫線１あたりのスケールを1に
それぞれ変更しています。結果は以下のようになります。
evince pdf/yeast1.pdf
Enrichmentが2倍付近まで到達している領域が、既に複製されたゲノム領域を表します。複製開始点を中心に、両側に複製フォークが広がっていく様子がわかります。
pdfファイルの作成（複数サンプル）
では、全サンプルを同時に可視化してみましょう。
dir=parse2wigdir drompa_draw PC_ENRICH \ -ars ARS-oriDB.txt \ -gt genome_table \ -i $dir/SRR398612-n2-k1-scer-raw-mpbl-GR,$dir/SRR398611-n2-k1-scer-raw-mpbl-GR,YST1019_Gal \ -i $dir/SRR398610-n2-k1-scer-raw-mpbl-GR,$dir/SRR398609-n2-k1-scer-raw-mpbl-GR,YST1019_Raf \ -i $dir/SRR398616-n2-k1-scer-raw-mpbl-GR,$dir/SRR398615-n2-k1-scer-raw-mpbl-GR,YST1053_Gal \ -i $dir/SRR398614-n2-k1-scer-raw-mpbl-GR,$dir/SRR398613-n2-k1-scer-raw-mpbl-GR,YST1053_Raf \ -i $dir/SRR398620-n2-k1-scer-raw-mpbl-GR,$dir/SRR398618-n2-k1-scer-raw-mpbl-GR,YST1076_Gal \ -i $dir/SRR398619-n2-k1-scer-raw-mpbl-GR,$dir/SRR398617-n2-k1-scer-raw-mpbl-GR,YST1076_Raf \ -i $dir/SRR398624-n2-k1-scer-raw-mpbl-GR,$dir/SRR398623-n2-k1-scer-raw-mpbl-GR,YST1287_Gal \ -i $dir/SRR398622-n2-k1-scer-raw-mpbl-GR,$dir/SRR398621-n2-k1-scer-raw-mpbl-GR,YST1287_Raf \ -p pdf/yeast2 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \ -if 3
各ペアを個別に -i で指定しています。また、今回の複製解析では遺伝子情報は必要ないので、 -gene オプションは省略しました。 結果は以下のようになります。
上のコマンドはちょっと長いので、変数を使って短くしてみましょう。得られる結果は yeast3 と同じです。
# サンプルの定義 dir=parse2wigdir post=-n2-k1-scer-raw-mpbl-GR # サンプルに共通するpostfixを格納 IP1_60="$dir/SRR398612$post" # YST1019 Gal 60min IP1_0="$dir/SRR398611$post" # YST1019 Gal 0min IP2_60="$dir/SRR398610$post" # YST1019 Raf 60min IP2_0="$dir/SRR398609$post" # YST1019 Raf 0min IP3_60="$dir/SRR398616$post" # YST1053 Gal 60min IP3_0="$dir/SRR398615$post" # YST1053 Gal 0min IP4_60="$dir/SRR398614$post" # YST1053 Raf 60min IP4_0="$dir/SRR398613$post" # YST1053 Raf 0min IP5_60="$dir/SRR398620$post" # YST1076 Gal 60min IP5_0="$dir/SRR398618$post" # YST1076 Gal 0min IP6_60="$dir/SRR398619$post" # YST1076 Raf 60min IP6_0="$dir/SRR398617$post" # YST1076 Raf 0min IP7_60="$dir/SRR398624$post" # YST1287 Gal 60min IP7_0="$dir/SRR398623$post" # YST1287 Gal 0min IP8_60="$dir/SRR398622$post" # YST1287 Raf 60min IP8_0="$dir/SRR398621$post" # YST1287 Raf 0min # サンプルペアの定義 s1="-i $IP1_60,$IP1_0,YST1019_Gal" s2="-i $IP2_60,$IP2_0,YST1019_Raf" s3="-i $IP3_60,$IP3_0,YST1053_Gal" s4="-i $IP4_60,$IP4_0,YST1053_Raf" s5="-i $IP5_60,$IP5_0,YST1076_Gal" s6="-i $IP6_60,$IP6_0,YST1076_Raf" s7="-i $IP7_60,$IP7_0,YST1287_Gal" s8="-i $IP8_60,$IP8_0,YST1287_Raf" # DROMPAの実行 drompa_draw PC_ENRICH \ -ars ARS-oriDB.txt \ -gt genome_table \ $s1 $s2 $s3 $s4 $s5 $s6 $s7 $s8 \ -p pdf/yeast3 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \ -if 3
サンプルを定義するのにだいぶ行数を使ってしまいましたが、このようにファイル名をわかりやすい変数名に変換することで、 誤字やサンプルの取り違えなどのミスをある程度予防することができます。
-nosig を消すと、Enrichment >= 2.0 の部分が赤でハイライトされます。 （ハイライトの閾値は -ethre オプションで指定可能です）
drompa_draw PC_ENRICH -gftype 3 -ars $ARS -ter $TER -gt $gt \ $s1 $s2 $s3 $s4 $s5 $s6 $s7 $s8 \ -p pdf/yeast4 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -sm 1000 \ -if 3
このようにすると、複製されている領域のサンプル間の差がより顕著になりますね。
ピーク抽出
サンプル間で複製状態を比較するために、 yeast4において赤でハイライトされた領域（ピーク領域）のリストが欲しい、という場合もあるでしょう。 リスト形式でピークを出力するには drompa_peakcall を使います。
#（変数の行は省略） drompa_peakcall PC_ENRICH $s1 -p YST1019_Gal -gt $gt -if 3
こうすると、 YST1019_Gal.bed と YST1019_Gal.xls の２つのファイルが生成されます。 YST1019_Gal.bed はBED形式のファイル、 YST1019_Gal.xls はより詳細な情報が書かれたファイルになります。
まとめ
今回は複製解析について紹介しました。 酵母の場合はChIP-seqでも同様に PC_ENRICHで可視化します。これはヒトやマウスに比べるとウィンドウあたりのリード数が十分多く、割り算してもnoisyにならないという理由によります。 PC_ENRICHはデフォルトで 統計検定はせず、シンプルにenrichmentとそのリード数だけを閾値に用います。有意差検定を用いたピーク抽出をする場合は通常の PC_SHARP コマンドを使ってください。
