macでインフォマティクス

macでインフォマティクス

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

ナノポアのロングリードを使うとアセンブルはどのくらい改善されるのか?

ハイクオリティなショートリードのデータに、ロングリード情報を混ぜ込むとどれくらいアセンブリは改善されるのか調べてみる。 

 

 

NがあってもgrepUCSC ゲノムブラウザで除くことができるが、それでは肝心の繰り返し領域の評価が曖昧になる。やはりNがないコンプリートな真核ゲノムを使いたいということで、100M以下の生き物で調べたところ、シゾンが完全にFinishしていた。そこでまずはシゾンのゲノムを例にアセンブルのパフォーマンスを調べてみる。Ensemblからfastaファイルをダウンロードする。

f:id:kazumaxneo:20171006194024j:plain

chromosomeが20あり、ゲノムサイズは16.7Mbである。クロロプラストゲノムとミトコンドリアゲノムは今回使わない。

 

 

ショートリードのシミュレーションデータ

Miseq v3をシミュレートできるARTを使う。

シェルスクリプト "Q30fatq_generation.sh"

# ./workflow.sh #シェルスクリプトでfastq作成、trimmingまで自動処理
#実行権をつけ、genome.fastaのあるパスに移動して./create.shでラン



#ランは "./Q30fatq_generation.sh"
cov='250'
insert='600'
size='210'
stdev='60'

#ループ処理 カレントディレクトリの’ファイル’のみを取得
for file in `\find * -maxdepth 1 -type f`; do
a=${file%.fasta}_raw_paired #rename
folder=${file%.fasta} #folder_name

mkdir $folder
#fastqを発生 後でサンプリングするため、多めに250x作る。(リンク
art_illumina -ss MSv3 -sam -i $file -p -l $size -f $cov -m $insert -s $stdev -o $folder/$a
fastqF=$folder/${a}1.fq
fastqR=$folder/${a}2.fq

#quality trimming(リンク
after.py -1 $fastqF -2 $fastqR -q 30 -s 25 -f 50 --no_overlap
mv good/* $folder
fastqFHQ=${fastqF%fq}good.fq
fastqRHQ=${fastqR%fq}good.fq
done

Miseq v3のプロファイルだと先頭35bpほどqaulityが低い。上のコードでは5'側を35-bp強制トリミングしている。

 

平均リードサイズが求まったので、ゲノムサイズから考えて必要なカバレッジになるようリードをランダムサンプリング。

seqkit sample -p $sampling $fastqFHQ > $folder/R1.fastq
seqkit sample -p $sampling $fastqRHQ > $folder/R2.fastq

ここではx100取り出す。

 $samplingは計算で求めた0~1の範囲の少数。ここでは分けて書いているが、実際は1つのシェルスクリプトにしてカバレッジ100でQ30以上のリードを取り出すまで自動処理した。

 

phread quality score 30以上のハイクオリティなリードが手に入った。リード長を確認。

user$ seqkit stats *fastq

file      format  type   num_seqs      sum_len  min_len  avg_len  max_len

R1.fastq  FASTQ   DNA   5,067,191  810,750,560      160      160      160

R2.fastq  FASTQ   DNA   5,067,191  810,750,560      160      160      160

250-bp発生させたが、QT後は平均160-bpまで縮んでいる。大きく減ったのは、illumina Miseq V3プロファイルでは5'側のアダプター直後が低クオリティなのが要因として大きい。

 

 

 

ナノポアのロングリードのシュミレーションデータ

OXTの2Dのデータ(R9.5)をシミュレートしてみる(リンク)。

wget ftp://ftp.bcgsc.ca/supplementary/NanoSim/yeast* #yeastのデータをダウンロード
read_analysis.py -i yeast_2D.fasta -r yeast_S288C_ref.fa #2Dのデータをもとにプロファイル作成
simulator.py linear -r Cyanidioschyzon_merolae.ASM9120v1.dna.toplevel.fasta -c training -n 60000 #しゾンゲノムをもとにロングリードを発生させる。リード数6万

やや多いが6万リード発生させた。 リード情報を確認。

user$ seqkit stats simulated_reads.fasta 

file                   format  type  num_seqs      sum_len  min_len  avg_len  max_len

simulated_reads.fasta  FASTA   DNA     60,000  288,693,949       35  4,811.6  214,880

平均4800-bp、max210-kbとなった。クオリティは調べていないが、2Dデータをもとに発生させているのでエラー率は10%以下と思われる。カバレッジはx17となる。

 

 

 De novo assembly

最近色々な方法論が発表されているが、まずはspadesでアセンブルしてみる (version 3.11)。

ショートリードのみアセンブル

spades.py -t 40 -k auto --careful -1 R1.fastq -2 R2.fastq -o short

 ショート+ロング。

spades.py -t 40 -k auto --careful -1 R1.fastq -2 R2.fastq --nanopore simulated_reads.fasta -o hybrid

miniasmでロングリードのみアセンブル

minimap -Sw5 -L100 -m0 -t8 simulated_reads.fasta simulated_reads.fasta | gzip -1 > reads.paf.gz
miniasm -f pacbio_filtered.fastq reads.paf.gz > reads.gfa
awk '/^S/{print ">"$2"\n"$3}' reads.gfa | fold > reads.fa

polishが手間なので今回は行なっていない。ショートリードについては、spadesからUnicyclerに切り替えてもよかった。今回は時間の関係でspadesにした。

 

Results

アセンブル結果はこんな感じ。長さtop50。

f:id:kazumaxneo:20171006205248j:plain

long read単独でもショートリードよりアセンブルは長くなっている。

 

bwa memを使い、contigをゲノムにアライメントしてIGVで様子を眺めてみる。FInishしたゲノムがなければできない方法だが、こうすることで切断部位がどうなっているか掴むことができる場合がある。

シェルスクリプト bwa-mem.shを各fastaに対して実行。

genome='Cyanidioschyzon_merolae.ASM9120v1.dna.toplevel.fasta'
contig='short_scaffolds.fasta'

bwa index -a is $genome
bwa mem -t 20 $genome $contig |samtools view -S -b - > short.bam
samtools sort -@ 20 short.bam > short_sorted.bam
samtools index short_sorted.bam

IGVを起動。

igv -g Cyanidioschyzon_merolae.ASM9120v1.dna.toplevel.fasta short_sorted.bam,hybrid_sorted.bam,long_sorted.bam,nanopore2D_sorted.bam,nanopore1D_sorted.bam

f:id:kazumaxneo:20171006212517j:plain

上からshortのみ、hybrid、longのみ、そしてシミュレートしたナノポアのリードをマッピングしたものとなる。一番下は同様のやり方で発生させた1Dのデータをマッピングしたものである(ナノポアのリードは-x on2Dをつけてbwa memを走らせている)。一番切れている部位が少ないのはhybrid-assemblyで作ったcontigである。longのみだとたまにアセンブルされていない領域があるようである。カバレッジがまだ飽和していないのかもしれない。

 

 

ACTでも比較してみる(リンク)。long readのみでアセンブルされたcontigはhybridと傾向が似ているので除く。

bwast.py -a short_scaffolds.fasta hybrid_scaffolds.fasta

 

f:id:kazumaxneo:20171006211123j:plain

一番上が公開されているゲノムで、boxがchr1~20まで並んでいる。真ん中がhybridのcontigで、一番下がshortのみのcontigである。hybridのcontigはchr全体とアライメントされているものが多い。一方、shortのみのcontigは線画あちこちのcontigから出ている。まだまだバラバラの状態ということになる。

 

例えばchr11を少し拡大してみる。

f:id:kazumaxneo:20171006211808j:plain

 黄色の部分がchr11だが、一番上の段のゲノムと真ん中の段のhybridのcontigは1:1で対応しているのに対し、hybridと一番下のshortの比較では shortのみのcontigは線があちこちに散らばっている。これはshortのみのcontigはロングより小さなcontigにまだ分かれているということである。

 

chr13をIGVでも見てみる。

f:id:kazumaxneo:20171006212948j:plain

shortのみのcontigでは何度も切れている部位があることがわかる。

 

もう1つ他のchromosomeを見てみる。下はchr3。

f:id:kazumaxneo:20171006214326j:plain

切れている領域を拡大。ショートリードで全然繋がらない部位もロングリードがあれば繋がっている。

f:id:kazumaxneo:20171006214330j:plain

 

 

hybridアセンブルの結果をbandageを使いgraphパスで確認する(やり方)。

f:id:kazumaxneo:20171014172614j:plain

クロモソームレベルのアセンブリができているためか、長いscaffoldは単独のものが多い。

 

今回は面倒なのでやっていないが、jellyfishなどでk-merを指定して配列を書き出し、それをゲノムに当てれば繰り返し領域の部分がカバレッジの肥大部位として可視化できます(ついでにwigに変換してもO.K)。切れた部位とリピートの関係を調べるために、k-merのサイズを変えてやってみると面白いと思います(どこまでk-merサイズをあげれば繋がりうるか推定できる)。

 

 

 

 

 

 

次はゲノムサイズ41MのZymoseptoria_tritici.MG2。やはりNがない完全にFinishしたゲノムになっていたので選んだ。

f:id:kazumaxneo:20171006214914j:plain

Ensembl fungiより。クロロプラストゲノムとミトコンドリアゲノムは今回使わない。

 

シゾント同じ手順でQ30以上のショートリードx100と、ロングリードを6万リードだけ発生させた。

user$ seqkit stats *

file                   format  type    num_seqs        sum_len  min_len  avg_len  max_len

R1.fastq               FASTQ   DNA   19,025,739  3,044,118,240      160      160      160

R2.fastq               FASTQ   DNA   19,025,739  3,044,118,240      160      160      160

simulated_reads.fasta  FASTA   DNA       20,000    164,947,284       90  8,247.4   59,524

 

このケースではゲノムサイズが40Mbなので、ロングリードのカバレッジはx4にしかならない。そのためかMiniasmではアセンブルできなかった(極端に短いcontigが少しだけできた)。つまりこのケースでは、ショートリードだけでおそらく完全なアセンブルはできないし、ロングリードだけでもアセンブルできないということである。カバレッジとしてはごくわずかなロングリードを混ぜることで、アセンブルはどこまで改善するだろうか?

 

以下がアセンブル結果をblastnで比較しACTで描画したものである。 一番上が公開されているゲノム、真ん中がhybridのcontig、一番下がshortのみのcontigである。

f:id:kazumaxneo:20171008234921j:plain

シゾンのケースほどhybrid assemblyは完全ではない。chr1を拡大してみる。

f:id:kazumaxneo:20171008235231j:plain

上のオレンジ部分が公開されているゲノムのchr1だが、それと下のhybrid assemblyの位置を対応させる黄色の線はあちこちから出ていることがわかる。

 

IGVでも確認してみる。上がhybrid assemblyのscaffold、真ん中がshort readのみのassembly、一番下がnanoporeの2Dリードである。scaffoldが切れている中央部位に2Dリードがないのがわかる。だが一方で、short readのみで得られたcontigよりcontiguityは改善している。

f:id:kazumaxneo:20171009001205j:plain

 

拡大してみる。

f:id:kazumaxneo:20171009002556j:plain

左のほうの部位は、ショートリードで切れている部位にロングリードが存在しており、そのため、hybrid assemblyでつながったと推測される。一方、真ん中の部位は、ロングリードがシーケンスされておらず、そのためhybrid assemblyでも切れていると推測される。hybrid assemblyでも切れている他の部位も確認したが、 ロングリードが見つからなかった。2Dのraw long readで4xというスループットはchromosomeレベルのアセンブリには不十分かもしれない。

 

14xのロングリードでやり直した。quastで評価する(使い方)。

f:id:kazumaxneo:20171012012659j:plain

 左端からillumina only、真ん中がillumina+2D_read(x3)、右端がillumina+2D_read(x14)。

f:id:kazumaxneo:20171012012816j:plain

 

が illumina only、がillumina+2D_read(x3)、がillumina+2D_read(x14)。

やはり14xの方が良い結果が出ている。このhybridアセンブリの結果をbandageを使いgraphパスで確認する 。

f:id:kazumaxneo:20171014173610j:plain

シゾンの結果と似ているが、シゾンよりは線が絡まったlong scaffoldが多い印象を受ける。long readがまだ足りないのかもしれない。例えば中央左寄りの真っ黒になったcontigを拡大してみる。

f:id:kazumaxneo:20171014174010j:plain

さらに拡大。lengthを表示。大半が500bp以下の極端に短いcontigである。

f:id:kazumaxneo:20171014223002j:plain

coverageに切り替えると、短いcontigのカバレッジが左右の長いcontigの倍以上あるものが多いことがわかる。

f:id:kazumaxneo:20171014223112j:plain

 

どうやらこのような繰り返し配列がクラスターを作った非常に複雑な領域があるらしい。そこでアセンブルが止まっていると推測される。

 

 

genomeがコンプリートにどれだけ近づいているかは、登録されている遺伝子がどれだけ出てくるかで1つ評価できる。FTPサイトにcomplete cDNAがあったので、cDNAが部分的でもscaffoldからいくつ見つかるかで評価してみる。方法としては不完全だが1つの指標とはなる。

cDNAセット 10967個

illumina only (x100) => 10966 hit

illumina (x100) + 60000 long read (3x) => 10966 hit

illumina (x100) + 200000 long read (3x) => 10966 hit

 

cDNAが部分的にヒットすればカウントしているので(blastn使用 <1-e-10)、かなり誤差はあるだろうが、illuminaのリードだけでも遺伝子ははほぼフルセットあるのかもしれない。つまり、目的がゲノムを決めることでなくwhole genome re-sequencingによる変異部位特定にあるなら、illuminaのデータセットだけでも勝負できる可能性が高いことを意味している。ただしこれはゲノムの複雑さによって変わってくるもので、ゲノムサイズが100-Mbほどの程よいサイズとしても、ヒトのようにリピートが非常に多くて、変異の多くがcopy umber variationとして効いてくる生き物ではショートリードのアセンブルだけでは難しくなると考えられる。ただし、それがどのくらい影響してくるのか判断するのはとても難しい。

 

 

 

 

遺伝子情報のない新規ゲノムならBUSCOなどの評価手法を使う。BUSCOは以前紹介しています。


 

 

 

補足

今回、1Dのデータではエラーが多くてminiasmやcanuでアセンブルはできないし、bwa mem -x 0n2Dでもほとんどアライメントできなかった。これはつまり、汚すぎてエラーコレクションかpolishでなんとかしないと使えないリードになっていることを意味する。

そのあたりの手法については徐々に論文で報告され初めてきています。今後紹介していきます。