Bowtie2 (A. thaliana, paired-end RNA-Seq) | Bowtie2 を利用した paired-end RNA-Seq データのマッピング

TopHat2 を使って、A. thaliana の paired-end RNA-Seq データを、ゲノム配列にマッピングしてみる。ここで使用するデータは PRJNA480638 でアーカイブされている。nsra nsrb double mutant と対照群それぞれ 3 サンプルからなる計 6 サンプルからなる。

RNA-Seq データ

データベースにアーカイブされた PRJNA480638 のデータを利用する。このデータは A. thaliana のデータで、nsra nsrb double mutant と対照群それぞれ 3 サンプルからなる計 6 サンプルからなる。サンプルは paired-end でシーケンシングされているので、リード情報が 1 サンプルに付き forward と reversed の 2 つのファイルに保存されている。この 6 サンプル 12 ファイルをローカルにダウンロードする。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR750/00${seqlib:9:10}/${seqlib}/${seqlib}_1.fastq.gz
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR750/00${seqlib:9:10}/${seqlib}/${seqlib}_2.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。フィルタリング後のファイルを *_1.clean.fastq.gz および *_2.clean.fastq.gz のような名前で保存する。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    java -jar trimmomatic-0.38.jar PE -phred33 -threads 4    \
    ${seqlib}_1.fastq.gz ${seqlib}_2.fastq.gz                \
    ${seqlib}_1.clean.fastq.gz ${seqlib}_1.unpaired.fastq.gz \
    ${seqlib}_2.clean.fastq.gz ${seqlib}_2.unpaired.fastq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:60
done

ゲノム配列データ

次に Ensembl Plants から A. thaliana のゲノム配列（dna.toplevel）をダウンロードする。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz

マッピングには必要としないが、マッピングされたリードを遺伝子領域ごとに集計するときにアノテーションファイルも必要である。このアノテーションファイルも（GFF または GTF）も合わせてダウンロードしておく。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.40.gff3.gz
gunzip Arabidopsis_thaliana.TAIR10.40.gff3.gz

Bowtie2 によるマッピング

RNA-Seq データ（*.fastq.gz）とゲノムの配列データ（Arabidopsis_thaliana.TAIR10.dna.toplevel.fa）を揃えたら、次にゲノムの配列データのインデックスを作成する。インデックスの名前は自由に付けられるが、ここでは TAIR10 と名付けている。ゲノムが大きいとき、1 時間近くかかる場合がある。

bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna.toplevel.fa TAIR10

インデックスを作成し終えると、そのフォルダには TAIR10.1.bt2、TAIR10.2.bt2、TAIR10.3.bt2、TAIR10.4.bt2、TAIR10.rev.1.bt2、TAIR10.rev.2.bt2 の 6 つのファイルが生成される。これらのファイルのプレフィックスが、インデックスの名前に指定した TAIR10 となっている。また、ヒトゲノムなどのように、ゲノムが大きい場合、インデックスの拡張子は .bt2 ではなく .bt2l になる場合がある。

次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。古いバージョンの Bowtie2 は圧縮された FASTQ ファイルに対応できないので、FASTQ を一度解凍してからマッピングする。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    gunzip ${sqelib}_1.clean.fastq.gz
    gunzip ${sqelib}_2.clean.fastq.gz
    bowtie2 -q -N 1 -p 4 -x HUMAN_INDEX -1${seqlib}_1.clean.fastq -2 ${seqlib}_2.clean.fastq -S ${sqelib}.sam
done

Bowtie2 のマッピング結果は SAM ファイルとして出力される。ここで得られた SAM ファイルを BAM ファイルに変換しソートした後に、 htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。

ls *.sam
## SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944