7/10 LAST コマンドミス修正
bwa memとLASTがナノポア向けにチューニングされたとナノポア公式ページでアナウンスされている。
https://nanoporetech.com/publications/bwa-and-last-have-been-tuned-work-nanopore-reads
bwa memはショートリード時代から1Mbpのリードのマッピングに対応していた有名なアライメントツールである。LASTはゲノム中の相同性の高い領域を探すために開発されたツールである。また、2016年にはGraphmapというロングリードのアライメントツールがnature communicationsに発表されている(ref.1)。この辺りを導入して、パフォーマンスを比較してみる。
LASTとbwaはbrewコマンドでmacに導入できる。（ない人だけ）。
brew install bwa
brew install LAST
１、bwa mem
bwa memを使ってnanopore readをリファンレスにマッピングするには、bwaのバージョン0.7.11から導入されたnanopore用のオプションを使うのが良いとされる。
bwa mem -x ont2d -t 20 ref.fa R1.fq R2.fq > nanopore.sam
- -t number of threads [1]
- -R read group header line such as '@RG id:foo SM:bar' [null]
-x ont2Dでナノポアリード向けのパラメータの一括調整を行なっている。このオプションをつけると-k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0と同じ意味になる。
- -k INT minimum seed length [19]
- -W INT discard a chain if seeded bases shorter than INT [0]
- -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [1.5]
- -A INT score for a sequence match, which scales options -TdBOELU unless overridden [1]
- -B INT penalty for a mismatch [4]
- -O INT[,INT] gap open penalties for deletions and insertions [6,6]
- -E INT[,INT] gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1]
- -L INT[,INT] penalty for 5'- and 3'-end clipping [5,5]
スコアの変更点を見るに、ミスマッチやギャップアライメントのペナルティーを下げ、最初のマッチングに使われるシード領域も短くしてアライメントが伸びるように工夫されたものと思われる。
Pacbioのリードは -x pacbioをつける。
詳細はこちらを確認してください。bwa memのパラメータチューニングをロングリードに向けどう変更していくかについて触れられています。Heng LiさんのBlogです。
BWA-MEM for long error-prone reads
2、LAST
LASTのチュートリアルHP LAST Tutorial
リファンレンスのindexを作る。
lastdb -uNEAR -R01 db ref.fa
- uNEAR selects a seeding scheme that makes it better at finding short-and-strong similarities. (It also changes the default scoring scheme.)
- -R01 単純リピートへのマッピング抑制
アライメントを行う。
lastal -Q1 -q1 -r1 -a1 -b1 -P0 db input.fq | last-split > myalns.maf
- -r: match score (2 if -M, else 6 if 0<Q<5, else 1 if DNA)
- -q: mismatch cost (3 if -M, else 18 if 0<Q<5, else 1 if DNA)
- -a: gap existence cost (DNA: 7, protein: 11, 0<Q<5: 21)
- -b: gap extension cost (DNA: 1, protein: 2, 0<Q<5: 9)
- -Q: input format: 0=fasta, 1=fastq-sanger, 2=fastq-solexa, 3=fastq-illumina, 4=prb, 5=PSSM (0)
- -P: number of parallel threads (1). #P0でmaximum
mafフォーマットからsamへの変換はmaf-convertを使う（マニュアル）。このツールでmafからaxt, blast, blasttab, html, psl, sam, tabへ変換できる（マニュアルリンク）。
maf-convert sam myalns.maf -r 'ID:X SN:Lambda_NEB LN:48502' > last-alignments.sam
リードグループも加えて変換する。ただし、できたsamのリードグループ情報２行目
@HD VN:1.3 SO:unsorted
@RG ID:X SN:Lambda_NEB LN:48502
のRGはSQに変えないとsam->bam変換でsamtools viewでエラーになってしまう（リードグループ）。vimかvi等で修正しておく。
アライメントのパラメーター設定であるが、bwaの戦略と方向性は同じで、ミスマッチペナルティやギャップペナルティを下げることでアライメント感度を高めている。特にgap existence costは標準の7-->1へと大きく下げている。このパラメータはONTリードをバクテリアにアライメントさせた論文（ref. 3）でも用いられている。
Both LAST alignment settings use a match score of 1 (-r1), gap opening penalty of 1 (option -a1) and a gap extension penalty of 1 (option -b1). Values of 1 (-q1) and 2 (-q2) were tried for the nucleotide substitution penalty parameters.
LASTの開発チームはONTリードをアライメントさせるパラメータについて特設ページを設けて考えを述べている。リンク=>
last-rna/last-long-reads.md at master · mcfrith/last-rna · GitHub
３、Graphmap
macにインストールするのは面倒そうだったので、cent OSにインストールした。
git clone https://github.com/isovic/graphmap.git
cd graphmap
make modules
make
make modulesでmake前に必要なものをダウンロードしている。
ビルドしてできた/graphmap/bin/Linux-x64/に本体が入っている。パスを通す。リファレンスへのアライメントは
graphmap align -t 12 -r ref.fa -d reads.fasta -o output.sam
-t: thread数
-d: シーケンスデータ（fasta、fastqに対応）。
-r: リファレンスファイル指定
-o: 出力
１−３の方法で作ったsamそれぞれをbamに変換してソート。
samtools view -@ 20 -bS input.sam |samtools sort -@ 20 -o sorted.bam
最初の方からパイプで繋いでも良い。
bamのindexをつける。
samtools index sorted.bam
検証
MNIONの1dモードでλコントロールを１時間ほど読んだデータを使う。データは
fasqに変換し、Porechopでトリムしてある。リファレンスのλゲノム配列 (47kb)はNCBIからダウンロードした。
０、アライメント状況
gifを貼ると重くなるので動画でリンクを貼る。λに貼り付けた時のエラーの出方はこんな感じになる(上半分がbwa mem、下半分がgraphmap）。
１、カバレッジ。
samtools depth dorted.bam > depth.txt
平均カバレッジ(read depth / bp)は以下のようになった。
- bwa mem (-x ont2dなし）1515 (1515)
- bwa mem (-x ont2dあり）1694 (1695)
- LAST 1772 (1774)
- graphmapデフォルト条件 1692 (1692)
( )はアダプタートリミング前のfastqを使った時の結果。
IGVでミスマッチの出方を見る。
Preferenceからtrack heightを300に変えて見やすくしている。上から、
- bwa mem (-x ont2dなし）
- bwa mem (-x ont2dあり）
- LAST
- graphmap
の順である。ぱっと見、一番ミスマッチが少なそうなのはgraphmapである。
２、mapされたリードの割合。
samファイルの３フィールド目の*(unmap)を数えて算出した。結果をまとめると
- bwa mem (-x ont2dなし） => 84.5 %（84.5）
- bwa mem (-x ont2dあり） => 82.3 %（82.3）
- LAST => 100 %（100 %）
- graphmap => 81.6 %（81.6）
( )はアダプタートリミング前のfastqを使った時の結果。
LASTはすべてのリードがλゲノムにアライメントされていた。ここで考えるべきは、LAST以外のマッパーでアライメントできなかった配列は何かである。
perlスクリプトを使いunmapのリードを抽出し、NCBIのblastにかけて調べてみることにした。まずmegablastとdiscontiguous megablastにかけたが、ヒットした配列が0になった。そこで断片的に似た配列を検索できるSomewhat similar sequences (blastn)を使ってみた。その結果、ヒットが出た配列は半分ほどになった。ただし、ほとんどがラムダに関係ないゲノムとの局所的なマッチだったので、Organismをλゲノムに限定してやり直した。そうすると、本数はわずかであるが、リードが部分的にλゲノムと一定の相同性を示すリードがあることがわかった。その画面を貼っておく。
黒いブロック１つは50bpほど。
赤いブロックのアライメントを貼っておく。
できればこうゆうリードはclippingしてアライメントして欲しいところ。
マッピングできなかったリードの大半は、よくわからないリードであった。信号がどこでおかしくなるのかは不明だが、fast5からの変換ツール次第で結果が変わる可能性もある。
アライメントされたリードの平均サイズ
Picardの CollectAlignmentSummaryMetrics を使う。
picard CollectAlignmentSummaryMetrics I=sorted.bam OUTPUT=aln_sum_metric.tsv
- bwa mem (-x ont2dなし） => 6073 bp（6107）
- bwa mem (-x ont2dあり） => 6073 bp（6107）
( )はアダプタートリミング前のfastqを使った時の結果。
graphmapとLASTで作ったbamはエラーになった。
Picardの代わりに、BAMstatsで４つのbamを比較する。順番は以下の通り、前と逆の順番にしてしまった。見にくくて申し訳ないです。
- graphmap
- LAST
- bwa mem (-x ont2dあり）
- bwa mem (-x ont2dあり）
アライメントされたリードの平均サイズ
マッピングクオリティ
今回試した中では、graphmapが一番マッピング率とエラー率の少なさのバランスが取れているように感じた (計算時間はbwa memよりは劣るが悪くはない）。
追記
他のマッパーも出てきています。
引用
1、Fast and sensitive mapping of nanopore sequencing reads with GraphMap
Ivan Sović, Mile Šikić, Andreas Wilm, Shannon Nicole Fenlon, Swaine Chen & Niranjan Nagarajan
Nature Communications 7, Article number: 11307 (2016) doi:10.1038/ncomms11307
2、LAST: Genome-Scale Sequence Comparison
LAST: genome-scale sequence comparison
３、A reference bacterial genome dataset generated on the MinION™ portable single-molecule nanopore sequencer
Joshua Quick,1,2 Aaron R Quinlan,3 and Nicholas J Loman Gigascience. 2014; 3: 22.
Published online 2014 Oct 20. doi: 10.1186/2047-217X-3-22
関連
