mrsFAST-UltraはSNPsに対応した次世代リードのアライメントツール。 mrsFASTの改良版となる。既知SNPsを許容しながら(ミスマッチとして扱わないためidentityが上がる)アライメントを行うことができる。indexファイルの軽量化にも成功しており、bowtie2でindexをつけるとヒトゲノムでは4GB近くの容量になってしまうが、mrsFAST-Ultraではおよそ2GBのindexで済む。また冗長なアライメントについて多様なオプションを持っており、bowtie2より高感度かつ高速(6倍)にアライメントを実行することができる。
公式ページ
mrsFAST: micro-read substitution-only Fast Alignment Search Tool
マニュアル
https://github.com/sfu-compbio/mrsfast//blob/master/README.md
インストール
https://github.com/sfu-compbio/mrsfast
git clone https://github.com/sfu-compbio/mrsfast
cd mrsfast
make
./mrsfast -h #ヘルプが表示されるか確認
mrsfastをパスの通ったディレクトリにコピーしておく。
ラン
サンプルデータをランする。
公式サイトからサンプルデータをダウンロード。解凍して中に入る。
http://sfu-compbio.github.io/mrsfast/
まずリファレンスのindexを作成する。
cd sfu-compbio-mrsfast-390c8eb/
mrsfast --index genome.fa
indexファイル(バイナリ)ができる。
アライメントを行う。シングルリード。
mrsfast --search genome.fa --seq reads.fq -o mapped.sam -u unmapped.fq --threads 8 -e 4
- -o Output the mapping record into file (default: output.sam)
- -u Output unmapped reads in file (default: output.nohit).
- --seq Set the input sequence to file. In paired-end mode, file should be used if the read sequences are interleaved. This file will be generated in all mapping mode.
- --threads 8 Use 1 number of cores for mapping the sequences (default: 1). Use 0 to use all the available cores in the system.
- -e error-threshold Allow up to error-threshold mismatches in the mappings.
デフォルトでは96 %同一の配列があればアライメントされる(-e 4)。
ペアリード。
mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8
- --seq1 Set the input sequence (left mate) to file. Paired-end option only.
- --seq2 Set the input sequence (right mate) to file. Paired-end option only.
- --pe Map the reads in Paired-End mode. min and max of template length will be calculated if not provided by corresponding options.
- --threads 1 Use 1 number of cores for mapping the sequences (default: 1). Use 0 to use all the available cores in the system.
インターレースのペアリードを指定するなら--seqを使う。
ペアリードのインサートサイズ(リードの外側からの距離)でアライメントを制限する。例えばインサートサイズが100以上3000以下に制限するなら
mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --min 100 --max 3000
100以下や3000異常など異常なインサートサイズの可能性が高いペアリード(discordant read pair)を排除できる。
minとmaxは--discordant-vhとともに使うと、正常なペアリードの範囲を指定するフラグとして働く。
mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --min 100 --max 1000 --discordant-vh
- --discordant-vh Map the reads in discordant fashion that can be processed by Variation Hunter / Common Law. Output will be generate in DIVET format.
- --max-discordant-cutoff m Allow m discordant mappings per read. Should be only used with --discordant-vh option.
ペアリードそれぞれについて先頭の40bpだけを使いアライメントを行う。
mrsfast --search genome.fa --seq1 mates1.fq --seq2 mates2.fq --pe -e 4 --threads 8 --crop 40
- --crop n Trim the reads to n base pairs from begining of the read.
5'末端 n-bpまでトリムではなく、3'末端を例えば50-bpをトリムする(捨てる)場合は--tail-crop 50とつける。
--その他--
- リファレンスにアライメントされたリードはsam形式で出力される。リファンレンスにアライメントされなかったリードはfastqのまま出力される。
- unmappedのfastqが必要なければ--disable-nohitsをつけることでunmapの出力を無効化できる。
- samのヘッダーが必要なければ--disable-sam-headerをつける。
- 巨大なfastqをアライメントする際は、--mem 4 などとつけることで、使用メモリー量を制限できる。この例では4GB以内でリードを読めるだけキャッシュし、それを繰り返すことでアライメントが進行する。
- gz圧縮されたfastqを入力とする時は、--seqcompをつける。
- 出力をgz圧縮したければ--outcompをつける。
- 1リードが複数アライメントされる場合、-nオプションでアライメントされる数を制限できる。-n 5なら最大で5カ所となる。
- --bestをつけるとベストマッピング(hamming distanceが最小)の箇所のみアライメントされる。
他にもdbSNPのデータベースを使い、既知SNPsでアライメントが正しく行えるようにするフラグも用意されています。その辺りは公式サイトを参考にしてください。
mrsfast/README.md at master · sfu-compbio/mrsfast · GitHub
引用
mrsFAST-Ultra: a compact, SNP-aware mapper for high performance sequencing applications
Faraz Hach Iman Sarrafi Farhad Hormozdiari Can Alkan Evan E. Eichler S. Cenk Sahinalp
Nucleic Acids Research, Volume 42, Issue W1, 1 July 2014, Pages W494–W500, https://doi.org/10.1093/nar/gku370
mrsFAST: a cache-oblivious algorithm for short-read mapping
Faraz Hach, Fereydoun Hormozdiari, Can Alkan, Farhad Hormozdiari, Inanc Birol, Evan E Eichler & S Cenk Sahinalp