HISAT2 による single-end RNA-Seq データのマッピング

HISAT2 (A. thaliana, single-end RNA-Seq)

データベースにアーカイブされた 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

References

  • Huang T, López-Giráldez F, Townsend JP, Irish VF. RBE controls microRNA164 expression to effect floral organogenesis. Development. 2012, 139(12):2161-9. DOI: 10.1242/dev.075069
  • Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015, 12(4):357-60. DOI: 10.1038/nmeth.3317