ハイクオリティなショートリードのデータに、ロングリード情報を混ぜ込むとどれくらいアセンブリは改善されるのか調べてみる。
NがあってもgrepやUCSC ゲノムブラウザで除くことができるが、それでは肝心の繰り返し領域の評価が曖昧になる。やはりNがないコンプリートな真核ゲノムを使いたいということで、100M以下の生き物で調べたところ、シゾンが完全にFinishしていた。そこでまずはシゾンのゲノムを例にアセンブルのパフォーマンスを調べてみる。Ensemblからfastaファイルをダウンロードする。
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。
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
上から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
一番上が公開されているゲノムで、boxがchr1~20まで並んでいる。真ん中がhybridのcontigで、一番下がshortのみのcontigである。hybridのcontigはchr全体とアライメントされているものが多い。一方、shortのみのcontigは線画あちこちのcontigから出ている。まだまだバラバラの状態ということになる。
例えばchr11を少し拡大してみる。
黄色の部分がchr11だが、一番上の段のゲノムと真ん中の段のhybridのcontigは1:1で対応しているのに対し、hybridと一番下のshortの比較では shortのみのcontigは線があちこちに散らばっている。これはshortのみのcontigはロングより小さなcontigにまだ分かれているということである。
chr13をIGVでも見てみる。
shortのみのcontigでは何度も切れている部位があることがわかる。
もう1つ他のchromosomeを見てみる。下はchr3。
切れている領域を拡大。ショートリードで全然繋がらない部位もロングリードがあれば繋がっている。
hybridアセンブルの結果をbandageを使いgraphパスで確認する(やり方)。
クロモソームレベルのアセンブリができているためか、長いscaffoldは単独のものが多い。
今回は面倒なのでやっていないが、jellyfishなどでk-merを指定して配列を書き出し、それをゲノムに当てれば繰り返し領域の部分がカバレッジの肥大部位として可視化できます(ついでにwigに変換してもO.K)。切れた部位とリピートの関係を調べるために、k-merのサイズを変えてやってみると面白いと思います(どこまでk-merサイズをあげれば繋がりうるか推定できる)。
次はゲノムサイズ41MのZymoseptoria_tritici.MG2。やはりNがない完全にFinishしたゲノムになっていたので選んだ。
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である。
シゾンのケースほどhybrid assemblyは完全ではない。chr1を拡大してみる。
上のオレンジ部分が公開されているゲノムのchr1だが、それと下のhybrid assemblyの位置を対応させる黄色の線はあちこちから出ていることがわかる。
IGVでも確認してみる。上がhybrid assemblyのscaffold、真ん中がshort readのみのassembly、一番下がnanoporeの2Dリードである。scaffoldが切れている中央部位に2Dリードがないのがわかる。だが一方で、short readのみで得られたcontigよりcontiguityは改善している。
拡大してみる。
左のほうの部位は、ショートリードで切れている部位にロングリードが存在しており、そのため、hybrid assemblyでつながったと推測される。一方、真ん中の部位は、ロングリードがシーケンスされておらず、そのためhybrid assemblyでも切れていると推測される。hybrid assemblyでも切れている他の部位も確認したが、 ロングリードが見つからなかった。2Dのraw long readで4xというスループットはchromosomeレベルのアセンブリには不十分かもしれない。
14xのロングリードでやり直した。quastで評価する(使い方)。
左端からillumina only、真ん中がillumina+2D_read(x3)、右端がillumina+2D_read(x14)。
赤が illumina only、青がillumina+2D_read(x3)、緑がillumina+2D_read(x14)。
やはり14xの方が良い結果が出ている。このhybridアセンブリの結果をbandageを使いgraphパスで確認する 。
シゾンの結果と似ているが、シゾンよりは線が絡まったlong scaffoldが多い印象を受ける。long readがまだ足りないのかもしれない。例えば中央左寄りの真っ黒になったcontigを拡大してみる。
さらに拡大。lengthを表示。大半が500bp以下の極端に短いcontigである。
coverageに切り替えると、短いcontigのカバレッジが左右の長いcontigの倍以上あるものが多いことがわかる。
どうやらこのような繰り返し配列がクラスターを作った非常に複雑な領域があるらしい。そこでアセンブルが止まっていると推測される。
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でなんとかしないと使えないリードになっていることを意味する。
そのあたりの手法については徐々に論文で報告され初めてきています。今後紹介していきます。