HISAT2 (T. aestivum, paired-end RNA-Seq) | HISAT2 を利用して小麦の paired-end RNA-Seq リードをマッピングする例

HISAT2（Kim et al, 2015）を使って、6 倍体のパン小麦（T. aestivum）の paired-end RNA-Seq データ（Ma et al, 2017）を、Chinese Spring (IWGSC RefSeq v1.0) のリファレンスゲノム（IWGSC, 2018）にマッピングしてみる。ここで、乾燥ストレス実験で取られた RNA-Seq データを利用する。データはは、PRJNA380841 でアーカイブされている。

RNA-Seq データ

データベースにアーカイブされた PRJNA380841 のデータを利用する。データは全部で 24 サンプルからなるが、HISAT2 の実行テストであれば、数サンプルをダウンロードすればよい。以下の例では、bash を使って 24 サンプル 48 ファイルをローカルにダウンロードしている。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR538/000/${seqlib}/${seqlib}_1.fastq.gz
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR538/000/${seqlib}/${seqlib}_2.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。フィルタリング後のファイルを *_1.clean.fastq.gz および *_2.clean.fastq.gz のような名前で保存する。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

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

ゲノムの配列データ

小麦（T. aestivum）Chinese Spring 系統のゲノム配列は URGI の FTP サイトで入手できる。wget コマンドで、ゲノム配列をダウンロードしてくる。ダウンロードした ZIP を展開すると、iwgsc_refseqv1.0_all_chromosomes フォルダが作られる。このフォルダの中に、小麦 Chinese Spring のゲノム配列が FASTA 形式で保存されている。

wget https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/iwgsc_refseqv1.0_all_chromosomes.zip
unzip iwgsc_refseqv1.0_all_chromosomes.zip
ls iwgsc_refseqv1.0_all_chromosomes
## 161010_Chinese_Spring_v1.0_pseudomolecules_AGP.tsv  161010_Chinese_Spring_v1.0_pseudomolecules_parts.fasta       161010_Chinese_Spring_v1.0_pseudomolecules_parts_to_scaffolds.bed  Agreement.txt
## 161010_Chinese_Spring_v1.0_pseudomolecules.fasta    161010_Chinese_Spring_v1.0_pseudomolecules_parts_to_chr.bed  161010_Chinese_Spring_v1.0_pseudomolecules_to_scaffolds.bed        README

また、マッピングには必要としないが、マッピングされたリードを遺伝子領域ごとに集計するときにアノテーションファイルも必要である。現在、Chinese Spring のゲノム配列 RefSeq v1.0 に対して、アノテーションファイルは v1.0 と v1.1 の 2 種類提供されているが、ここでは新しいバージョンをダウンロードしておく。ダウンロードした ZIP を展開すると、アノテーションファイル（GFF/GTF）が展開されたフォルダに保存されている。

wget https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.1/iwgsc_refseqv1.1_genes_2017July06.zip
unzip iwgsc_refseqv1.1_genes_2017July06.zip
ls iwgsc_refseqv1.1_genes_2017July06
## IWGSC_v1.1_HC_20170706_cds.fasta  IWGSC_v1.1_HC_20170706.idmapping.csv      IWGSC_v1.1_LC_20170706_cds.fasta  IWGSC_v1.1_LC_20170706_pep.fasta          IWGSC_v1.1_LC_removed.ids
## IWGSC_v1.1_HC_20170706.gff3       IWGSC_v1.1_HC_20170706_pep.fasta          IWGSC_v1.1_LC_20170706.gff3       IWGSC_v1.1_LC_20170706_transcripts.fasta  README.txt
## IWGSC_v1.1_HC_20170706.gtf        IWGSC_v1.1_HC_20170706_transcripts.fasta  IWGSC_v1.1_LC_20170706.gtf        IWGSC_v1.1_LC.correspondanceTEHC.txt

HISAT2 によるマッピング

まず hisat2-build を使用して、小麦のゲノム配列のインデックスを作成する。ここではインデックスの名前を IGWSCv1 としてインデックスを作成する。上でダウンロードしたゲノム配列には、161010_Chinese_Spring_v1.0_pseudomolecules.fasta と 161010_Chinese_Spring_v1.0_pseudomolecules_parts.fasta の 2 種類存在する。前者のゲノム配列は、HISAT2 や STAR などのマッピングリファレンスとして使用することが推奨されていない。原因として、アセンブリ配列に含まれている 1 contig （1 染色体）の塩基数が長すぎて、マッピング後の処理が正常に行われなかったりする。そのため、マッピング用として、後者の方の FASTA を利用する。また、小麦のゲノム配列が非常に大きいので、インデックスの作成に非常に時間がかかる。

hisat2-build 161010_Chinese_Spring_v1.0_pseudomolecules_parts.fa IGWSCv1

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

次に、作成したインデックスを利用してリードをゲノム配列にマッピングしていく。マッピング率などは画面に出力される。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

for seqlib in ${SEQLIBS[@]}; do
    hisat2 -x IGWSCv1 -1 ${seqlib}_1.clean.fastq.gz -2 ${seqlib}_2.clean.fastq.gz -p 4 -S ${seqlib}.sam 
done
## 24729208 reads; of these:
##   24729208 (100.00%) were paired; of these:
##     2515871 (10.17%) aligned concordantly 0 times
##     20025563 (80.98%) aligned concordantly exactly 1 time
##     2187774 (8.85%) aligned concordantly >1 times
##     ----
##     2515871 pairs aligned concordantly 0 times; of these:
##       226564 (9.01%) aligned discordantly 1 time
##     ----
##     2289307 pairs aligned 0 times concordantly or discordantly; of these:
##       4578614 mates make up the pairs; of these:
##         2625062 (57.33%) aligned 0 times
##         1541246 (33.66%) aligned exactly 1 time
##         412306 (9.01%) aligned >1 times
## 94.69% overall alignment rate
## ...

HISAT2 が出力しているログを確認すると、uniquely mapped リードが全体の 80.98% であった。6 倍体小麦で同一遺伝子のコピー数が多いにもかかわらず、uniquely mapped リードの割合が大きかった。

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

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

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 SRR5383650.sorted.bam | grep 'NH:i:1' | wc -l

References

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
Ma J, Li R, Wang H, Li D, Wang X, Zhang Y, Zhen W, Duan H, Yan G, Li Y. Transcriptomics Analyses Reveal Wheat Responses to Drought Stress during Reproductive Stages under Field Conditions. Front Plant Sci. 2017, 8:592. DOI: 10.3389/fpls.2017.00592
IWGSC. Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science. 2016, 361(6403):eaar7191. DOI: 10.1126/science.aar7191