次は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)。
チュートリアル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も使えます。