生命情報科学

STAR による RNA-Seq リードの高速マッピング

STAR

STAR は RNA-Seq リードをリファレンス配列にマッピングするプログラムである。STAR が発表された当時、他のマッピングプログラムに比べ、STAR のマッピング速度は非常に高速であると知られていた。現在でも、マッピングプログラムの中で高速な方に分類される。

STAR のマッピングアルゴリズムは、 seed searching と clustering / stitching / scoring の 2 ステップからなる（Dobin et al, 2013）。seed searching のステップにおいて、あるリードが与えられると、STAR はそのリードの先端の塩基から、リファレンス配列と照合していき、もっとも長くマッピングできる場所を探す。論文中では、リードの塩基配列のうち、リファレンス配列とマッピングできた部分を MMP (Maximal Mappable Prefixes) と呼んでいる。もしリード全体がリファレンス配列にマッピングできなかった場合、マッピングできなかった部分（MMP でない部分）を取り出して、同様な作業により再度マッピングしていく。このように、仮に 1 つのリードが 2 つのエクソンから構成されていても、STAR はこれを正確にマッピングできる。また、リードとリファレンス配列の間に、ミスマッチが存在しする場合、STAR はミスマッチを容認してマッピングすることができる。さらに、リードの末端部分にアダプター配列、ポリ A 配列、あるいは品質が極端に低い配列などが含まれる場合、STAR はその直前でマッピングを止めて、リードの残りの部分をトリムする。

マッピングにより、リード上の MMPs が決定されると、次に clustering / stitching / scoring のステップにおいて、各リードの MMPs を評価していく。例えば、あるリードの最初の MMP とそれに続く MMP の間の距離が指定された長さ（maximum intron size）以下かどうかをチェックしたり、ミスマッチやギャップ数が指定された数以下かどうかをチェックしたりして、最終的にそのリードのマッピング位置を決定する。また、STAR は、1 つのリードの先端と後半について、その距離が長すぎたり、両者が異なる染色体にマッピングされたりあるいは異なる DNA 鎖にマッピングされたりするキメラリードもつけることができる。

STAR の使い方

STAR を使うとき、リファレンス配列からインデックスを作成し、RNA-Seq データをインデックスに対してマッピングするという 2 ステップからなる。

インデックスの作成

リファレンス配列からインデックスを作成するとき、--runMode genomeGenerate で STAR を実行する。リファレンス配列は、--genomeFastaFiles を利用して FASTA ファイルを指定する。作成するインデックスの名前は --genomeDir で指定する。STAR はメモリを大量に消費するので、メモリエラーが起こる場合、メモリサイズ --limitGenomeGenerateRAM を増やす必要がある。

mkdir index_name
STAR --runMode genomeGenerate \
     --genomeDir index_name \
     --genomeFastaFiles genome.fa \
     --limitGenomeGenerateRAM 3400000000

マッピング

single-end リード

STAR --runThreadN 4 \
     --genomeDir index_name \
     --readFilesIn SRR0000.fastq \
     --genomeLoad NoSharedMemory \
     --outFilterMultimapNmax 1

paired-end リード

STAR --runThreadN 4 \
     --genomeDir index_name \
     --readFilesIn SRR0000_1.fastq SRR)))_2.fastq \
     --readFilesCommand gunzip -c \
     --genomeLoad NoSharedMemory \
     --outFilterMultimapNmax 1

STAR のマッピングオプション

STAR ソフトウェアには多くのオプションが用意されている。ここでは、よく使えそうなものだけをリストアップした。完全なオプションは STAR の配布サイト（GitHub）から入手できる。

--readFilesCommand gunzip -c FASTQ ファイルが圧縮されている場合、このオプションを与えることによって、STAR は自動的に解凍しながらファイルを読み込む。*.gz の場合は --readFilesCommand gunzip -c、*.bz2 の場合は --readFilesCommand bunzip2 -c のように使用する。
--outFilterMultimapNmax 10 1 つにリードがマッピングできる位置の最大値。指定された数値よりも多くの箇所にマッピングされる場合は、そのリードをマッピングできなかったとみなす。Uniquely mapped reads のみを解析対象としたい場合は、--outFilterMultimapNmax 1 と指定する。何も指定しないときは 10。
--outFilterMismatchNmax 許容できる最大ミスマッチ数。マッピング結果に指定された数よりも多くのミスマッチを含む場合、そのリードをマッピングできなかったとみなす。何も指定しないときは 10。
--outFilterMismatchNoverLmax 0.3 マッピングできた部分の長さに含まれるミスマッチの許容率。マッピングできた部分の長さに含まれるミスマッチの割合が指定された数値を超える場合、そのリードをマッピングできなかったとみなす。何も指定しないときは 0.3。
--outFilterMismatchNoverReadLmax 1.0 マッピング結果に含まれるミスマッチのリード長に対する許容率。マッピング結果に含まれているミスマッチの、リード全長に対する割合が、指定された割合を超える場合、そのリードをマッピングできなかったとみなす。何も指定しないときは 1.0。
--chimOutType
--chimOutType SeparateSAMold キメラリードの出力ファイルを指定する。キメラリードマッピング結果専用のファイルに保存するには SeparateSAMold を指定し、他の通常のマッピング結果と同じファイルに保存するには WithinBAM を指定する。デフォルトは SeparateSAMold。
--chimSegmentMin 0 キメラリードの場合、リードの前半と後半で異なる染色体などにマッピングされる。リードの前半あるいは後半の長さのうち短い方が、指定された長さ以上のとき、そのリードをキメラリードとみなして出力する。paired-end リードの場合は、ペアを合計して考える。例えば --chimSegmentMin 20 のとき、2 x 75 bp の paired-end リードの forward リードの 75 bp と reversed リードの最初の 55 bp が連続的にマッピングされ、reversed リードの後半 20 bp が異なる染色体にマッピングされた場合、これをキメラリードとみなす。前半（forward 75 bp + reversed 55 bp）のマッピングクオリティが十分に高い場合、これは通常のマッピング結果としてもみなされる。つまり、1 つのリードが通常のマッピング結果とキメラリードのマッピング結果の両方に含まれる場合がある。また、例えば reversed リードの後半の 15 bp だけが異なる染色体にマッピングされた場合は、これをキメラリードと見なさない。デフォルトは 0 で、つまりキメラリードを検出しない。

STAR を用いた RNA-Seq マッピング例

SATR マッピングログ

STAR によるマッピングログは、Log.final.out ファイル中に記述される。このログファイルの中で「Uniquely mapped reads number」の項目には uniquely mapped リードの数が記載されている。例えば、次の例では、uniquley mapped リードは 12928073 本あることになる。STAR では、paired-end リードを 1 つのリードとして扱うため、この Uniquely mapped reads number には paired-end リードを 1 リードとして数えられている。また、Uniquely mapped reads number には paired-end のうち一方だけが uniquely mapped したリードも含まれている。

                          Number of input reads |	26305337
                      Average input read length |	233
                                    UNIQUE READS:
                   Uniquely mapped reads number |	12928073
                        Uniquely mapped reads % |	49.15%
                          Average mapped length |	217.10
                       Number of splices: Total |	9813105
            Number of splices: Annotated (sjdb) |	9213614
                       Number of splices: GT/AG |	9681255
                       Number of splices: GC/AG |	110284

STAR が出力する SAM/BAM ファイルから Uniquely mapped reads number を計算するためには、次のように A および B の数を先に計算し、次に A/2 + B で計算する（参考）。

samtools view -c -q 255 -f 3 Aligned.out.refsort.bam
## A

samtools view -c -q 255  -f 9 -F 4  Aligned.out.refsort.bam
## B

STAR では uniquely mapped リードの MAPQ を 255 としているため、-q 255 で uniquely mapped リードだけを出力できるようにする。また、-f 3 で paired-end mapped リードを取得でき、-f 9 -F 4 で paired-end のうちどちらか一方だけが uniquely mapped したリードを取得できる。

References

  • Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013, 29(1):15-21. DOI: 10.1093/bioinformatics/bts635
  • Uniquely mapped reads number. rnastar (Google groups)