macでインフォマティクス

macでインフォマティクス

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

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を開く。

f:id:kazumaxneo:20170614112124j:plain

凡例を拡大表示。

f:id:kazumaxneo:20170614134504j:plain

 

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

 

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

 -z: ズーム倍率の指定

f:id:kazumaxneo:20170614134122j:plain

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

 

 

 

チュートリアル7へ

 

 

 

#参考

比較する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の1、2行目だけ残し、>の名前も特殊文字を含まないものに変える必要がある。また、拡張子は.faは受け付けず、.fnaにする必要がある。一方、特に改行は必要ない。

>read_1

ATGC....

>read_2

GTGC....

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

 

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

 

備考

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

 

追記

seqkitも使えます。