データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ（PRJNA153493）をサンプルとして、Bowtie2 で TAIR10 のゲノム配列へマッピングする。このデータセットは 8 サンプルを含み、4 サンプルが DEX-treated 35S で、他の 4 サンプルは mock-treated 35S である（Huang et al, 2012）。

RNA-Seq データ

データベースにアーカイブされた single-end RNA-Seq （PRJNA153493） のデータを ENA からダウンロードする。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR444/${seqlib}/${seqlib}.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)

for seqlib in ${SEQLIBS[@]}; do
    java -jar trimmomatic-0.38.jar SE -phred33 -threads 4  \
    ${seqlib}.fastq.gz ${seqlib}.clean.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

• Ensembl にはゲノムの全配列データの他に、cDNA（coding DNA）や CDS（cDNA から UTR などを取り除いたデータ）などの配列データも提供されている。ここでは、ゲノムの全配列をダウンロードしたいので、DNA (.toplevel.fa.gz) をダウンロードしてくる。
• ゲノムの配列データとして *_sm.toplevel.fa.gz も配布されているが、これを利用しないので、ダウンロードする必要がない。
• Ensembl にゲノム配列に対応したアノテーションデータ（GFF/GTF）も配布されている。このアノテーションデータは .toplevel.fa.gz の配列に対応している。

Bowtie2 によるマッピング

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

bowtie2-build -f Arabidopsis_thaliana.TAIR10.22.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 を使って RNA-Seq データをゲノム上にマッピングする。サンプルが 8 個あるので、ここでは bash の for を使って、各サンプルに対して Bowtie2 マッピングを行う。

SEQLIBS=(SRR444595 SRR444596 SRR444597 SRR444598 SRR444599 SRR444600 SRR444601 SRR444602)

for seqlib in ${SEQLIBS[@]}; do
    gunzip ${seqlib}.clean.fastq.gz
    bowtie2 -q -N 1 -p 8 -x TAIR10 -U ${seqlib}.clean.fastq -S ${seqlib}.sam
done

ls *.sam
## SRR444595.sam SRR444596.sam SRR444597.sam SRR444598.sam SRR444599.sam
## SRR444600.sam SRR444601.sam SRR444602.sam

Bowtie2 によるマッピング後、マッピング結果として SAM ファイルが出力される。samtools を使ってこの SAM を BAM ファイルに変換して、ソートする。ソートした BAM を featureCounts などのプログラムに与えれば、各遺伝子領域にマッピングされたリード数を数えることができる。