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

解析手法

はじめに

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

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

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

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


使用するデータ

Pal, A., Leung, J. Y., Ang, G. C. K., Rao, V. K., Pignata, L., Lim, H. J., Hebrard, M., Chang, K. T., Lee, V. K., Guccione, E., & Taneja, R. (2020). EHMT2 epigenetically suppresses Wnt signaling and is a potential target in embryonal rhabdomyosarcoma. eLife, 9, e57683. https://doi.org/10.7554/eLife.57683

著者の実行環境
  • PC: MacBook Pro 14インチ、2021
  • チップ: Apple M1 Pro
  • メモリ: 32 GB
  • OS: macOS Sonoma 14.5
  • ストレージ: 512 GB
  • ターミナル: zsh

STARのインストール

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

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

コンパイル(compile)とは、プログラムをコンピュータが実行できる形式に変換することです。STARはC++言語で書かれており、C言語やC++言語などのプログラムは、実行する前にコンパイルする必要があります。

# カレントディレクトリの場所を確認
% 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.7.11b/bin/MacOSX_x86_64/」はエラーで動かすことができませんでした。修正版STARのコンパイルを試みるも、M1 Macではコンパイルできないようで、筆者は前のversionのSTAR2.7.10a.(STAR_MacOSX_x86_64.zip)をダウンロードして使用しました。

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

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のサイトにはリファレンスゲノムおよびアノテーションデータファイルが複数あり、それぞれ中身が少しづつ異なっているので、実験の目的に沿ったファイルを選択しましょう。今回選択したデータは下記の通りです。


リファレンスゲノム配列
  • ファイル保存名:Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
  • FASTA形式(.fa または .fasta)のデータで、配列情報が保存されている
  • gzip形式(.gz)で圧縮されている
  • GRCh38バージョンのヒトゲノム配列
  • ゲノム配列中の繰り返し領域(リピート領域)や低複雑性領域をマスキングしたリファレンスゲノムも公開されているが、今回はマスキングなしのリファレンスゲノムを選択した

アノテーションデータ
  • ファイル保存名:homo_sapiensHomo_sapiens.GRCh38.112.chr.gtf.gz
  • GTF(.ftf)形式のデータで、遺伝子のアノテーション情報が保存されている
  • gzip形式(.gz)で圧縮されている
  • GRCh38バージョンのヒトゲノムアセンブリに対する遺伝子アノテーションが含まれている
  • 遺伝子、トランスクリプト、エクソン、CDS(コーディングシーケンス)、UTR(非翻訳領域)などのアノテーションが含まれている

それでは、使用するデータを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時間、マッピングに約18時間かかりました。

# 作業するファイル(ディレクトリ)をホームディレクトリの下に作成
% 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

引数の説明
  • –runMode genomeGenerate:インデックスの作成モードを指定。
  • –genomeDir:インデックスを保存するディレクトリを指定。
  • –genomeFastaFiles:リファレンスゲノム配列のFASTAファイルを指定。
  • –sjdbGTFfile:アノテーションデータのGTFファイルを指定。
  • –runThreadN:使用するスレッド数(適宜変更する)を指定。
  • –sjdbOverhang:RNA-seqデータのリード長から1引いた値を指定(この値によってスプライスジャンクションの予測が改善する)。今回の場合、2×150 PE(ペアエンド)リードですので、150-1=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

引数の説明
  • –genomeDir:事前に作成されたリファレンスゲノムのインデックスディレクトリを指定。
  • –runThreadN:使用するスレッド数(適宜変更する)を指定。
  • –outFileNamePrefix:出力ファイル名の接頭辞を指定。
  • –quantMode: RSEMなどで発現量を算出するために、「TranscriptomeSAM」でBAM形式データの生成を指定。
  • –outSAMtype:出力するSAM/BAMファイルの形式を指定。「BAM SortedByCoordinate」とすると、BAM形式で出力し、リードを座標順にソートする。
  • –readFilesIn:マッピングするfastqファイルを指定。ペアエンドリードの場合は2つのファイルを指定する。
  • –readFilesCommand:入力ファイルを解凍するためのコマンドを指定。入力ファイルが圧縮されている場合、このコマンドを使って解凍する。今回の場合、fastqファイルはgzip形式で圧縮されていたため、「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 %:一意にマッピングされたリードの割合(マッピング率)を示す。
  • Mismatch rate per base, %:ミスマッチ率が低いほど、シーケンスデータの質が高く、正確にマッピングされていることを示す。
  • Average mapped length:マッピングされたリードの平均的な長さを示す。今回の場合、ペアエンドリード(PE150)なので、値が300bpに近いほどリードの大部分がリファレンスゲノムにうまくマッピングされていることを示す。
  • % of reads mapped to multiple loci:複数の位置にマッピングされたリードの割合が高いと、リファレンスゲノムの反復領域が多いことや、シーケンスの質が低いことを示す可能性がある。
  • % of reads unmapped: too many mismatches/too short/other:マッピングされなかったリードの割合が低いほど、シーケンスの質やリファレンスゲノムへのマッピングがうまくいっていることを示す。

この中でも特に 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」という遺伝子周辺を拡大しています。

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

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

3番目の「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遺伝子のオープンクロマチン化が促進され、その転写活性が増強される」と考察することができるでしょう。

ATAC-seqは染色体上のオープンクロマチン領域を特定する実験技術です。オープンクロマチン領域では、転写因子や他のタンパク質がDNAに容易に結合し、遺伝子の転写が活発に行われることが多いです。具体的には、ATAC-seqではトランスポゼースを使ってオープンクロマチン領域を切り出し、その領域のDNAシーケンスを行い、その位置をマッピングします。

まとめ

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