macでインフォマティクス

読者です 読者をやめる 読者になる 読者になる

IGVの使い方

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

 

例を挙げて説明する。

 

 

1、Split view 表示

下のようにウィンドウを分けて表示できる。本来はmate-pairの相方リード表示用だが、複数領域を同時に比較することにも利用できる。

 

やり方はMiなどのテキストエディタにあらかじめタブで区切った二つのポジションを書いておき、それをウィンドウ内にペーストするだけ。例えばchrの471900-474900と519,700-522,701を同時にみたいなら以下のように記載する。(中央はタブ)

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

 

3以上ももちろんいける。以下は4つ表示した例。

f:id:kazumaxneo:20170324173457j:plain

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

 

 

2、gatkやbed toolsのコール部位をIGVで表示

gatkの結果を読み込ませるだけ。

f='/Users/user/test.fa'

p='R1.fastq'

q='R2.fastq'

 

bwa index -p $f

bwa mem -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" -t 12 $f $p $q > R1R2.sam

samtools view -@ 20 -bS R1R2.sam > R1R2.bam

samtools sort -@ 20 R1R2.bam >  R1R2_sorted.bam

samtools index R1R2_sorted.bam

samtools faidx $f

 

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

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

 

 

できた.vcfファイルを読み込ませるだけ。上の文全体をperl形式やシェルスクリプト形式で保存しておけば、ワンライナーでgatkのコールまでできて便利。gatkのlow quality callはあらかじめフィルタリングしておくとよい。

f:id:kazumaxneo:20170324173521j:plain

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

 

 

 

 


2、gatkやbed toolsのコール部位をIGVで表示

gatkの結果を読み込ませるだけ。

f='/Users/user/test.fa'

p='R1.fastq'

q='R2.fastq'

 

bwa index -p $f

bwa mem -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" -t 12 $f $p $q > R1R2.sam

samtools view -@ 20 -bS R1R2.sam > R1R2.bam

samtools sort -@ 20 R1R2.bam >  R1R2_sorted.bam

samtools index R1R2_sorted.bam

samtools faidx $f

 

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

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

 

 

できた.vcfファイルを読み込ませるだけ。上の文全体をperl形式やシェルスクリプト形式で保存しておけば、ワンライナーでgatkのコールまでできて便利。gatkのlow quality callはあらかじめフィルタリングしておくとよい。

f:id:kazumaxneo:20170324173546j:plain

図 gatkコールの表示(真ん中の赤と青の縦線がgatkコール)

 

 

 

 

 

 

 

 

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 test100bps.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="\t"; OFS="\t"}

{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に読み込ませる。

 

 

 

 

 

 

 

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

 

IGV tools -> 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から読み込ませる。バクテリアならテキストでも動くだろうが。

 

 

以上、ややマニアックな3つを紹介した。

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

参考ページ

 

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