超初心者向け！！RNA-seq解析シリーズ④HISAT2でマッピングする

誰でもできるRNA-seq解析シリーズ！
今回はRNA-seq解析のメインとも言えるマッピングを行なっていきます！

HISAT2のインストール

まずはHISAT2をインストールします。
すでにhomebrewをインストールしましたので、
簡単です。
brew install プログラム名でインストールできます。
installの前後は半角スペースです。

brew tap brewsci/bio
brew install hisat2

condaを使っている場合は、bioconda チャンネルから入れられます。

conda install -c bioconda hisat2

問題なくインストールされたら、以下のコマンドを入力してみましょう。

hisat2 -h

HISAT2のヘルプが表示されたらOKです。
もし、hisat2 command not foundが出たら、
インストールがうまくできていません。

リファレンスゲノムの取得

HISAT2のダウンロードが完了したら、次はリファレンスに用いるゲノム配列を用意します。
今回はマウスのRNA-seqデータでしたので、マウスのゲノム配列が必要になります。
各自解析するRNA-seqデータの動物種を予め確認しましょう。

HISAT2のサイトからリファレンスゲノムをダウンロード

さて、今回はマウスのリファレンスゲノムを取得していきます。
HISAT2でマッピングするには、リファレンスゲノムを取得した後、
マッピングの速度を上げるために、リファレンスゲノムをindex化する必要があります。

よく使われるリファレンスについては、HISAT2のサイトにindex化されたものが既にありますので、
そこからindex済のリファレンスをダウンロードしていきます。

外部リンクHISAT2 公式ダウンロードページ（index一覧） ↗https://daehwankimlab.github.io/hisat2/download/

その中から、今回はマウスのリファレンスゲノムとしてよく使われるmm10をダウンロードします。

外部リンクmm10 index ダウンロードはココから ↗https://genome-idx.s3.amazonaws.com/hisat/mm10_genome.tar.gz

ダウンロード後に、tar -zxvfコマンドでtar.gzを解凍します。

#tar.gzファイルの解凍
tar -zxvf mm10_genome.tar.gz

解凍すると mm10_genome というフォルダができ、その中に genome.1.ht2 〜 genome.8.ht2 の8個のファイルが入っています。この8個セットで1つのindex で、HISAT2に渡すときは拡張子や番号を除いた mm10_genome/genome という「共通の頭の部分」を指定します。

自分でindexを作る場合（hisat2-build）

公式サイトにindexが用意されているのは、ヒト・マウスなど主要な生物種だけです。それ以外の生物種を解析する場合は、ゲノムFASTAから自分でindexを作ります。

# genome.fa からindexを作成（出力の頭の部分を genome と指定）
hisat2-build -p 4 genome.fa genome

実行すると genome.1.ht2 〜 genome.8.ht2 が生成されます。あとはダウンロードしたindexと同じように -x genome で指定できます。

スプライス部位を考慮したindex（公式サイトで _tran と付いているもの）を作りたい場合は、アノテーションGTFから情報を抽出してから渡します。

hisat2_extract_splice_sites.py annotation.gtf > splicesites.txt
hisat2_extract_exons.py annotation.gtf > exons.txt
hisat2-build -p 4 --ss splicesites.txt --exon exons.txt genome.fa genome_tran

マッピング

では、マッピングしていきます。
その前に、先ほど作成されたリファレンスゲノムのフォルダの位置を移動させます。
リファレンスゲノムのフォルダ（mm10_genome）と、解析するRNA-seqデータを同じフォルダに入れてください。

そして、ターミナルの作業ディレクトリをRNA-seqとリファレンスの入ったディレクトリに移動させましょう。

#ディレクトリの移動
cd /RNA-seqデータとリファレンスゲノムフォルダのあるディレクトリ
ls
mm10_genome       SRR1571967_1.fastq.bz2       SRR1571967_2.fastq.bz2

lsコマンドでmm10_genome, SRR1571967_1.fastq.bz2, SRR1571967_2.fastq.bz2がある状態ならOK！

マッピングの実行

以下のコマンドを実行し、マッピングを行います。

hisat2 -t -p 4 -x mm10_genome/genome  -1 SRR1571967_1.fastq.bz2 -2 SRR1571967_2.fastq.bz2 -S SRR1571967.sam

コマンドの説明

hisat2 hisat2の実行を命令
-t 時間経過を表示
-p 4 使用するCPUのスレッド数
-x mm10_genome/genome indexの「共通の頭の部分」を指定（.1.ht2 などの拡張子は付けない）
-1 SRR1571967_1.fastq.bz2 pair-endの１つ目のファイルを指定
-2 SRR1571967_2.fastq.bz2 pair-endの2つ目のファイルを指定
-S SRR1571967.sam samファイルで出力することと出力ファイル名を指定

結果の確認

#マッピング終了後
Time loading forward index: 00:00:07
Time loading reference: 00:00:01
Multiseed full-index search: 00:44:05
34523088 reads; of these:
  34523088 (100.00%) were paired; of these:
    3020396 (8.75%) aligned concordantly 0 times
    28397973 (82.26%) aligned concordantly exactly 1 time
    3104719 (8.99%) aligned concordantly >1 times
    ----
    3020396 pairs aligned concordantly 0 times; of these:
      147825 (4.89%) aligned discordantly 1 time
    ----
    2872571 pairs aligned 0 times concordantly or discordantly; of these:
      5745142 mates make up the pairs; of these:
        3508209 (61.06%) aligned 0 times
        1930162 (33.60%) aligned exactly 1 time
        306771 (5.34%) aligned >1 times
94.92% overall alignment rate
Time searching: 00:44:07
Overall time: 00:44:14

SRR1571967のマッピングの結果です。
最後の overall alignment rate がマッピング率で、この例では94.92% です。これを全てのサンプルに対して同様に実行します。

よくある質問

Q. hisat2 command not found になる インストールが完了していないか、パスが通っていません。Homebrewの場合は brew tap brewsci/bio を先に実行したか確認してください。tapを忘れると brew install hisat2 自体が失敗します。condaで入れた場合は、目的の環境を conda activate してから実行しているかを確認します。which hisat2 でパスが表示されるかが切り分けの起点です。

Q. condaでHISAT2を入れるには？ conda install -c bioconda hisat2 で入ります。依存関係の解決に失敗する場合は、チャンネルの優先順位を設定してから試してください。

conda config --add channels bioconda
conda config --add channels conda-forge
conda install hisat2

Q. indexは自分で作る必要がある？ ヒト・マウスなど主要な生物種であれば、公式サイトの作成済みindexをダウンロードするだけで済みます。index作成は時間もメモリも消費する ので、あるものは使うのが基本です。公式に無い生物種の場合のみ hisat2-build で作成します。

Q. -x には何を指定すればいい？ フォルダ名でもファイル名でもなく、indexファイルの「共通の頭の部分」 を指定します。mm10_genome/genome.1.ht2 〜 genome.8.ht2 があるなら -x mm10_genome/genome です。.1.ht2 まで書いたり、フォルダ名だけを書いたりするとエラーになります。

Q. genome と genome_tran はどう違う？ genome はゲノム配列だけから作ったindex、genome_tran はアノテーションのスプライス部位・エクソン情報を組み込んだindexです。既知のスプライスジャンクションをまたぐリードのマッピング精度が上がるため、アノテーションが整備された生物種では _tran が有利です。ただしファイルサイズが大きくなります。

Q. ヒトゲノム（GRCh38）で解析したい 手順はマウスと同じで、ダウンロードするindexを差し替えるだけです。公式ダウンロードページの GRCh38 の項目から genome または genome_tran を取得し、-x grch38/genome のように指定します。後の工程で使うアノテーションGTFも、必ず同じアセンブリ（GRCh38）のものに揃えてください。ゲノムとアノテーションのバージョンがずれると座標が合いません。

Q. マッピング率が低いときはどこを見る？ まず生物種とアセンブリが合っているかを疑います。次にアダプター配列やクオリティの低い塩基が残っていないか（トリミング前のデータではないか）、そしてstrand指定やpair-endの -1 -2 の指定が正しいかを確認します。rRNAが大量に混入している場合もマッピング率とは別に下流の結果が歪むので、overall alignment rate だけでなく、後の工程で得られる発現量分布も見るようにします。

次回は、StringTieの使い方を説明いたします。

関連記事Nextflowとnf-core/rnaseqでbulk RNA-seq解析を自動化する：QC・STAR・Salmonを再現可能に実行 →

関連記事

関連記事超初心者向け！！RNA-seq解析シリーズ③公共データベースからRNA-seqデータをダウンロード →

関連記事超初心者向け！！RNA-seq解析シリーズ⑤StringTieの使い方 →