今回はHISAT2について調べてみる。HISAT2は最もよく使われるRNA-seqのマッピングツールなのではないかと思う。
インデックスの作成は、
hisat2-build -p 8 GRCh38.primary_assembly.genome.fa GRCh38.primary_assembly.genome.fa のようにすれば良いのだけど、ちゃんとしたインデックスを作るために、 –snp, –ss, –exon といったオプションを付けると、数百GBのメモリが必要となる。HISAT2のHPにインデックスが存在するなら、それを使うほうが無難かも。
マッピングは前回のSTARの時と同じように、まずは入力順をランダムにシャッフルしたFASTQファイルを3つ使ってみる。
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_1.fastq -2 SRR7189551_100k_2.fastq -S test1.sam
#結果の一部
100000 reads; of these:
100000 (100.00%) were paired; of these:
6737 (6.74%) aligned concordantly 0 times
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_s23_1.fastq -2 SRR7189551_100k_s23_2.fastq -S test2.sam
#結果の一部
100000 reads; of these:
100000 (100.00%) were paired; of these:
6708 (6.71%) aligned concordantly 0 times
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_s55_1.fastq -2 SRR7189551_100k_s55_2.fastq -S test3.sam
#結果の一部
100000 reads; of these:
100000 (100.00%) were paired; of these:
6718 (6.72%) aligned concordantly 0 timesマッピングの結果は、上記の通りFASTQの順番を入れ替えただけなのだけど、なぜかマッピングできないリード数が変化している。念のため、同じ入力ファイルでもう一度実行してみる。
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_1.fastq -2 SRR7189551_100k_2.fastq -S test1-2.sam
#結果の一部
100000 reads; of these:
100000 (100.00%) were paired; of these:
6737 (6.74%) aligned concordantly 0 timesすると、大まかな結果は同じで、ちゃんと中を比較しても同じ結果が得られていた。HISAT2は同じ入力ファイルであれば、同じ結果を返すようだ。
入力ファイルをシャッフルした場合の影響を詳しく調べてみるために、出力のSAMをリード名順に並び変えて、FASTQをシャッフルした影響のあったリード数を調べてみた。
diff <(cat test1.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n|cut -f 1-6) <(cat test2.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n|cut -f 1-6)|grep "<"|cut -f 1|sort|uniq|wc -lすると、10万リード中約5,000リードほど影響があった。そもそもマッピングされた場所が違うリードも多く、違っていてもマッピングクオリティは最大の60がつけられていた。HISAT2を使う場合、スコアが最大でも、本当にマッピングされた場所由来のリードなのかというと、FASTQの入力順で変わるくらいなので、あまり信頼しないほうが良いのかも？
さて、出力されるリードのヒット数の最大値を調べると
cat test1.sam|grep -v "^@"|cut -f 1|uniq -c|sort -nr|head
10 SRR7189551.99968
10 SRR7189551.99945
10 SRR7189551.99891
10 SRR7189551.99833
10 SRR7189551.99545
10 SRR7189551.994795 x 2 (paired-end) = 10ヒットとなって、5 hitが上限の様子。マニュアルで見てみると、-kと–max-seedsあたりが上限を決めている感じ。そこで、それらのオプションを大きい値(1万)にしてみる。
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_1.fastq -2 SRR7189551_100k_2.fastq -S test1-3-k10000-max-seeds10000.sam -k 10000 --max-seeds 10000
cat test1-3-k10000-max-seeds10000.sam|grep -v "^@"|cut -f 1|uniq -c|sort -nr|head
20000 SRR7189551.78592
20000 SRR7189551.32156
1826 SRR7189551.42636
1752 SRR7189551.49875
1732 SRR7189551.4850
1672 SRR7189551.39053
すると、上限よりも少ないヒット数が多いのだけど、まだ1万か所以上にマッピングされるリードがある模様。HISAT2のマニュアルには、kを大きく設定できるように設計されていないと書かれており、確かに10万などに設定すると100倍以上は遅くなった。
hisat2 -q -p 8 --dta -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_1.fastq -2 SRR7189551_100k_2.fastq -S test1-3-k100000-max-seeds100000.sam -k 100000 --max-seeds 100000
cat test1-3-k100000-max-seeds100000.sam|grep -v "^@"|cut -f 1|uniq -c|sort -nr|head
120448 SRR7189551.78592
30080 SRR7189551.32156
1826 SRR7189551.42636一応、10万個以下のヒットで収まっている様子。
ここまでの結果を見ると、HISAT2は同一スコアの中からランダムに-kの数だけ出力する挙動に見えます（ただし本当に最善の場所にマッピングされているのかは不明）。リピート配列を対象とする場合、マニュアルに書かれていますが、HISAT2は不向きです。
最後にStringTie用の–dtaオプションを外した場合、
hisat2 -q -p 8 -x GRCh38.primary_assembly.genome.fa -1 SRR7189551_100k_1.fastq -2 SRR7189551_100k_2.fastq -S test1-4-nodta.sam
100000 reads; of these:
100000 (100.00%) were paired; of these:
2781 (2.78%) aligned concordantly 0 times何やらマッピングできない数が3％ほど減りました。diffを取って比べてみると…
diff <(cat test1.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n|cut -f 1-6) <(cat test1-4-nodta.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n|cut -f 1-6)|grep "<"|cut -f 1|sort|uniq|wc -l10万リード中約1万リードが影響を受けるようです。HISAT2のマニュアルには、「このオプションを使用すると、HISAT2はスプライスサイトの新規発見のために長いアンカー長を必要とします」とあり、確かに–dtaを付けた場合にはイントロンをまたいだマッピングが減っているようではあります。