STARの2-pass mappingについて

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	=

コメントする