GATK ver3時代にRNA-seqのデータからSNPをコールする際にSTARの2-passを使っていたときは、インデックスを2回作らないといけなかったのだけど、今マニュアルを見ていたら2.4.1a以降はマッピング時のパラメータに「–sjdbFileChrStartEnd /path/to/sj1.tab」を追加するだけで良いとなっている。早速2-pass mappingの効果を見てみる。マッピングして、比較を下記のように行う。
#mapping
/suikou/tool/STAR-2.7.3a/bin/Linux_x86_64_static/STAR --genomeDir . --runThreadN 8 --outFileNamePrefix SRR7189551_unmapped_2pass --outSAMunmapped Within KeepPairs --sjdbFileChrStartEnd SRR7189551_unmappedSJ.out.tab --readFilesIn SRR7189551_100k_1.fastq SRR7189551_100k_2.fastq
#比較
diff <(grep -v "^@" SRR7189551_unmappedAligned.out.sam|cut -f 1-7) <(grep -v "^@" SRR7189551_unmapped_2passAligned.out.sam|cut -f 1-7) |more
その結果、例えば次のような箇所で違いがある。上が1-passで、下が2-passの結果。5列目に例えば23041Nなどと(数字)+Nが追加されているところは、23,041塩基のイントロンをまたいでマッピングできたという意味。イントロンをまたいだマッピングは困難なので、2-passにするとより精度良くイントロンをまたいだマッピングができているのがわかる。
#read SRR7189551.91
< SRR7189551.91 147 chr16 81077136 255 12S88M =
---
> SRR7189551.91 147 chr16 81054083 255 12M23041N88M =
#read SRR7189551.2673
< SRR7189551.2673 147 chr7 51156252 255 2S98M =
---
> SRR7189551.2673 147 chr7 51136329 255 3M19921N97M =
#read SRR7189551.5612
< SRR7189551.5612 163 chr20 17659132 255 1S99M =
---
> SRR7189551.5612 163 chr20 17659072 3 1S16M30N83M =