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

TopHat によるマッピング

RNA-Seq データ（*.fastq.gz）とゲノムの配列データ（Arabidopsis_thaliana.TAIR10.dna.toplevel.fa）を揃えたら、次にゲノムの配列データのインデックスを作成する。TopHat2 では Bowtie2 のインデックスを利用する。ここで Bowtie2 のコマンド bowtie2-build を使用してインデックスを作成する。インデックスの名前は自由に付けられるが、ここでは 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 になる場合がある。

次に TopHat2 を使って RNA-Seq データをゲノム上にマッピングする。サンプルが 8 個あるので、ここでは bash の for を使って、各サンプルに対して Bowtie2 マッピングを行う。

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

for seqlib in ${SEQLIBS[@]}; do
    tophat2 -p 4 -o ${seqlib}_tophat2 TAIR10 ${seqlib}.clean.fastq.gz
done

TopHat2 によるマッピング結果は -o で指定したフォルダ（例：SRR7508939_tophat2 など）に保存される。accepted_hits.bam がリードをゲノム配列にマッピングした結果である。マッピング率などは align_summary.txt に保存されている。

ls SRR444595_tophat2
## accepted_hits.bam  align_summary.txt  deletions.bed  insertions.bed  junctions.bed  logs  prep_reads.info  unmapped.bam

cat SRR444595_tophat2/align_summary.txt
## Reads:
##           Input     :  38088331
##            Mapped   :  33200505 (87.2% of input)
##             of these:   1377354 ( 4.1%) have multiple alignments (0 have >20)
## 87.2% overall read mapping rate.

ここで得られた BAM ファイル（accepted_hits.bam）に変換しソートした後に、 htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。SAM ファイルをソートしながら BAM ファイルに変換するには、samtools を利用する。

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

for seqlib in ${SEQLIBS[@]}; do
    samtools sort -O bam -o ${seqlib}.sorted.bam ${seqlib}.bam
    samtools index ${seqlib}.sorted.bam
done