STAP細胞関係のゲノムデータを解析してみた

本記事の目的と注意

注意! 私は、NGS については amplicon sequencing の解析経験(しかも半年)しかない。本記事は、データを解析して、STAP論文(Obokata et al, Nature 2014. ArticleLetter)に対して何らかの結論を導くのが目的ではない。これだけリード数が少なくて、しかもサンプルがポリクローナルな混合物であることを考えると、ここから何かを結論するのは極めて慎重にならないといけないと思う。したがって、結果の「解釈」には立ち入らない(し、その能力もない)。本記事は、「ネットで話題になっているデータを、自分も解析してみたい!」と、「行為」そのものに魅力を感じる方のために、私が行った操作の流れを紹介するものである。

私は当初 RNA-seq のデータを解析しようとしたが、リファレンス・トランスクリプトームに存在しない再構成後の T cell receptor 領域の適切な解析法が分からず頓挫していた。そこへ3月5日頃から slashdot の kaho 氏が ChIP-seq の input を利用した解析(その1, その2, その3)を発表した。「このアプローチなら自分の知識でもなんとかなりそう」と思い、早速やってみたのがこれである。解析内容(着眼点)は kaho 氏のアイデアである。具体的な解析方法(ソフトの組み合わせ方やコード)は氏が公開していないので、自分で考えた。たぶんもっと適切なやり方があるはずだ。自信がないところは分からないとはっきり書いておいたので、アドバイスいただけると感謝する。

重要なキーワード、ソフトウェア名などはリンク付きで紹介しておいたので、これをきっかけに次世代シーケンサ界隈に興味を持ったり、勉強を始める人が増えればよいと思う。全体的な入門としては、緒方法親氏の「お家でできるMac Bookでやる次世代シーケンスデータ解析」が大変分かりやすいチュートリアルである。チャプタ1 (28ページ) までの知識で今回の解析は十分可能だ。

基礎知識としては、LinuxMac OS のターミナルが使える方を想定している。ハードウェア的には HDD が 100GB、メモリが 8GB あり、一晩計算させておく時間の余裕があれば十分である。自作PCファンは、円周率ばかり計算させていないで、たまにはこういう事にCPUを使ってみるのもいかが?

SRAから配列データを取ってくる

NCBI Bioproject の該当エントリから入る。"Low pH treated CD45 positive Cells; ChIPSeq_input" など、 ChIP-input と書いてあるものが全ゲノムのDNA配列に対応するのでこれを集める。リンクを押すと、size 1.1G などと書いてあるリンクがあるので、それを押すとデータファイルがダウンロードできるフォルダ(SRR1171584)が開く。 他に、STAP-SC, CD45 positive Cells, Embryonic Stem Cells をもらってきた。各データは1-2GBくらい。以下では、SRR1171584 での処理を説明していく。他の3つも全く同じだ。

これらのファイルはSRA という独自形式で圧縮されているので、FASTQ 形式に戻す必要がある。FASTQ 形式は、FASTA 形式にクオリティ・スコアなどの情報をたせるようにしたフォーマットである。

NCBI SRA Toolkit をもらってきて bin にパスを通す。

fastq-dump SRR1171584.sra 

のようにすると、SRR1171584.fastq ができる。ファイルサイズが約5-10倍になる。展開処理が重いようで、10分くらいかかる。

このセクションの内容をさらに詳しく知りたい方には、g86 has gone.を勧める。

リファレンス・ゲノムの準備

上のデータは、ショートリードの配列そのものなので、まず、リファレンスゲノムにアライメント(マッピング・貼り付け)しなければならない。

マウスのリファレンスゲノムを UCSC から入手する。http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/ に入って、"chromFa.tar.gz" をもらってくる。なお、kaho 氏が使っているのは mm9 というバージョンであり、本稿で使っているのは mm10 という1つ新しいバージョンなので座標が若干ズレていることに注意してほしい。

解凍すると、染色体 chromosome ごとに別々のファイルになっている。*_random.fa のようなファイルが何なのかは不明。とりあえず chr??.fa のみ集めて

cat *.fa > chr_all.fa 

として1つのファイルにマージする(CHECK: 真面目に解析するときは、こういういい加減な態度ではもちろんダメですよ)。
以降のステップで時間やメモリを使いすぎるのが嫌なら、話題の TCR が入っている chr6 と chr14 だけで進めるのもよいかもしれない。

インデックスを作る

http://bio-bwa.sourceforge.net/ から bwa 0.7.7 をもらってきてコンパイルする(CHECK: 今回の用途では Bowtie のほうが良い?)。また、データファイルの変換などに用いる samtools も最新版 0.1.19 をもらってきてインストールする。どちらも make するだけだ。

出来上がった bwa と samtools にパスを通して、

bwa index -a bwtsw chr_all.fa

とすると、入力ファイルと同じ場所にインデックスが作られる。3時間くらいかかる。ここが一番メモリを使うステップ。たぶん 8GB あれば大丈夫。

ここでは、情報系クラスタで一時期流行になった Burrows-Wheeler Transform とかをして、FM-index を構築しているのである。関連する資料としては、Burrows Wheeler TransformとLF mapping | Preferred Researchウェーブレット木の世界 などが分かりやすいかな。

アライメント

いよいよ実際のアライメントを行う。上で作ったリファレンスと、リードの fastq ファイルを指定して、bwa aln コマンドを実行する。結果は標準出力に吐き出されるので、リダイレクトして保存する。

bwa aln -t 4 ../ref/bwa_index/chr_all.fa ../data/SRR1171584.fastq > SRR1171584.sai

オプションの -t 4 は並列数を指定している。4並列で90分弱かかった。続けて、

bwa samse ../ref/bwa_index/chr_all.fa SRR1171584.sai ../data/SRR1171584.fastq | samtools view -bS -h - | samtools sort - SRR1171584
samtools index SRR1171583.bam

のようにして完成。このステップも30分くらいかかる。ディスク IO に時間がかかるのでパイプでつないでいるが、メモリがやばいときは中間ファイルを作ったほうがいいかも。

つまり、

bwa samse ../ref/bwa_index/chr_all.fa SRR1171584.sai ../data/SRR1171584.fastq > SRR1171584.sam

で、マッピング後の配列を sam ファイルに保存する。sam ファイルは、アライメントを表すテキスト形式のファイルだ。テキスト形式は人間が hack するには都合がよいが、パフォーマンス上は不利なのでバイナリの .bam 形式に変換する。

samtools view -bS -h SRR1171584.sam > SRR1171584-tmp.bam 

この段階ではリードが処理されたままの順序になっているので、染色体上の座標で sort する(最後の引数に拡張子はいらない)

samtools sort SRR1171584-tmp.bam SRR1171584

これで、bam ファイルができた。圧縮が効くので 2-3GB ほどである。最後の samtools index は、このファイルに索引を作って、「染色体○の○○塩基目」といった座標データから、目的の配列をすぐに得られるようにしている。

ここまで来たら、.fastq, .sai, .sam ファイルは消してよい。どれも10GB くらいあるのでかなりムダである。

これで、IGV(無償だが、メールアドレスの登録が必要)などのゲノムビューアで閲覧できるようになった。ほとんどのゲノムビューアは、sort して index がついた BAM ファイルでないと開けないので注意しよう。 bioinformatics には by-eye informatics *1も大事。

T cell receptor 領域の解析例

ゲノム全体の .bam ファイルは大きすぎて取り回しが悪いので、まず関心領域だけ切り出すことにする。
関心領域は、UCSC Genome Browser で、genome として "mouse" を選択したあと、 Tcra と Tcrb (T cell receptor alpha と beta)を検索して、適当にズームアウトして(zoom out x10 などを押す)、周辺 3.2Mbp くらいの範囲を選んだ。ここでは、

chr14:51666755-54899329 TCRa 周辺
chr6:39526641-42799973 TCRb 周辺

を使用する。この範囲を samtools view コマンドに引数として渡して切り出しを行い、その範囲をソート・インデックスしなおす(これをしないと、次のステップでエラーになる)。

samtools view -b SRR1171584.bam chr6:39526641-42799973 chr14:51666755-54899329 | samtools sort - lowph-in
samtools index lowph-in.bam

次に、この領域の各残基がそれぞれ何回読まれているかを調べる。それには samtools mpileup が便利である。-f オプションでリファレンスを指定しているが、省略しても可。

samtools mpileup -f ../ref/bwa_index/chr_all.fa lowph-in.bam > lowph-in.mpileup

この出力はタブ区切りのテキストファイルであり、染色体・座標・リファレンスの塩基(省略した場合はNになる)・リード数……が場所ごとに並んでいる。

この出力には chr6 と chr14 が混ざっているので grep で区別したあと、座標と塩基と depth を cut で切り出す(塩基は本来必要ないが、確認用)。

grep chr6 lowph-in.mpileup | cut -f 2,3,4 > lowph-chr6.lst
grep chr14 lowph-in.mpileup | cut -f 2,3,4 > lowph-chr14.lst

あとは R や Python + matplotlib などでプロットして遊ぶ。その際、mpileup はリードがない部分は出力しないので、移動平均などの連続区間に対する操作には注意が必要である。

以下では2つのゲノム間でリード数の差を取って、その平均が0になるように調節したあと移動平均をプロットしているが、おそらくこの処理は不適切である(処理範囲が狭すぎると、過度の正規化を起こしてしまうから)。実際、kaho 氏がアップロードしている結果とも完全には一致しない。ゲノム全体のリード数(あるいはマップできたリード数?)に対する割合で正規化すべきだろうか? また、summary(depth1) の出力を見れば分かるように、depth はかなり低く、リード数 0 の領域もそれなりにある。つまり、2つのゲノム間でリード数に差があったとしても、それが統計的に有意かを判断するのは難しそうである。このようなデータから、ゲノムの deletion の有無を判断できるのか、私は知らない。本来は類似の研究例を調べるべきだが、そこまでの余裕はない。お決まりの方法があるなら、教えてほしい。

setwd('~/prog/stap')

# install.packages('TTR')
library(TTR)

data1 <- read.table('esc-chr14.lst')
data2 <- read.table('lowph-chr14.lst')

# mpileup はリードがない部分は省略されてしまうので、補う必要がある
chr14_start <- 51666700
chr14_end <- 54900000
pos <- seq(from=chr14_start, to=chr14_end)
depth1 <- rep(0, chr14_end - chr14_start + 1)
depth2 <- rep(0, chr14_end - chr14_start + 1)
depth1[data1$V1 - chr14_start + 1] <- data1$V3
depth2[data2$V1 - chr14_start + 1] <- data2$V3

summary(depth1)

# 注意: おそらくこの処理は不適切
diff <- depth2 - depth1
diff <- diff - mean(diff)

window <- 2500
plot(pos, SMA(diff, window), type='l', xlim=c(53500000, 55000000))
abline(h = 0)

Y染色体領域の解析例

細胞がオス由来かメス由来か調べるために、Y染色体にマップしたリードの本数を数えてみよう。

上でも使った samtools view コマンドは -b オプションを付けずに実行すると、テキストファイルである SAM 形式で切り出しを行える。この形式では、リード1本が1行になるので、wc -l で行数を数えればリードの本数はすぐ分かる。

samtools view SRR1171584.bam chrY | wc -l

あとは、全リード数に対する割合を求めるとか、それが細胞ごとに違うのか(割合に対するカイ二乗検定)調べてみよう。

Oct4 遺伝子周辺の解析例 (3/6 夜(英国時間)追記)

kaho 氏が解析報告の第3段を公開した。細胞が Oct4-GFP transgenic mouse に由来することを確認するのが目的のようだ。確実な判断をするには transgene がゲノムに組み込まれている箇所を探すのが一番良いのだろうが、繋ぎ目の部分はリファレンス・ゲノムにないので、ここまでの解析の流れのままでは難しそうだ。kaho 氏は、transgene の存在によって Oct4 のコピー数が増えることを調べている。これをやってみよう。

といっても、流れは TCR 領域の解析と同じ。Oct4 遺伝子周辺を切り出してきて、mpileup で塩基ごとの depth を抽出、R でプロットするだけ。

samtools view -b SRR1171584.bam chr17:35449000-35845000 | samtools sort - lowph-oct4
samtools index lowph-oct4.bam
samtools mpileup lowph-oct4.bam | cut -f2,4 > lowph-oct4.lst

という感じで、興味ある細胞ごとに処理したあと、R に読み込んでプロット。

esc <- read.table('esc-oct4.lst')
stapsc <- read.table('stapsc-oct4.lst')
lowph <- read.table('lowph-oct4.lst')
cd45 <- read.table('cd45-oct4.lst')
summary(esc$V1)

# mpileup の足りない場所を埋める
region_start <- min(c(esc$V1, cd45$V1, lowph$V1, stapsc$V1))
region_end <- max(c(esc$V1, cd45$V1, lowph$V1, stapsc$V1))
pos <- seq(from=region_start, to=region_end)
depth <- matrix(0, ncol=4, nrow=region_end-region_start+1)
rownames(depth) <- pos
colnames(depth) <- c('ESC', 'CD45+', 'low pH', 'STAP-SC')
depth[esc$V1 - region_start + 1, 'ESC'] <- esc$V2
depth[cd45$V1 - region_start + 1, 'CD45+'] <- cd45$V2
depth[stapsc$V1 - region_start + 1, 'STAPSC'] <- stapsc$V2
depth[lowph$V1 - region_start + 1, 'low pH'] <- lowph$V2

# mm10 での Oct4 は mm9 よりも座標が手前なので注意
plot(0, 0, xlab='Position (chr17, mm10)', ylab='Depth', type='n', xlim=c(35490500, 35518000), ylim=c(0, 33))
lines(pos, depth[, 1], type='l', col=1)
lines(pos, depth[, 2], type='l', col=2)
lines(pos, depth[, 3], type='l', col=3)
lines(pos, depth[, 4], type='l', col=4)
legend("topleft", colnames(depth), col=1:4, lty=1, y.intersp=5)

これで、氏が UCSC genome browser にアップロードしているのと同じようなプロットができるはずだ。なお、kaho 氏が使っているマウスのリファレンス・ゲノムは mm9 であり、本稿で使っている mm10 とは遺伝子の座標が若干ズレていることに注意してほしい。

今回は両端のリード数はほぼ同じで、明らかに Oct4 周辺だけ飛び抜けているので結果は明らかだが、真面目にコピー数変化を議論するには、run 間でリード数が違う分なども考慮しなければならないはずである。

transgene のさらなる解析例 (3/7 追加)

kaho 氏の解析報告第3段とそのコメント欄に、CAG-GFP 細胞株が混ざっているのではないかという指摘があった。その配列を探してみよう。

CAG プロモータには、サイトメガロウィルスのエンハンサーやニワトリのβアクチンなど、マウスゲノムに存在しないはずの配列が含まれている、したがって CAG プロモータに由来するリードは現在 unmapped な状態になっているはずだ。

bwa がリファレンスに map できなかった配列は 4 というフラグが立つ。そこで、samtools view に -f 4 オプションをつけると map できなかったリードの抽出、-F 4 とすると map されたリードの抽出ができる。

samtools view -f4 -b ../map/SRR1171584.bam | samtools sort - lowph-unmapped
samtools index lowph-unmapped.bam

CAG配列は、http://www.cdb.riken.jp/pcs/protocol/vector/vector_top.html で提供されている。この配列はGenBank 形式であり、配列のほかにアノテーションを含んでいる。これを FASTA 形式にするには、GenBank to FASTA を利用すればよい。といっても、ORIGIN と書いてある行から始まるブロックを切り出して左端の塩基番号を削除し、一行目に配列名を > に続けて書くだけなので、テキストエディタawk でもすぐできる。

この配列に対して、上で抽出した unmapped リードをアライメントする。blatblast+ などいろいろなツールがあって適不適が異なるので、何を使うかは考慮すべきだが、今回は最初に用意した bwa で片付けてしまおう。やり方はゲノムへのマッピングと全く同じ。

bwa index -a bwtsw cag.fa

でインデックスをつくった後、アライメントして、

bwa aln -t 8 -b cag.fa lowph-unmapped.bam | samtools view -bS - | samtools sort - lowph-cag
samtools index lowph-cag.bam

CAG 配列にマッチした部分を抽出する。bwa は FASTQ のほか BAM ファイルも入力として利用できるのだ(-b オプション)。今回は結果を直接目で見てみよう。

samtools view -F4 lowph-cag.bam | wc -l
samtools view -F4 lowph-cag.bam | less -S

sam ファイルは、タブ区切りのテキストファイルで、リード名、フラグ(16はリファレンスと逆向きの鎖ということ)、配列名、マップした場所、マップの品質スコア、CIGARコード(アライメントのギャップなどを表す)……と続く。

という具合に、マッピングは簡単にできるのだが、ここからの解釈は難しい。

  • 操作に誤りはなかったか。ファイル名の取り違え、上書きなどは起こりやすい
  • ソフトウェアやパラメータの選択は適切か

という技術的な問題、

  • シーケンサーの読み間違い、アライメントのミスは確実に生じる。何本マッチしたら「有意に検出された」と言えるのか
  • 陽性・陰性対照を用意すべき

という統計的な問題、

  • どの配列が検出されたら CAG-GFP 株だと言えるのか。今回リファレンスにはベクターの配列全体を使ったが、ゲノムに組み込まれるのは全体ではない。
  • 逆に、CAG-GFP 以外の transgene コンストラクトにも共通して使われている配列もある

という生物学上の問題(分子生物学やトランスジェニック動物の知識が必要)、

  • 仮に CAG-GFP 株だったとして、論文の主張に重大な問題を来すのか

という論理上の問題がある。

データ処理の結果を解釈して何らかの結論を導くには、これらの点1つ1つをきちんと確認する必要がある。データ処理は始まりに過ぎないのだ。本記事で解釈に踏み込まないのは、私にはそれだけのスキルや知識もなければ、責任も持てないからである。結果の図を載せないのも、図だけが転載されて、これらの注意事項を離れて一人歩きし、騒ぎの材料に使われるのを危惧するからである。これらの点に留意した上で、各自、データ解析を楽しんでいただきたい。

*1:こちらで流行してるジャーゴン。必ずしも否定的なニュアンスではない。統計値だけでなくて、生データを見ることは大切