RNA-seq | トランスクリプトーム解析

Bowtie(シロイヌナズナ single-end)

Bowtie2を利用したマッピングの流れ

Bowtie はマッピングを行うプログラムで、リードが比較的に短い(〜 50bp)ときに利用される。1 個のリードがリファレンス上複数箇所にマッピングできるような場合にも対応できる。また、オプションを指定して、複数箇所にマッピングされるようなリードを無視するといったこともできる。

Bowtie を利用してリードをリファレンスゲノムへマッピングするには、おおよそ左図のような流れで行う。まず、公共データベースからリファレンスとなるゲノム配列をダウンロドして、インデックス化する。次に、FASTQ ファイル中のリードをリファレンスゲノムに対してマッピングを行う。マッピングするときに、インデックスを利用する。

このページでは、サンプルデータとして SRP011435 に登録されている 8 つの RNA-seq データを利用する。これらのデータはシロイヌナズナの RNA-seq データであり、single-end リードである。SRP011435 に登録されている 8 つの RNA-seq データの Accession Number はそれぞれ、SRR444595、SRR444596、SRR444597、SRR444598、SRR444599、SRR44460、SRR44461、SRR444602 である。これらのファイルを予めダウンロードして FASTQ 形式に変換しておく。

  1. 公共データベースからリファレンスゲノムをダウンロードする
  2. インデックスを作成する
  3. Bowtie2 を利用してマッピングを行う

公共データベースからリファレンスゲノムをダウンロードする

シロイヌナズナの全ゲノムデータは公共データベースである Ensembl Plants でダウンローできる。しかし、全ゲノムをダウンロードしても、どの領域が遺伝子でどの領域がそうではないかのアノテーション情報が含まれていない。そのため、(マッピングには必要はないが、)全ゲノムの他にアノテーションデータ(GFF または GTF)も合わせてダウンロードするとあとの解析が便利になる。

例えば、シロイヌナズナ関連のデータは Ensembl Plants FTP の Arabidopsis thaliana の項でダウンロードできる。

Ensembl FTP ダウンロードするファイル

上図の赤い枠で囲まれたファイルをダウンロードする。

  • Ensembl FTP ページのシロイヌナズナの項目から、DNA(FASTA) 列にある FASTA ファイルをダウンロードする。その他に、cDNA、CDS や Protein sequence など様々な FASTA ファイルが用意されているが、ダウンロードするのはあくまで DNA(FASTA) の FASTA である。該当リンクをクリックすると、複数のファイルがリストアップされる。このうち、.toplevel.fa.gz で終わるファイルをダウンロードする。このファイルが全ゲノム配列である。なお、※_sm.toplevel.fa.gz は利用しないので、ダウンロードする必要はない。
  • 同様に、Ensembl FTP ページの Gene sets 列にある GTF リンクをクリックすると、アノテーションデータのダウンロードページに移動できる。.gtf.gz で終わるファイルをダウンロードする。(このファイルは、マッピング時に利用しないが、後々の解析がスムーズに進めるように予めダウンロードしておくとよい。)

これらは gzip 形式で圧縮されているため、ダウンロード後にこれらを展開する。

gunzip Arabidopsis_thaliana.TAIR10.22.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.22.gtf.gz

展開すると、リファレンスとなる全ゲノムファイル(Arabidopsis_thaliana.TAIR10.22.dna.toplevel.fa)と アノテーションファイル (Arabidopsis_thaliana.TAIR10.22.gtf)が得られる。

インデックスを作成する

リファレンスとなる全ゲノムのインデックスを作成する。このとき、bowtie-build コマンドを利用する。-f の後に、ダウンロードした全ゲノムのファイル場所を指定し、その後にインデックスの名前を付ける。インデックス名前は任意で構わないが、ここでは TAIR10 とする。インデックス作成は非常に時間がかかる。パソコンの性能によっては数時間以上及ぶ場合がある。

bowtie-build -f Arabidopsis_thaliana.TAIR10.22.dna.toplevel.fa TAIR10

インデックスを作成し終えると、ディレクトリには 6 つのファイルが生成される。これら合わせて Bowtie のインデックスとなる。

  • TAIR10.1.ebwt
  • TAIR10.2.ebwt
  • TAIR10.3.ebwt
  • TAIR10.4.ebwt
  • TAIR10.rev.1.ebwt
  • TAIR10.rev.2.ebwt

Bowtie によるマッピング

Bowtie によるマッピングは以下のように bowtie コマンドを利用する。 サンプルデータは 8 つであるから、bowtie コマンドをそれぞれに対して実行する必要がある。

bowtie コマンドにオプション、インデックス、入力リードファイル、マッピング結果の出力ファイルの順で指定する。オプションに関しては、-v 2 はミスマッチを 2 個まで許容する、-p 4 は 4 スレッドを利用してマッピングを行う、-m 1 はリファレンス上の 1 カ所のみにマップされるリードを採用する、-S はマッピング結果を SAM 形式で出力することを意味する。

bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444595.fastq SRR444595.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444596.fastq SRR444596.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444597.fastq SRR444597.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444598.fastq SRR444598.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444599.fastq SRR444599.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444600.fastq SRR444600.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444601.fastq SRR444601.sam
bowtie -v 2 -p 4 -m 1 -S TAIR10 SRR444602.fastq SRR444602.sam

カウントデータの取得

次に HTSeq の htseq-count コマンドを利用してカウントデータの取得する。この際に、マッピング結果である SAM 形式のファイルと、アノテーションファイルである Arabidopsis_thaliana.TAIR10.22.gtf を必要とする。また、Arabidopsis_thaliana.TAIR10.22.gtf ファイルをそのまま HTSeq に与えると、エラーが生じるのでこれを Python スクリプトで変換する。ダウンロードした GTF ファイル中の9 列目のデータの値の部分に「;」が含まれているために、HTSeq ではこれをうまく取り扱えない。そのため、Python スクリプトで「;」を「_」に変換する。

python refine-gtf.py Arabidopsis_thaliana.TAIR10.22.gtf > Arabidopsis_thaliana.TAIR10.22.mod.gtf

次に、修正した GTF を利用してカウントデータの計算する。

htseq-count SRR444595.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444595.counts.txt
htseq-count SRR444596.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444596.counts.txt
htseq-count SRR444597.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444597.counts.txt
htseq-count SRR444598.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444598.counts.txt
htseq-count SRR444599.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444599.counts.txt
htseq-count SRR444600.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444600.counts.txt
htseq-count SRR444601.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444601.counts.txt
htseq-count SRR444602.sam Arabidopsis_thaliana.TAIR10.22.mod.gtf > SRR444602.counts.txt

実行が完了すると、SRR444595.counts.txt などのファイルが生成される。これらのファイルの中にカウントデータが記載されている。

カウントデータが複数のファイルにまたがっているため、これらを join コマンドで結合する。また、上述のカウントデータが記載されいてるファイルのなかに、遺伝子のカウントデータの他に、アノテーションがないリード数やマップされてなかったリード数なども保存されている。前者の場合は ENSG から始まり、後者の場合は「__no_feature」や「__not_aligned」などとなっている。そのため、grep 文を利用して ENSG だけから始まるデータのみを取得する。

# カウントデータのみを tmp.count.txt (一時ファイル)に保存する
join -j 1 SRR444595.count.txt SRR444596.count.txt | \
    join -j 1 - SRR444597.count.txt | \
    join -j 1 - SRR444598.count.txt | \
    join -j 1 - SRR444599.count.txt | \
    join -j 1 - SRR444600.count.txt | \
    join -j 1 - SRR444601.count.txt | \
    join -j 1 - SRR444602.count.txt > tmp.count.txt

# ヘッダーを all.count.txt に書き入れる
echo "EnsemblID    SRR444595    SRR444596    SRR444597    SRR444598    SRR444599    SRR444600    SRR444601    SRR444602" > all.count.txt

# カウントデータを all.count.txt に書き入れる
grep "^ENSG" tmp.count.txt >> all.count.txt

# 一時ファイルを削除する
rm tmp.count.txt