bedGraph wiggleファイルの特徴領域の抜粋

#bedGraph wiggleファイル
リファレンスゲノム上にマッピングされた塩基対単位のリードスコアデータ
fixedStepとかvariableStepとかいろいろあるみたいだけどここでは割愛。詳しくはここらへん参照。
筆者はだいたい「.bigwig」からbigWigToWigで染色体ごとに分けて使ってる。

devide_each_chr.sh

#ヒト染色体の場合
for i in `seq 22`
do
bigWigToWig -chrom=chr$i hoge.bigwig hoge_chr$i.wig
done

shell楽
それでwiggleファイルの中身がこれ

hoge.wig

#bedGraph section chr1:0-840194
chr1	0	56898	0
chr1	56898	56998	1.02889
chr1	56998	91465	0
chr1	91465	91468	1.02889
chr1	91468	91469	2.05778
chr1	91469	91476	3.08639
chr1	91476	91487	4.11528
chr1	91487	91565	1.32123
chr1	91565	91568	0.9909

一行目でコメントアウトされてるのがデータの情報。
第1染色体の0塩基対目から840194塩基対目まで。
このコメントアウトが随所に見られるけど割と次のコメントアウトと前のコメントアウトで書かれてる範囲が繋がってたりする。IGVで見るためかも。
二行目から先は肝心のデータ部分

染色体	スタート位置	エンド位置	スコア

になってる。
スタート位置からエンド位置までスコアがずっと並んでるわけですね。

#何がしたいの
特定の閾値より高いor低いのが一定長以上続くような場所（特徴領域）を抜き出したい。
早い安いうまい古いで評判のPerlを使った。
抜き出した領域もbed形式風にまとめる。

search_region_on_wig.pl

use Statistics::Lite qw(:all);

my ($INPUT ,$OUTPUT, $border_u,$border_t, $length_u, $length_t, ) = ("","",0,0,0,3000000000);

($INPUT ,$OUTPUT, $border_u,$border_t, $length_u, $length_t, ) = @ARGV;

open IN,"$INPUT"or die;
open OUT,"> $OUTPUT"or die;
my $start = -1;
my $end = -1;
my @score = ();
my $chr;
while(my $line = <IN>){
    chomp $line;
    
    if($line !~ /#/){
        my @data = split /\s/,$line;
        if(($border_u<=$data[3])&&($data[3]<=$border_t)){#ボーダー範囲内に入っている
            if($start==-1){#スタート位置を覚えていなければ現在位置の記憶と終わりの記憶
                $start = $data[1];
                $end = $data[2];
                $chr = $data[0];
                @score = ();
                push @score,$data[3]*($data[2]-$data[1]+1);
            }
            else{
                if($end >= $data[1]){#一つ前のデータの終わりから繋がっているか
                    $end = $data[2];
                    push @score,$data[3]*($data[2]-$data[1]+1);
                }
                else{
                    if(($length_t>=$end-$start)&&($end-$start>=$length_u)){
                        print OUT "$chr $start $end ",$end-$start," ",sum(@score)/($end-$start),"\n";
                    }
                    #現在位置と終わりの記憶を忘れる
                    $start = -1;
                    $end = -1;
                }
            }
        }
        else{
            if($start != -1){#スタート位置に記憶が残っていたら、前のデータとの繋がりが切れると書き出す
                if(($length_t>=$end-$start)&&($end-$start>=$length_u)){
                    print OUT "$chr $start $end ",$end-$start," ",sum(@score)/($end-$start),"\n";
                }
                $start = -1;
            }
        }
    }
}
close OUT;
close IN;

これで

pickup.sh

for i in `seq 22`
do
perl search_region_on_wig.pl hoge_chr$i.wig output.txt 0 1.5 100 150
done

とすれば全染色体で探索ができる。

#まとめ
リードスコアが低い領域あつめたーいというためだけに書いたけどここから特に可視化ができるわけではない
IGVにマーカー機能があれば少し役に立つかも。