Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

DROMPA3: その11 複製解析(出芽酵母)

久しぶりのDROMPAのエントリです。今回は出芽酵母(S. cerevisiae)の複製解析を行います。複製開始後のDNAを複製前のDNAで割り算することで、ゲノム上のどこでどの程度複製フォークが進んでいるかを可視化することができます。 データは以下の論文のものを用います。

Origin association of Sld3, Sld7, and Cdc45 proteins is a key step for determination of origin-firing timing. - PubMed - NCBI

ここで用いるサンプルは、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

手元の環境ではダウンロード完了まで4時間ほどかかりました。 私の場合は通常、ファイルサイズを削減するためにダウンロードしてから pigz で圧縮するのですが、今回は実行時間節約のため圧縮せずにいきます。(pigzsudo apt install pigz でインストール可能です。)

ファイルを確認しましょう。

$ ls -lh SRR*
-rw-r--r-- 1 rnakato rnakato 7.3G 116 11:40 SRR398609.fastq
-rw-r--r-- 1 rnakato rnakato 6.7G 116 11:49 SRR398610.fastq
-rw-r--r-- 1 rnakato rnakato 7.4G 116 12:06 SRR398611.fastq
-rw-r--r-- 1 rnakato rnakato 5.5G 116 12:17 SRR398612.fastq
-rw-r--r-- 1 rnakato rnakato 8.2G 116 12:32 SRR398613.fastq
-rw-r--r-- 1 rnakato rnakato 6.0G 116 12:42 SRR398614.fastq
-rw-r--r-- 1 rnakato rnakato 8.2G 116 12:55 SRR398615.fastq
-rw-r--r-- 1 rnakato rnakato 6.8G 116 13:07 SRR398616.fastq
-rw-r--r-- 1 rnakato rnakato 6.3G 116 13:19 SRR398617.fastq
-rw-r--r-- 1 rnakato rnakato 7.2G 116 13:32 SRR398618.fastq
-rw-r--r-- 1 rnakato rnakato 6.8G 116 13:43 SRR398619.fastq
-rw-r--r-- 1 rnakato rnakato 6.0G 116 13:54 SRR398620.fastq
-rw-r--r-- 1 rnakato rnakato 7.0G 116 14:19 SRR398621.fastq
-rw-r--r-- 1 rnakato rnakato 6.4G 116 14:35 SRR398622.fastq
-rw-r--r-- 1 rnakato rnakato 7.1G 116 14:54 SRR398623.fastq
-rw-r--r-- 1 rnakato rnakato 7.1G 116 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行で1つのリードを表します。 @ で始まる行がリード名、次の T で始まる行が リード配列、4行目は読まれた各塩基のクオリティ(信頼度)を表します。 (このデータでは、3行目は意味を持ちません。)

リード配列が 0,1,2 の数字で表されていますね。 これが color space データ で、【2つの塩基の組み合わせ】を数字で表現したものです。 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 という名前のディレクトリを作成し、マップファイルをその中に出力するようにしています。

\ は1行のコマンドを複数行に分割して書く時に使う記号です。

-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軸の罫線の数を3に
-ystep 14 で罫線1あたりの長さ(pixel)を14に(デフォルトは20)
-scale_ratio 1 で罫線1あたりのスケールを1に
それぞれ変更しています。結果は以下のようになります。

evince pdf/yeast1.pdf

f:id:rnakato:20191102181354j:plain
Repliseq1

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 オプションは省略しました。 結果は以下のようになります。

f:id:rnakato:20191102181527j:plain
Repliseq2

上のコマンドはちょっと長いので、変数を使って短くしてみましょう。得られる結果は 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

f:id:rnakato:20191102182154j:plain
Repliseq3

このようにすると、複製されている領域のサンプル間の差がより顕著になりますね。

ピーク抽出

サンプル間で複製状態を比較するために、 yeast4において赤でハイライトされた領域(ピーク領域)のリストが欲しい、という場合もあるでしょう。 リスト形式でピークを出力するには drompa_peakcall を使います。

#(変数の行は省略)
drompa_peakcall PC_ENRICH $s1 -p YST1019_Gal -gt $gt -if 3

こうすると、 YST1019_Gal.bed と YST1019_Gal.xls の2つのファイルが生成されます。 YST1019_Gal.bed はBED形式のファイル、 YST1019_Gal.xls はより詳細な情報が書かれたファイルになります。

まとめ

今回は複製解析について紹介しました。 酵母の場合はChIP-seqでも同様に PC_ENRICHで可視化します。これはヒトやマウスに比べるとウィンドウあたりのリード数が十分多く、割り算してもnoisyにならないという理由によります。 PC_ENRICHはデフォルトで 統計検定はせず、シンプルにenrichmentとそのリード数だけを閾値に用います。有意差検定を用いたピーク抽出をする場合は通常の PC_SHARP コマンドを使ってください。