データベースにアーカイブされた A. thaliana の single-end RNA-Seq データ（PRJNA153493）をサンプルとして、HISAT2（Kim et al, 2015）で 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 の配列に対応している。

HISAT2 によるマッピング

リファレンスとなる全ゲノムのインデックスを作成する。このとき、hisat2-build コマンドを利用する。インデックス名前は任意で構わないが、ここでは TAIR10 とする。インデックス作成は非常に時間がかかる。パソコンの性能によっては数時間及ぶ場合がある。

hisat2-build Arabidopsis_thaliana.TAIR10.22.dna.toplevel.fa TAIR10

インデックスを作成し終えると、そのフォルダには TAIR10.*.ht2（* は 1 から 8 の数字）のようなファイルが 8 個生成される。これらのファイルは HISAT2 が使うインデックスとなる。

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

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

for seqlib in ${SEQLIBS[@]}; do
    hisat2 -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

HISAT2 によるマッピング後、マッピング結果として SAM ファイルが出力される。samtools を使ってこの SAM を BAM ファイルに変換して、ソートする。ソートした BAM を 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

samtools で BAM ファイルあるいは SAM ファイルから uniquely にマッピングできたリード数を数えるには次のようにする。（ただし、paired-end の場合は、-f 0x2 で paired-end リードを取得してから grep する必要がある。）

samtools view -h SRR444595.sorted.bam | grep 'NH:i:1' | wc -l