CGView Comparison Toolによるゲノム比較 実践編 -大量のバクテリアゲノムの同時比較

インストールは以下で説明しています。

チュートリアルの総仕上げとして、CCTのコマンドfetch_all_refseq_bacterial_genomes.shを使って、登録されているバクテリアのrefseq配列全てを自動ダウンロードして、リファレンスゲノムと比較してみることにする。リファレンスとして、まだゲノムが登録されていない未知バクテリアのドラフトゲノムを用意した。ドラフトゲノムはRASTでアノテーションをかけてgbkファイル化してある。これを1148.gbkとする。

プロジェクト構築までの流れは以前と同じである。

1148.gbkのプロジェクトをビルドする。

build_blast_atlas.sh -i 1148.gbk

1148ディレクトリができる。

fetch_all_refseq_bacterial_genomes.shコマンドを使い、登録されているバクテリアゲノムをダウンロード。

fetch_all_refseq_bacterial_genomes.sh --min 1000000 -o 1148/comparison_genomes/

プラスミドゲノム除去のため、--min 1000000をつけて実行したが、それでも合計8272ファイル見つかった（--minなしだと16680）。ただし検索中にhitしなかった配列がかなりの割合で出現し、最終的には3142gbkファイルダンロードされた。

 追記 2018 05実行時は11385ゲノムダウンロードされた。

比較する。

build_blast_atlas.sh -p 1148 --memory 10000m --max_blast_comparisons 3200 
 --map_size x-large --custom 'width=40000 height=40000 backboneRadius=16000 
 featureThickness=240 rulerFontSize=200 rulerPadding=400  
 tickThickness=30 tickLength=80 draw_divider_rings=F  
 _cct_blast_thickness=4.0 labelFontSize=100'

CGView Comparison Toolによるバクテリアのゲノム比較7 -大量のミトコンドリアゲノムの同時比較

インストールは以下で説明しています。

公式ページのチュートリアル7を実践していく。

人のミトコンドリアゲノムを、他の生物のミトコンドリアゲノムと比較する。

ミトコンドリアゲノムNC_012920をダウンロードする。

fetch_genome_by_accession.sh -a NC_012920 -o ./

NC_012920のプロジェクトをビルドする。

build_blast_atlas.sh -i NC_012920.gbk

NC_012920ディレクトリができる。

他の生物のミトコンドリアゲノム全てをダウンロードする。それには専用のfetch_all_refseq_mitochondrial_genomes.shコマンドを使う。

fetch_all_refseq_mitochondrial_genomes.sh -o NC_012920/comparison_genomes/

comparison_genomes/に740ファイルダウンロードされた（５時間かかった）。

比較を行う。

build_blast_atlas.sh -p NC_012920 --memory 2500m --max_blast_comparisons 2500 
 --map_size x-large --custom 'width=20000 height=20000 backboneRadius=8000 
 featureThickness=120 rulerFontSize=100 rulerPadding=200 
 tickThickness=15 tickLength=40 draw_divider_rings=F 
 _cct_blast_thickness=2.0 labelFontSize=200'

巨大なファイルを扱うため、いくつかのオプションを使っている。

--memory 2500m

--map_size x-large: mapをx-largeのみ出力（小さい図は情報が潰れて見えない）。

--max_blast_comparisons 2500: 標準ではtop100ゲノムのみプリントされる。100以上の結果を出力するために使う。

--customコマンドでrulerFontSize、rulerPadding、tickThickness、tickLength、draw_divider_rings、_cct_blast_thickness、labelFontSizeを指定している。

• rulerFontSize; Size of font used for position numbers
• rulerPadding; Space between innermost circle and position numbers
• tickThickness; Thickness of tickmarks
• tickLength; Length of tickmarks
• draw_divider_rings; Draw thin rings between BLAST rings and other rings
• tickLength; Length of tickmarks
• _cct_blast_thickness; Thickness of BLAST rings when '-cct' option used
•  labelFontSize; Size of font used for feature labels

出力は以下のようになった。 

DNA vs DNA

参考までに、customのオプションをサイズ以外指定しないで描画したのが下の図になる。

build_blast_atlas.sh -p NC_012920 --memory 2500m --max_blast_comparisons 800 
 --map_size x-large --custom 'width=20000 height=20000 backboneRadius=8000 '

DNA vs DNA 

1-7まで紹介したが、公式ページのチュートリアルでは、特定orfに注目して解析を行う方法や、高度な検索を行うコマンドも紹介しています。確認してみてください。

次に、手に入るバクテリアゲノム全部を使った比較を実行してみる。

CGView Comparison Toolによるバクテリアのゲノム比較6 -次世代リードのアライメント

次はCCTを使って次世代データをリファレンスに当てて、リードの張り付きをビジュアル化するチュートリアルを見ていく（公式ページチュートリアル6）。

CCTのインストールは以下で説明しています。

 今までと流れが微妙に異なっているので注意する。

francisellaプロジェクトをビルドする。

cgview_comparison_tool.pl -p francisella

francisellaディレクトリができる。

Francisella tularensis (野兎病菌) のゲノムをダウンロードする。

fetch_genome_by_accession.sh -a NC_006570 -o francisella/reference_genome

NC_006570.gbkがダウンロードされる。

454シーケンサーのファイルをダウンロードする。#ペアリードの1

wget http://stothard.afns.ualberta.ca/cgview_server/examples/example_2/files/454_reads_1.fast -O francisella/comparison_genomes/454_reads_1.fna 

#ペアリードの2

wget http://stothard.afns.ualberta.ca/cgview_server/examples/example_2/files/454_reads_2.fasta -O francisella/comparison_genomes/454_reads_2.fna 

中身を見てみる。

>head -3 francisella/comparison_genomes/454_reads_1.fna
>read_1
taaatttttgattttttttactttttatatgcttgatatatctcatataatccttgcgcc
tcaagtattgatatattaatatcaggatctacttcagca

NGSに対応しているが、リードはfastaフォーマットになってないとダメらしい（リードのフォーマットについては一番下に書いています）。

解析するにはコンフィグファイル(francisella/project_settings.conf)の設定を修正する必要がある。

query_source = none

のnoneをnucleotideに修正。

database_source = none

のnoneをdnaに修正。

できたらランする。

cgview_comparison_tool.pl -p francisella -t --cct

 --cct: つけるとリードのアライメント相同性により色が変わる。

-t; リードを整頓して相同性が一番高いリードを一番外側に並べ換える。

mapのmedium.pngを開く。

凡例を拡大表示。

描画結果の赤い部分は454のリードが高密度にアライメントされた領域である。ただし稀に白くなっている領域がある。左下の白の部分を25倍拡大してみる。その場合は、チュートリアル2で出てきたcreate_zoomed_maps.sh コマンドを使う。

create_zoomed_maps.sh -p francisella -c 1080000 -z 25

 -z: ズーム倍率の指定

orfの注釈が小さければxmlを編集して拡大する（チュートリアル4）。 

チュートリアル７へ

#参考

比較するNGSのデータは、fasta形式になっていないといけない。またfastaの名前に特殊文字が入っているとblast実行時にエラーになってしまう。すなわち

>1
CCCTATGCCGCCAGCAAACGCATGGTGGAACAAATTTTGGCTGACTTTGACCAAGCTTACGTCTTCCAATCAGTTATTTTTCGCTACTTTAATGCTTCCGGCGCACACCCCCAAGGACTCTTAGGAGAAGACCATAACCCCGAAACCCATCTAATTCCCCTTCCCCTGTTTACAGCTTTGAAACAACGTCCCCAACTTTCCTTTTTTGGTACTGATTATGACACCCTCGACGTCACCGCCTTACGGGATTACATCCATTTCTGTGCCCTCGCGATCGCCCATGTATTAGTCTTGCAATACC
>2
ACAGGGGACTGTTGAATCAGAATTTGATGGGTTTGCCAAAAACTGCGGTATTTGGGGCCGTAGCCCAAAAATTGCGTCGTTTGGGGCAGGGTTTTTCTGACCGCCTTTCTGTAGAGGCTTTCACGAATATATTTTCCTTCCCCGTCTATATAGCCCGGACCCTGATCCGCCACGGCCATCGCCACCGTTTTTCTTTCCTGTTCTTCCACATTTACCTTCGCCGGAAATAACTTTCTAGTTACTTCTACTCTCATTCTGGCCGGATTTTTTGCCAGTTCATCCCCCAAGGCCTGTGTAGCA 

は比較できるが、

>M02077:18:000000000-A88N7:1:1101:20381:1495 1:N:0:1 
CCCTATGCCGCCAGCAAACGCATGGTGGAACAAATTTTGGCTGACTTTGACCAAGCTTACGTCTTCCAATCAGTTATTTTTCGCTACTTTAATGCTTCCGGCGCACACCCCCAAGGACTCTTAGGAGAAGACCATAACCCCGAAACCCATCTAATTCCCCTTCCCCTGTTTACAGCTTTGAAACAACGTCCCCAACTTTCCTTTTTTGGTACTGATTATGACACCCTCGACGTCACCGCCTTACGGGATTACATCCATTTCTGTGCCCTCGCGATCGCCCATGTATTAGTCTTGCAATACC
>M02077:18:000000000-A88N7:1:1101:12126:1509 1:N:0:1
ACAGGGGACTGTTGAATCAGAATTTGATGGGTTTGCCAAAAACTGCGGTATTTGGGGCCGTAGCCCAAAAATTGCGTCGTTTGGGGCAGGGTTTTTCTGACCGCCTTTCTGTAGAGGCTTTCACGAATATATTTTCCTTCCCCGTCTATATAGCCCGGACCCTGATCCGCCACGGCCATCGCCACCGTTTTTCTTTCCTGTTCTTCCACATTTACCTTCGCCGGAAATAACTTTCTAGTTACTTCTACTCTCATTCTGGCCGGATTTTTTGCCAGTTCATCCCCCAAGGCCTGTGTAGCA 

は比較できない。fastqを比較したいなら、スクリプトを書いて、fastqの１、２行目だけ残し、>の名前も特殊文字を含まないものに変える必要がある。また、拡張子は.faは受け付けず、.fnaにする必要がある。一方、特に改行は必要ない。

>read_1

ATGC....

>read_2

GTGC....

このように、１行目ヘッダー、２行目長い配列でも受け付けてくれる（foldコマンドで折りたたむ必要は特にない）。

NGSのアライメントをblastで行うとやはり相当な時間がかかってしまう。NGSデータを使いたいなら、データをランダムサンプリングしてカバレッジを減らしておいた方が良い。

備考

ランダムサンプリングならこの方の記事が参考になると思います。

追記

seqkitも使えます。 

CGView Comparison Toolによるバクテリアのゲノム比較5 - ゲノム比較とモニタージュ合成

葉緑体、ミトコンドリアの次は、CCTを使って複数ゲノムを比較するチュートリアルを見ていく（公式ページチュートリアル5）。

CCTのインストールは以下で説明しています。

全体比較は手順が異なる。まずbuild_blast_atlas_all_vs_all.shコマンドを使い、新しいプロジェクトをビルドする。ここではmontage_projectを作る。

build_blast_atlas_all_vs_all.sh -p montage_project

montage_projectディレクトリができる。

Bordetellaクロモソームをcomparison_genomes/にダウンロードする。

fetch_refseq_bacterial_genomes_by_name.sh -n "Bordetella*" --min 1000000 -o montage_project/comparison_genomes/

--min 1000000でゲノムサイズ1Mbp以下を除く。

321ファイルダウンロードされた。

全体比較には再びbuild_blast_atlas_all_vs_all.shを打つ。

build_blast_atlas_all_vs_all.sh -p montage_project

このランは非常に時間がかかる。

エラーになり、最後のモニタージュのデータが出力されるまでランできなかった。修正でき次第投稿します。このページは情報共有のため残しておきます。

チュートリアル６へ

CGView Comparison Toolによるゲノム比較4 - ミトコンドリアゲノムの比較

葉緑体ゲノムに続き、CCTを使ってミトコンドリアゲノムを比較するチュートリアルを見ていく（公式ページチュートリアル4）。

CCTのインストールは以下で説明しています。

前半は以前のクロモソーム、プラスミドと同じなので簡潔に説明する。

まずはドブネズミのミトコンドリアゲノムをダウンロードする。

fetch_genome_by_accession.sh -a AC_000022 -o ./

上記のコマンドを打つと、カレントディレクトリにAC_000022がダウンロードされる。

blastの準備をする。以下のコマンドを打つ。

build_blast_atlas.sh -i AC_000022.gbk

終わるとAC_000022ディレクトリができる。

他のミトコンドリアゲノムをダウンロードする。

fetch_all_refseq_mitochondrial_genomes.sh -o AC_000022/comparison_genomes/

740ファイルダウンロードされた。

デフォルト条件でゲノムを比較するなら以下のように打つ。

build_blast_atlas.sh -p AC_000022

DNA vs DNA 

 top100が表示されている。

cDNA vs cDNA

  top100が表示されている。

しかしこのままだとorfの注釈が非常に小さい。

小さすぎる！。 

そこで公式チュートリアルには、見栄えを良くするためxmlファイルのパラメータを変更する方法を挙げている。それには、cDNAのフォルダかDNAのフォルダを開き、xmlを編集する方法がある。例えば、cDNAのフィルダ(AC_000022/maps_for_dna_vs_dna/cgview_xml/)にはx-largeとlarge、mediumそれぞれの図のxmlファイルがある。例えばx-largeの文字を変えるなら、下のサイズに変えると随分大きくなる。

x-largeのxmlファイルをMiなどのテキストエディタで開き、下の４つを検索する。

labelFont="SansSerif, plain, 15"

のサイズを15から200に変更。

rulerFont="SansSerif, plain, 60"

のサイズを60から100に変更。

tickLength="20"

を20から40に変更。

tickDensity="0.0416666666666667"

を0.0416666666666667から0.08んい変更。

このサイズはx-largeでは程よいが、largeやmediumでは大きくなりすぎてしまう（largeなら上の数値の半分以下、mediumで1/4以下）。

AC_000022のルート直下に戻りランを再実行する。ただし前のコマンドを打つとxmlが上書きされてしまう。改変したxmlを使って実行するなら以下のように打つ。

redraw_maps.sh -p AC_000022

修正後のDNA vs DNAのpngを開いてみる。

orfの注釈が見やすく大きくなっている。  

ベクターグラフィックスのsvg形式で出力してみる。

redraw_maps.sh -p AC_000022 --format svg

 DNA vs DNA

ベクターグラフィックなので、いくら拡大してもドット化しない。

svg形式はadobe系のソフトと必ずしも互換性があるわけではないが、CCTで出力したsvgファイルはイラストレーターで開くことができることを確認した。

今回は公式チュートリアルに習いxmlファイルをいじったが、--customでフォントサイズを変えても同じことはできる。文字を今回の修正版のように大きくするならは以下のようなコマンドになる。

build_blast_atlas.sh -p AC_000022 --custom 'labelFontSize=100 backboneRadius=3000'

出力

チュートリアル５へ

CGView Comparison Toolによるゲノム比較3 - 葉緑体ゲノムの比較

プラスミド、クロモソームに続き、CCTを使って葉緑体ゲノムを比較するチュートリアルを見ていく（公式ページチュートリアル3）。

CCTのインストールは以下で説明しています。

前半は以前のクロモソーム、プラスミドと同じなので簡潔に説明する。

 まずはPorphyra purpureaのクロロプラストゲノムをダウンロードする。

fetch_genome_by_accession.sh -a NC_000925 -o ./

上記のコマンドを打つと、カレントディレクトリにNC_000925.gbkがダウンロードされる。 

プロジェクトをビルドする。

build_blast_atlas.sh -i NC_000925.gbk

終わるとNC_000925ディレクトリができる。

他のクロロプラストゲノムをダウンロードする。

fetch_all_refseq_chloroplast_genomes.sh -o NC_000925/comparison_genomes/

157ファイルダウンロードされた。

デフォルト条件でゲノムを比較するなら以下のように打つ。

build_blast_atlas.sh -p NC_000925 --map_size x-large --custom 'labelFontSize=100 backboneRadius=3000'

 --custtomで文字が潰れないよう工夫している。

• labelFontSizeはfeature tableの文字サイズ
• backboneRadiusは背景の円のサイズ。

--map_size x-large: x-largeのファイルのみ出力。

DNA vs DNA

cDNA vs cDNA

チュートリアル４へ

CGView Comparison Toolによるバクテリアのゲノム比較2 - クロモソームの比較

インストールは以下で説明しています。

プラスミドに続き、E.coliのゲノムを他のE.coliゲノムと比較してみる（公式ページのチュートリアル2）。

ゲノムをダウンロード。

fetch_genome_by_accession.sh -a CP001855 -o ./

CP001855.gbkがダウンロードされる。CP001855というgenbankのアクセッション番号で指定している。アクセッション番号はgenbankファイルのtop付近にある。

NCBIのNucleotideデータベース。青で強調しているのがアクセッション番号。

プロジェクトをビルドする。

build_blast_atlas.sh -i CP001855.gbk

CP001855ディレクトリができ、その中に複数のサブディレクトリとファイルができる。

他のE.coliゲノムをダウンロードする。

fetch_refseq_bacterial_genomes_by_name.sh -n "Escherichia*" --min 1000000 -o CP001855/comparison_genomes/

--min: ゲノムサイズ1Mbp以下は除く。（E.coliのプラスミドを排除するため）

343ファイルダウンロードされた。

デフォルト条件でゲノムを比較するなら以下のように打つ。

build_blast_atlas.sh -p CP001855

ゲノムになると、解析にかなりの時間がかかる。焦らず待つ。

一晩放置していると終わった。以下のフォルダができる。

CP001855/maps_for_dna_vs_dna

CP001855/maps_for_cds_vs_cds

 maps_for_dna_vs_dna/DNA vs DNAのpngを開いてみる。 

一部の領域を拡大して見るなら、create_zoomed_maps.shを使う。

create_zoomed_maps.sh -p CP001855 -c 4450000 -z 15

-z: ズームサイズ。

-c: 大まかな領域の指定。

 maps_for_dna_vs_dna/DNA vs DNAのpngを開いてみる。 

中央の領域が他のE.coliゲノムにはないことがわかる。

さらにズームしてみる。

create_zoomed_maps.sh -p CP001855 -c 4450000 -z 60

DNA vs DNAのpngを開くと

このようになった。全体マップを描いて興味がある領域が出てくれば、create_zoomed_maps.shコマンドを使うことで欲しい領域だけ取ってくることができる。

 create_zoomed_maps.shコマンドは--customオプションが存在しないので、フォントサイズをいじる場合はxmlファイルを直接編集する（xmlを編集する例はチュートリアル4参照)か、ゲノム全体を描画する時に前もって

--custom 'labelFontSize=100'

をつけてランする。"20"だとゲノム全体を描画すると見えないが、特定の領域だけ拡大するzoomコマンドを走らせると、x-large、large、mediumの３出力のうちmediumの出力でちょうど良いフォントサイズになる。例えば下の図は、20のフォントサイズで全体を描画してから、create_zoomed_maps.shでx15拡大して再計算しmediumの図を開いたものになる。

ちょうど良いフォントサイズ。

x-largeだとこのようにフォントは見えない。

つまりx-large出力を使うなら、20の数倍の値にしておく必要がある。

チュートリアル３へ