macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

Oxford Nanoporeリードのマッピング

 

  

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

 

1、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のパラメータチューニングをロングリードに向けどう変更していくかについて触れられています。bwaの中の人の記事と思われます。

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のリードグループ情報2行目

@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

 

 

 

 

3、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: 出力

 

  

1−3の方法で作ったsamそれぞれをbamに変換してソート。 

samtools view -@ 20 -bS input.sam |samtools sort -@ 20 -o sorted.bam

最初の方からパイプで繋いでも良い。

bamのindexをつける。

samtools index sorted.bam

 

 

 検証

MNIONの1dモードでλコントロールを1時間ほど読んだデータを使う。データは

fasqに変換し、Porechopでトリムしてある。リファレンスのλゲノム配列 (47kb)NCBIからダウンロードした。

 

0、アライメント状況

gifを貼ると重くなるので動画でリンクを貼る。λに貼り付けた時のエラーの出方はこんな感じになる(上半分がbwa mem、下半分がgraphmap)。

 

1、カバレッジ

samtools depth dorted.bam > depth.txt

平均カバレッジ(read depth / bp)は以下のようになった。

  1. bwa mem (-x ont2dなし1515  (1515)
  2. bwa mem (-x ont2dあり)1694  (1695)
  3. LAST           1772  (1774)
  4. graphmapデフォルト条件  1692 (1692)

( )はアダプタートリミング前のfastqを使った時の結果。

 

IGVでミスマッチの出方を見る。

f:id:kazumaxneo:20170620142250j:plain

Preferenceからtrack heightを300に変えて見やすくしている。上から、

  1. bwa mem (-x ont2dなし
  2. bwa mem (-x ont2dあり)
  3. LAST
  4. graphmap

の順である。ぱっと見、一番ミスマッチが少なそうなのはgraphmapである。

 

2、mapされたリードの割合。

samファイルの3フィールド目の*(unmap)を数えて算出した。結果をまとめると

  1. bwa mem (-x ont2dなし => 84.5 %(84.5)
  2. bwa mem (-x ont2dあり) => 82.3 %(82.3)
  3. LAST            => 100 %(100 %)
  4. graphmap         => 81.6 %(81.6)

( )はアダプタートリミング前のfastqを使った時の結果。

 

 

 

 

LASTはすべてのリードがλゲノムにアライメントされていた。ここで考えるべきは、LAST以外のマッパーでアライメントできなかった配列は何かである。

perlスクリプトを使いunmapのリードを抽出し、NCBIのblastにかけて調べてみることにした。まずmegablastとdiscontiguous megablastにかけたが、ヒットした配列が0になった。そこで断片的に似た配列を検索できるSomewhat similar sequences (blastn)を使ってみた。その結果、ヒットが出た配列は半分ほどになった。ただし、ほとんどがラムダに関係ないゲノムとの局所的なマッチだったので、Organismをλゲノムに限定してやり直した。そうすると、本数はわずかであるが、リードが部分的にλゲノムと一定の相同性を示すリードがあることがわかった。その画面を貼っておく。

f:id:kazumaxneo:20170620180259j:plain

 黒いブロック1つは50bpほど。

 

f:id:kazumaxneo:20170620180626j:plain

 赤いブロックのアライメントを貼っておく。

f:id:kazumaxneo:20170620180719j:plain

 できればこうゆうリードはclippingしてアライメントして欲しいところ。 

 

マッピングできなかったリードの大半は、よくわからないリードであった。信号がどこでおかしくなるのかは不明だが、fast5からの変換ツール次第で結果が変わる可能性もある。

 

 

 

 

アライメントされたリードの平均サイズ

 Picardの CollectAlignmentSummaryMetrics を使う。

picard CollectAlignmentSummaryMetrics I=sorted.bam OUTPUT=aln_sum_metric.tsv
  1. bwa mem (-x ont2dなし => 6073 bp(6107)
  2. bwa mem (-x ont2dあり) => 6073 bp(6107)

( )はアダプタートリミング前のfastqを使った時の結果。

graphmapとLASTで作ったbamはエラーになった。

 

Picardの代わりに、BAMstatsで4つのbamを比較する。順番は以下の通り、前と逆の順番である(見にくくて申し訳ないです)。

  1. graphmap
  2. LAST
  3. bwa mem (-x ont2dあり)
  4. bwa mem (-x ont2dあり)

アライメントされたリードの平均サイズ

f:id:kazumaxneo:20170620162746j:plain

 

カバレッジ

f:id:kazumaxneo:20170620162540j:plain

 

マッピングクオリティ

f:id:kazumaxneo:20170620163056j:plain

 

今回試した中では、graphmapが一番マッピング率とエラー率の少なさのバランスが取れているように感じた (計算時間はbwa memよりは劣るが悪くはない)。

 

 

 

bcftolosやBAMstatsなどの分析ツールは以下のエントリーでも紹介しています。

 

 

引用 

 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

 

 

3、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

A reference bacterial genome dataset generated on the MinION™ portable single-molecule nanopore sequencer