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.gzTopHat2 によるマッピング
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 になる場合がある。
次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。
SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)
for seqlib in ${SEQLIBS[@]}; do
tophat2 -p 4 -o ${seqlib}_tophat2 TAIR10 ${seqlib}_1.clean.fastq.gz ${seqlib}_2.clean.fastq.gz
doneTopHat2 によるマッピング結果は -o で指定したフォルダ（例：SRR7508939_tophat2 など）に保存される。リードをゲノム配列にマッピングした結果は accepted_hits.bam ファイルに保存されている。また、マッピング率などの情報は align_summary.txt ファイルに保存されている。
ls SRR7508939_tophat2
## accepted_hits.bam align_summary.txt deletions.bed insertions.bed junctions.bed logs prep_reads.info unmapped.bam
cat SRR7508939_tophat2/align_summary.txt
## Left reads:
## Input : 26526105
## Mapped : 26472702 (99.8% of input)
## of these: 851888 ( 3.2%) have multiple alignments (1788 have >20)
## Right reads:
## Input : 26526105
## Mapped : 26461614 (99.8% of input)
## of these: 851272 ( 3.2%) have multiple alignments (1794 have >20)
## 99.8% overall read mapping rate.
##
## Aligned pairs: 26411814
## of these: 849434 ( 3.2%) have multiple alignments
## 207876 ( 0.8%) are discordant alignments
## 98.8% concordant pair alignment rate.ここで得られた BAM ファイル（accepted_hits.bam）に変換しソートした後に、 htseq-count または featureCounts などのプログラムに与えれば、遺伝子のリードカウントデータを得ることができる。