STARの挙動調査

by

いろんなマッピングツールがあって、バイオインフォやれば誰でも一つはマッピングツールを使うはず。私がマッピングツールを使う際に気にする点を列挙すると、

  • スプライスジャンクションを考慮するタイプかどうか　(for genome or transcriptome)
  • ショートリード向けか、ロングリード向けか
  • 同一スコアで複数箇所の候補がある場合、全部出力するのか、一つだけなのか、そもそも出力されないのか
    • 一つだけの場合は、どのようなランダムさで出力されるのか
  • 一つのリードの前半と後半で別の染色体にヒットする場合はどう出力されるのか

特に、複数箇所にヒットする場合の挙動は把握しておかないと、その後どのような解析を行うにしても間違った解釈をする危険がある。個人的にはゲノム解析をやることが多いので、ランダムに1箇所だけ出力してくれる挙動が好きだけど、トランスクリプトームをやる場合は下流の解析手法を考慮すると全部出力してほしいと思う。

下記のコマンドでSTARのインデックスを作成しておく。

/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --runMode genomeGenerate --genomeDir . --runThreadN 8 --genomeFastaFiles GRCh38.primary_assembly.genome.fa --sjdbGTFfile gencode.v32.annotation.gtf

まずはSTARはFASTQの中の順番に影響されて結果が変わるかどうかを調べる。seqkitを使って並びをシャッフルする。シード値を23, 55の2つ使って2通りのシャッフルしたファイルを作る。

seqkit shuffle -s 23 SRR7189551_100k_1.fastq > SRR7189551_100k_s23_1.fastq
seqkit shuffle -s 23 SRR7189551_100k_2.fastq > SRR7189551_100k_s23_2.fastq

seqkit shuffle -s 55 SRR7189551_100k_1.fastq > SRR7189551_100k_s55_1.fastq
seqkit shuffle -s 55 SRR7189551_100k_2.fastq > SRR7189551_100k_s55_2.fastq

できたファイルをマッピングして違いがあるか見てみる。

#マッピング
/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --genomeDir . --runThreadN 8 --outFileNamePrefix SRR7189551 --readFilesIn SRR7189551_100k_1.fastq SRR7189551_100k_2.fastq
/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --genomeDir . --runThreadN 8 --outFileNamePrefix SRR7189551_s23 --readFilesIn SRR7189551_100k_s23_1.fastq SRR7189551_100k_s23_2.fastq
/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --genomeDir . --runThreadN 8 --outFileNamePrefix SRR7189551_s55 --readFilesIn SRR7189551_100k_s55_1.fastq SRR7189551_100k_s55_2.fastq

#比較
diff <(cat SRR7189551Aligned.out.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n) <(cat SRR7189551_s23Aligned.out.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n)
diff <(cat SRR7189551_s23Aligned.out.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n) <(cat SRR7189551_s55Aligned.out.sam|grep -v "^@"|sort -k1,1V -k3,3 -k 4,4n)

入力ファイルがシャッフルされているので、出力をソートし直してからdiffで比較すると、完全一致していた。STARにはランダムさは全くなさそう。

次に複数箇所にマッピングされる場合どうなっているのか調べてみる。まずは入力のリード数は10 x 2=20万リード抜き出していたところで、マッピング結果のアライメント数を調べると

cat SRR7189551Aligned.out.sam|grep -v "^@"|wc -l
#223638 と表示される

アライメント数は22万なので、複数箇所を出力していそうである。一番ヒット数の多いリードを抜き出すと

cat SRR7189551Aligned.out.sam|grep -v "^@"|cut -f 1|uniq -c|sort -nr|head
     20 SRR7189551.41671
     20 SRR7189551.32557
     20 SRR7189551.24826
     20 SRR7189551.2234
     20 SRR7189551.19276
     20 SRR7189551.13760
     18 SRR7189551.99479
     18 SRR7189551.99196
     18 SRR7189551.95845
     18 SRR7189551.91358

となり、1つのペアエンドあたり10箇所まで出力されていそうではある。マニュアルで確認すると、デフォルトだと「–outFilterMultimapNmax」が「10」に設定されていて、10箇所まで出力されている。

ただ注意点として、「–outFilterMultimapNmax」を例えば「5」などに設定すると、上記のアライメントは一つも出力されない。

/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --genomeDir . --runThreadN 8 --outFileNamePrefix SRR7189551_n5 --outFilterMultimapNmax 5 --readFilesIn SRR7189551_100k_1.fastq SRR7189551_100k_2.fastq

cat SRR7189551_n5Aligned.out.sam|grep SRR7189551.41671
#何も表示されない。

つまり、STARはデフォルトだと10箇所以上の複数箇所にマッピングされてしまうリードは出力しない挙動になります。リピート配列が解析対象に入る場合や、small RNAを対象としているようなときは、outFilterMultimapNmaxを極端に大きくしておいたほうが良さそうです。

それから、STARはデフォルトではマッピングできなかったリードは出力しないようです。もし必要な場合は、「–outSAMunmapped Within KeepPairs」をつける必要があるようです。

個人的には、「–outSAMunmapped Within KeepPairs –outFilterMultimapNmax 1000」あたりを付けたほうが汎用的かなと思いますが、下流のツール次第ではありますよね。

コメントする