macでインフォマティクス

macでインフォマティクス

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

IGVのtips 分割表示やGC変化の表示

 

banファイルを読み込んでリード情報を可視化できるフリーソフト。多様なオプションが用意されており、使いこなすにはかなり勉強が必要。   例を挙げて説明する。 

 

1、Split view 表示

下のようにウィンドウを分けて表示できる。本来はmate-pairの相方リード表示用だが、複数領域を同時に比較することにも利用できる。   やり方はMiなどのテキストエディタにあらかじめタブで区切った二つのポジションを書いておき、それをウィンドウ内にペーストするだけでOK。

例えばchrの471900-474900と519,700-522,701を同時にみたいなら、間をタブで区切り以下のように記載する。

chr:471900-474900 chr:519,700-522,701  

二つの領域が表示される。以下は4つ表示した例。f:id:kazumaxneo:20170324173457j:plain

split viewを解除するには、上のウィンドウの空間でダブルクリックする。

 

 

2、変異のコール結果のVCFファイルをIGVに取り込ませる

 

 gatkの検出結果を読み込ませてみる。例えば

 

gatk -T UnifiedGenotyper -R $f -I R1R2_sorted.bam -o UnifiedGenotyper.vcf

 

出力されるvcfファイルを読み込ませる。

f:id:kazumaxneo:20170324173521j:plain

indel検出ソフトの出力はVCF形式に対応しているものが多い。変異をコールして気になる部位があればIGVで見てみればいい。主観なので危険性は伴うが、偽陽性かどうか推測できる場合も多い。

 

 

 

 

3、ゲノムのGC含量をIGVで表示。

例えば100bpのウィンドウサイズで描画したいとするなら、以下のように行う(1-bpのサイズは意味がないので、ある幅のウィンドウサイズで計算する)

 

 

解析に必要なもの

 

・samtools 

・bed tools

gawk (macOSXで計算する時) 

 

macOSXへのgawkの導入は

brew install gawk

(homebrewがない人は前もって導入しておく。ターミナルでbrewと打ってコマンドが見つからなければ入ってない)

gawkと同様にsamtoolsとbedtoolsの導入は

brew install samtools

brew install bedtools

 


ゲノム(test.faとする)のGC%を100bpのウィンドウサイズで計算してIGVで可視化するには、以下のようなフローで行なう。

 

1、fastaのindex作成

samtools test.fa

 

2、カラム1、2を抽出

cut -f 1,2 test.fasta.fai > test.sizes

 

3、GC計算幅を指定するファイルを作成

bedtools makewindows -g test.sizes -w 100 > test_100bps.bed

 

(w=100ならこんなのができる)

chr 0 100

chr 100 200

chr 200 300

chr 300 400

chr 400 500

….

 

4、ATGC分布を計算。

bedtools nuc -fi test.fasta -bed test_100bps.bed > test_nuc.txt

 

こんなのができる

head -5 test_100bps.bed

#1_usercol 2_usercol 3_usercol 4_pct_at 5_pct_gc 6_num_A 7_num_C 8_num_G 9_num_T 10_num_N 11_num_oth 12_seq_len

chr 0 100 0.390000 0.610000 20 25 36 19 0 0 100

chr 100 200 0.560000 0.440000 20 18 26 36 0 0 100

chr 200 300 0.580000 0.420000 31 15 27 27 0 0 100

chr 300 400 0.500000 0.500000 23 23 27 27 0 0 100

 

 

 

5、必要なカラムを編集しながら抽出。gawkを使っている。

gawk -v w=100 'BEGIN{FS=" "; OFS=" "}

{if (FNR>1) {print $1,$2,$3,"GCpc_"w"bps",$5}}' test_nuc.txt > test_100bp.igv

 

 

成功するとこんなファイルができる。

head -4 test_100bp.igv

chr     0       100     GCpc_100bps     0.610000

chr     100     200     GCpc_100bps     0.440000

chr     200     300     GCpc_100bps     0.420000

chr     300     400     GCpc_100bps     0.500000

 

 

6、test_100bp.igvをIGVに読み込ませる。

 

 

 

 

 

 

 

おまけ;
IGVtoolsでバイナリー化しておけば軽くなる。

 

IGVtools -> Run igvtools

f:id:kazumaxneo:20170324173727j:plain

Command -> toTDF を選択

input file Browse -> .igvファイルを選択。

Runをクリック

 

f:id:kazumaxneo:20170324173641j:plain

test_100bp.igvができる。Load from genomesから読み込ませる。igvtoolsは、この例のようにIGVのGUIから呼び出して計算させることができる。

 

 

 

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

参考ページ

 

http://wiki.bits.vib.be/index.php/Create_a_GC_content_track