HISAT2（Kim et al, 2015）を使って、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

HISAT2 によるマッピング

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

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

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

次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。マッピング率などは画面に出力される。以下に 1 つ目のサンプルの出力例を示した。この出力結果を見ると、リードペアの 97.49% がゲノム上どこか一箇所だけにマッピングされる uniquely mapped reads であり、1.04% が複数箇所にマッピングされるリードであった。全体として 98.53% のリードがマッピングされた。

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

for seqlib in ${SEQLIBS[@]}; do
    hisat2 -x TAIR10 -1 ${seqlib}_1.clean.fastq.gz -2 ${seqlib}_2.clean.fastq.gz -p 4 -S ${seqlib}.sam 
done
## 26526105 reads; of these:
##   26526105 (100.00%) were paired; of these:
##     390228 (1.47%) aligned concordantly 0 times
##     25860774 (97.49%) aligned concordantly exactly 1 time
##     275103 (1.04%) aligned concordantly >1 times
##     ----
##     390228 pairs aligned concordantly 0 times; of these:
##       350005 (89.69%) aligned discordantly 1 time
##     ----
##     40223 pairs aligned 0 times concordantly or discordantly; of these:
##       80446 mates make up the pairs; of these:
##         7484 (9.30%) aligned 0 times
##         64289 (79.92%) aligned exactly 1 time
##         8673 (10.78%) aligned >1 times
## 99.99% overall alignment rate

HISAT2 によるマッピング結果は SAM ファイルに保存される。

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

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

SEQLIBS=(SRR7508939 SRR7508940 SRR7508941 SRR7508942 SRR7508943 SRR7508944)

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 にマッピングできたリード数を数えるには次のようにする。samtools view では -f 0x2 オプションで paired-end リードの両方がマッピングされているリードを取得できる。また、HISAT2 では uniquely mapped リードには NH:i:1 が定義されているので、これを含む行を取得すれば、uniquely mapped リードを取得できるようになる。

samtools view -h -f 0x2 SRR7508939.sorted.bam | grep 'NH:i:1' | wc -l