ONTのロングリードのアセンブリとポリッシュ

2019 7/28 テスト結果追記

　2018年２月のNature CommunicationsにシロイヌナズナのゲノムをONTのロングリードを使ってアセンブリした論文（ PCR-free paired-end readsでpolish）が出ている（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5803254/）。ONTの MinION R9.4 flowcell (FLO-MIN106)で48時間シーケンスし、平均サイズ11.4kbのリードが合計3.4Gb得られたらしい（およそx30）。それからminiasmを使い一般的な4 Cores, 16 Gb RAMハードウエア（15インチのmacbook proのようなマシン）でアセンブルして、62 contigs（N50 12.3 Mb）得た、となっている。良好なパフォーマンスである。ただしエラーを許容するminiasmでアセンブリすると大量に誤りが残ってしまうため、著者らは、それからRaconを３回かけてコンセンサス配列を出力し、さらに最後にMiseqのPCR free ショートリード250x2を使ってPilonでポリッシュして、最終出力としている。十分なcontiguityのcontigが得られるため、オーサーはIntroductionの最後に、ロングリードでドラフトを作り構造変化（SV）を検出するのが一番シンプルでstraightforwardな手法であろうと主張してる。

Oxford nanoporeの公式レポジトリでも、このフローの精度が調べられている。

nanoporetech/ont-assembly-polish

https://github.com/nanoporetech/ont-assembly-polish

こちらも参考になります。

Assembly using miniasm+racon

http://inf-biox121.readthedocs.io/en/2016/Assembly/practicals/07_Assembly_using_minasm+racon.html

ONTのリードだけでアセンブルする手法は上記の論文以外に様々報告されており、精度を上げるアイデアも提案されている（リンク）。ここでは、上記の論文のようなminiasm（紹介）とracon（紹介）を使ったアセンブリとポリッシュのワークフローをまとめておく。Raconを複数回ランして精度が上がるように工夫している。

１、SRA_toolkitを使い上記のシーケンスデータ(ONT & illumina)をダウンロードする。

#ONTリードのダウンロード
prefetch ERR2173373 #このデータ
fastq-dump ~/ncbi/public/sra/ERR2173373.sra -O ONT #fastqに解凍

#illumina PCR free short read(250x20
prefetch ERR2173372 #このデータ
fastq-dump --split-files ~/ncbi/public/sra/ERR2173372.sra -O illumina #fastqに解凍

上記の論文では、この時点でナズナのリファレンスゲノムにリードをアライメントし、アライメントされなかったリードはlow qualityで、GC含量が大きく異なるようなリードであると見出している。

２、ダウンロードされたリードの分析

ショートリードはseqkit（紹介）で簡易チェック。

seqkit stats simulated.fa

$ seqkit stats ONT/ERR2173373.fastq

file format type num_seqs sum_len min_len avg_len max_len

ONT/ERR2173373.fastq FASTQ DNA 300,071 3,421,779,258 5 11,403.2 269,087

ロングリードはNanofiltを使いサイズとクオリティを分析（紹介）。

NanoPlot --fastq ONT/ERR2173373.fastq --loglength -t 12

f:id:kazumaxneo:20180324212546j:plain

f:id:kazumaxneo:20180324212549j:plain

数千bp付近に妙なピークがある。コピー数が多いクロロプラストの配列かもしれない（該当するサイズを抽出し、データベース検索して集計すれば調べられる）。

ここではアダプター除去や、クオリティトリミングは考えない。de novoアセンブリを実行する。

アセンブリ

３、minimap（紹介）でアセンブリ。

minimap -Sw5 -L100 -m0 -t 16 ONT/ERR2173373.fastq ONT/ERR2173373.fastq |gzip -1 > minimap_out.paf.gz #論文と同じパラメータ

miniasm -f ONT/ERR2173373.fastq minimap_out.paf.gz > miniasm_out.gfa

awk '/^S/{print ">"$2"\n"$3}' miniasm_out.gfa | fold > raw_assembly.fasta

トリミングなしでアセンブルすると、11,425,724bpのcontigが１つだけできた（染色体は5本ある、半数体ゲノムサイズはおよそ1.2億bp）。

ポリッシュ

４、Racon 1回目（Racon紹介）。Raconのランにはall reads vs all reads のマッピング部位を示したPAFファイル(Pairwise mApping Format) が必要なので、３でアセンブリしたcontigを使いminimapを動かす。

minimap raw_assembly.fasta ONT/ERR2173373.fastq > racon1.paf

#Raconでコンセンサス配列を出力。
racon -t 16 ONT/ERR2173373.fastq racon1.paf raw_assembly.fasta racon_polished1.fasta

５、Racon 2回目。

minimap racon_polished1.fasta ONT/ERR2173373.fastq > racon2.paf
racon -t 20 ONT/ERR2173373.fastq racon2.paf racon_polished1.fasta racon_polished2.fasta

６ --Racon 3回目。

minimap racon_polished2.fasta ONT/ERR2173373.fastq > racon3.paf
racon -t 16 ONT/ERR2173373.fastq racon3.paf racon_polished2.fasta racon_polished3.fasta

７、Pilon（紹介）を使いilluminaのハイクオリティリードでポリッシュする。まずbamを作る。

#マッピングとsort
bwa mem -t 16 racon_polished3.fasta illumina/ERR2173372_1.fastq illumina/ERR2173372_2.fastq |samtools sort -@ 16 -bS - > aln.bam

#index
samtools index aln_sorted.bam

８、ポリッシュを実行。

#pilonの実行ファイルをダウンロード
wget https://github.com/broadinstitute/pilon/releases/download/v1.22/pilon-1.22.jar

java -Xmx16G -jar pilon-1.22.jar --genome racon_polished3.fasta --frags aln_sorted.bam --vcf --tracks --changes --verbose --threads 16 --outdir output_dir

-----改装中------

結果