【RNA-seq解析 Part 3】マッピング（アライメント）のやり方について詳しく解説！

はじめに

前回はFATSQファイルの前処理・クオリティチェックを行いました。今回は下記の記事の続きとなります。

リードの配列と既知の配列（リファレンス配列）を比較し、リードがどこのDNA領域に由来するかを特定するプロセスをマッピングまたはアライメントと言います（図1）。現在、RNA-seqのマッピングには、HISAT2またはSTARのどちらかが使用されている場合が多いです。STARは大量のメモリを消費しますが、HISAT2と比較してマッピング率、マッピング速度、正確性に優れるという報告があり、今回はSTARを用いたマッピング方法について詳しく解説していきます。

HISAT2も広く使用されており、使用することに全く問題はないですが、「パソコンのスペックが低い場合にHISAT2の使用する」という選択が良いかと思います。

マッピングが完了すると、BAM形式のファイルが生成されます。このBAMファイルを用いて発現量を計算します。発現量の計算については、次回の記事で紹介します。

STARのインストール

STARをインストールします。最新版は、STARの公式サイトからダウンロードしてください。STARの詳しい使い方は、公式のドキュメントに記載がありますのでそちらを参照ください。

公式サイトでは自分でコンパイルする方法も書かれていますが、今回は「STAR-2.7.11b/bin/MacOSX_x86_64/」に保存されたコンパイル済みの実行ファイルを使用します。

# カレントディレクトリの場所を確認
% pwd

# カレントディレクトリにSTARをインストール
% wget https://github.com/alexdobin/STAR/archive/2.7.11b.tar.gz
% tar -xzf 2.7.11b.tar.gz

# 解凍したディレクトリに移動し、STARバイナリを適切な場所に配置
# ここにファイルを置くことで、システム全体でSTARコマンドを実行できるようになる
% sudo cp STAR /usr/local/bin

# バイナリの存在を確認
% ls /usr/local/bin | grep STAR

#バージョンの確認
% STAR --version

「#バージョンの確認」でSTARのバージョンが確認できればインストールは完了です。

インデックスファイルの作成

STARは高速なマッピングを実現するために、事前にゲノムのインデックスを作成する必要があります。STAR用のインデックスを作成するために必要なファイルは下記の2つで、これらのファイルを使用して、STARはゲノムのインデックスを作成します。

• リファレンスゲノム配列（FASTA形式）: ゲノム配列のFASTA形式のファイル。
• アノテーションデータ（GTFまたはGFF形式）: 遺伝子やトランスクリプトの位置、構造、属性などの情報を含むアノテーションデータ。一般的には、GTF（Gene Transfer Format）またはGFF（General Feature Format）形式のファイルが使用される。

リファレンスゲノムとアノテーションデータはいくつかのサイトで入手することができますが、今回はEnsembl（https://asia.ensembl.org/Homo_sapiens/Info/Index）というサイトからダウンロードします。Ensemblのサイトにはリファレンスゲノムおよびアノテーションデータファイルが複数あり、それぞれ中身が少しづつ異なっているので、実験の目的に沿ったファイルを選択しましょう。今回選択したデータは下記の通りです。

それでは、使用するデータをwgetコマンドでダウンロードします。

# 作業するファイル（ディレクトリ）をホームディレクトリの下に作成
% mkdir /Users/usrname/GSE142975
% cd /Users/usrname/GSE142975

# カレントディレクトリの場所を確認
% pwd

# カレントディレクトリにリファレンス配列とアノテーションデータをダウンロード
# リファレンスゲノム配列のダウンロード
% wget https://ftp.ensembl.org/pub/release-112/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz        

# アノテーションデータのダウンロード
% wget https://ftp.ensembl.org/pub/release-112/gtf/homo_sapiens/Homo_sapiens.GRCh38.112.chr.gtf.gz

# gunzipコマンドで、gzip形式で保存されているデータを解凍
% gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
% gunzip Homo_sapiens.GRCh38.112.chr.gtf.gz

続いて、STAR用のインデックスファイルを作成します。STARの詳細な使い方は、公式サイトで確認してください。インデックスファイルはゲノムの配列情報を効率的に格納し、検索可能な形式に変換します。これにより、STARは迅速にリードをゲノムにマッピング（アライメント）することができます。「finished successfully」と表示されたら、インデックスファイルの作成完了です。

筆者の環境ではインデックスファイルの作成に約2時間、マッピングに約5時間かかりました。

# 作業するファイル（ディレクトリ）をホームディレクトリの下に作成
% mkdir /Users/usrname/GSE142975
% cd /Users/usrname/GSE142975

# カレントディレクトリの場所を確認
% pwd

# STAR用インデックスファイルの保存先（ディレクトリ）を作成
% mkdir star_index

# STAR用インデックス作成
% STAR \
     --runMode genomeGenerate \
     --runThreadN 8 \
     --genomeDir /Users/usrname/GSE142975/star_index \
     --genomeFastaFiles /Users/usrname/GSE142975/Homo_sapiens.GRCh38.dna.primary_assembly.fa \                   
     --sjdbGTFfile /Users/usrname/GSE142975/Homo_sapiens.GRCh38.112.chr.gtf \
     --sjdbOverhang 149

リファレンス配列へのマッピング（アライメント）

マッピングを行います。マッピングでは、インデックスを利用して、各リードをリファレンスゲノム上の最も適切な位置に配置していきます。

マッピングが完了すると、BAMファイル（.bam）とSAMファイル（.sam）が生成されます。BAMファイルはマッピングされたリードの配列情報をバイナリ形式で圧縮して保存します。一方、SAMファイルは同じ情報をテキスト形式で保持しており、人間が読み取りやすい特性がありますが、ファイルサイズが大きくなります。この解析ではSAMファイルは使用しないため、生成を省略します。

BAMファイルは、可視化、変異検出、発現量の計算などに使用していきます。

# STARでマッピングを行う
% STAR \
     --genomeDir /Users/usrname/GSE142975/star_index \
     --runThreadN 8 \
     --outFileNamePrefix SRR10822355 \
     --quantMode TranscriptomeSAM \
     --outSAMtype BAM SortedByCoordinate \
     --readFilesIn /Users/usrname/GSE142975/SRR10822355_filt_1.fastq.gz       /Users/usrname/GSE142975/SRR10822355_filt_2.fastq.gz \
     --readFilesCommand gunzip -c

マッピングが終了すると、出力ファイルの一つに「SRR10822355Log.final.out」があります。このファイルには、マッピング実行結果の最終的な統計情報が記録されており、このファイルを確認することでマッピングのパフォーマンスや結果を確認することができます。

SRR10822355のマッピング結果（一部を抜粋）は以下の通りでした。

  % cat SRR10822355Log.final.out

# 結果を一部抜粋して表示しています
Uniquely mapped reads % |	95.32%
Mismatch rate per base, % |	0.27%
Average mapped length |	295.68
of reads mapped to multiple loci |	2.51%
of reads unmapped: too many mismatches |	0.00%
of reads unmapped: too short |	2.00%
of reads unmapped: other |	0.13%

この中でも特に Uniquely mapped reads %（マッピング率）をチェックすることが多いです。どの程度のマッピング率であれば良いかという明確な評価基準はありませんが、今回の場合は「95.32%」という高いマッピング率であり、マッピングに問題はないと考えられます。また、他の評価指標からも、マッピングに問題なさそうです。

マッピングの可視化

最後に、IGVというツールを用いて、マッピング結果を可視化します。

% cat SRR10822355Log.final.out
# samtoolsのインストール
% brew install samtools

# IGVのインストール
conda install -c bioconda -y igv 
conda install -c bioconda -y igvtools

# BAMファイルのインデックスファイル（.bai）の作成
samtools index SRR10822355Aligned.sortedByCoord.out.bam

# IGVを開き、マッピングを可視化する
igv -g hg38 SRR10822355Aligned.sortedByCoord.out.bam

今回は、解析に使用している論文で注目している「EHMT2」という遺伝子周辺を拡大しています。

１番目の「SRR10822355Aligned.sortedByCoord.out.bam Coverage」は、ゲノム上のリードのカバレッジ（深さ）を示します。特定のゲノム領域にマッピングされたリードの数をグラフ化して表示しています。

２番目の「SRR10822355Aligned.sortedByCoord.out.bam Junctions」は、リードのジャンクション（スプライスジャンクション）を示します。これは、RNA-seqデータなどで特に重要です。

３番目の「SRR10822355Aligned.sortedByCoord.out.bam」は、BAM形式で保存されたリードのマッピング（アライメント）結果を示します。具体的には、リードの配列がゲノム上のどの位置にマッピングされたかを表示します。

このように、IGVではマッピングの結果を詳細に確認することができます。

最後に、IGVによる解析例（Xu Z et al., Front Nutr, 2022から抜粋）を紹介し、その結果から考察できることを考えてみます。

この論文では、IMF（筋肉内脂肪）含有量の異なる豚の背最長筋に対してATAC-seqとRNA-seqを実施し、HOXA9遺伝子周辺の領域をIGVで確認しています。

ATAC-seqのマッピング結果を見ると、高IMF群ではHOXA9遺伝子領域がオープンクロマチンになっていることがわかります。特にHOXA9遺伝子の上流で、高IMF群は顕著なオープンクロマチン領域があります。

一方、RNA-seqのマッピング結果からは、高IMF群でHOXA9遺伝子の発現が増加し、低IMF群ではその発現がほとんどないことがわかります。これにより、オープンクロマチンの形成とHOXA9遺伝子の発現量との間に相関があることが示唆されます。

以上の結果から、「IMF含有量が高いとHOXA9遺伝子のオープンクロマチン化が促進され、その転写活性が増強される」と考察することができるでしょう。

まとめ

• RNA-seqのマッピング（アライメント）は、リードをリファレンスゲノムに配置するプロセスであり、HISAT2またはSTARが主に使用される。
• BAMファイルはマッピングにより生成され、リードの位置情報などのマッピング結果が含まれる。
• BAMファイルはマッピング結果の可視化、変異検出、発現量計算などに使用する。
• SAMファイルは、人間が読み取りやすいようにBAMファイルをテキスト形式に変換したもので、サイズが大きいため通常は生成しない。
• IGVを使用すると、BAMファイルからマッピング結果を可視化することができる。