genome graphに遺伝子のスタート座標を乗せよう。

bedgraphでカスタムトラックに追加してから、

インポートでgenome graphへ。


##----------------------------------------##
#
# note batch file for gb2bedgraph.pl
#
##----------------------------------------##


usage :

perl all_chrgb2bed.pl xxxx.gb outdir


for ex.
[hoge@boss5 rn4]# perl all_chr.gb2bed.pl rat_uniq_gene.gb rn4start_bedgraph
[hoge@boss5 rn4start_bedgraph]# ls
..........
[hoge@boss5 rn4start_bedgraph]# head chr1.start_gene.bedgraph
chr1 1306646 1306647 1
chr1 2394947 2394948 1
chr1 5975207 5975208 1
chr1 6051558 6051559 1
...............
[hoge@boss5 rn4start_bedgraph]# cat chr*.start_gene.bedgraph > all_chr.start_gene.bedgraph

genomegraph uploaded as a custom track in the UCSC genome browser
coverage select









#!/usr/bin/perl
#
# bacth file for bed and bedgraph

$| = 1;

@chrarr = ('1', '2', '3', '4', '5',
'6', '7', '8', '9', '10',
'11', '12', '13', '14', '15',
'16', '17', '18', '19', '20',
'X');

if ($#ARGV < 0) {
print STDERR "\n";
print STDERR "Usage: all_chr.gb2bed(graph).pl xxx.gb outdir\n\n";
exit;
} else {
$gbfile = $ARGV[0];
$outdir = $ARGV[1];

if (-d $outdir) {
print STDERR "\n";
print STDERR "ERROR: outdir, $outdir, alredy exists!!\n\n";
exit;
} else {
system("mkdir $outdir");
}
}

#$gb2bed = "gb2bed.pl";
$gb2bedgraph = "gb2bedgraph.pl";
$space = "\[\[\:space\:\]\]";

foreach $i (0..$#chrarr) {
$chr = "chr$chrarr[$i]";
$start_genefile = "$outdir/$chr.start_gene.bed";
#system("\grep \"$space$chr$space\" $gbfile | perl $gb2bed - > $start_genefile");
system("\grep \"$space$chr$space\" $gbfile | perl $gb2bedgraph - > $start_genefile");
}

[hoge@boss5 rn4]# perl all_chr.gb2bed.pl rat_uniq_gene.gb rn4start_bedgraph

##----------------------##
#
# gb2bedgraph.pl
#
##---------------------##
#!usr/bin/perl
#
##----------------------------------------##
#
# note batch file for gb2bedgraph.pl
#
##----------------------------------------##


usage :

perl all_chrgb2bed.pl xxxx.gb outdir


for ex.
[hoge@boss5 rn4]# perl all_chr.gb2bed.pl rat_uniq_gene.gb rn4start_bedgraph
[hoge@boss5 rn4start_bedgraph]# ls
chr1.start_gene.bedgraph chr15.start_gene.bedgraph chr20.start_gene.bedgraph chr8.start_gene.bedgraph
chr10.start_gene.bedgraph chr16.start_gene.bedgraph chr3.start_gene.bedgraph chr9.start_gene.bedgraph
..........
[hoge@boss5 rn4start_bedgraph]# head chr1.start_gene.bedgraph


このように

chr  start(0始まりのgene) start+1  高さ1 (適当に全部同じ高さ)


chr1 1306646 1306647 1
chr1 2394947 2394948 1
chr1 5975207 5975208 1
chr1 6051558 6051559 1
...............


こんなフォーマットで作ってみた。

[hoge@boss5 rn4start_bedgraph]# cat chr*.start_gene.bedgraph > all_chr.start_gene.bedgraph

genomegraph uploaded as a custom track in the UCSC genome browser
coverage select


genome graphに遺伝子のスタート座標を乗せよう。

bedgraphでカスタムトラックに追加してから、

インポートでgenome graphへ。

①普通のゲノムブラウザにaddカスタムトラックする。

みなのブログ

②ゲノムブラウザに飛ぶとこんな感じ。

gene startに|||||| ||||| || ||   と並ぶ。

みなのブログ

③先ほどのカスタムトラックに入れたのをゲノムグラフでインポート

Custum tracks rn4gene

table: ct_rn4gene_5451(勝手に番号が付いているみたい)

・ depth を選択。

1000000000 base と入力。

submit する。

みなのブログ



④ゲノムブラウザにrn4gene_011を選択。

みなのブログ

⑤なんかdepthまあ違ったっぽい。


みなのブログ

⑥importへ。

みなのブログ

⑦import table to genome graphsのページで今度はcoverageを選択する。


みなのブログ

⑧ゲノムブラウザでこれを見てみると、

^^^^   ^^^ ^ このような形でcoverageが示される。


みなのブログ



genome graphでもやろう。
[hoge@boss5 rn4]# cd rn4start_bed/
[hoge@boss5 rn4start_bed]# ls
chr1.start_gene.bed chr15.start_gene.bed chr20.start_gene.bed chr8.start_gene.bed
chr10.start_gene.bed chr16.start_gene.bed chr3.start_gene.bed chr9.start_gene.bed
chr11.start_gene.bed chr17.start_gene.bed chr4.start_gene.bed chrX.start_gene.bed
chr12.start_gene.bed chr18.start_gene.bed chr5.start_gene.bed
chr13.start_gene.bed chr19.start_gene.bed chr6.start_gene.bed
chr14.start_gene.bed chr2.start_gene.bed chr7.start_gene.bed
[hoge@boss5 rn4start_bed]# cat chr*.bed > all_chr.start_gene.bed winscpにてmydocにコピーして、 ucscブラウザに行こう。
設定忘れた。。
帰ろう。。
バッチファイルで一括自動処理 systemコマンドを作る時は foreach $i (0..$#chrarr) {
$chr = "chr$chrarr[$i]";
$start_genefile = "$outdir/$chr.start_gene.bed";
print "\grep \"$space$chr$space\" $gbfile | perl $gb2bed - > $start_genefile\n";
} systemのかわりにprint入れて確認すると良い。
また、
$start_genefile = "$outdir/$chr.start_gene.bed"; $outdir/$chr_start_gene.bed
    アンダーバーにすると [hoge@boss5 rn4]# perl all_chr.20100322.plrat_uniq_gene.gbout
grep "[[:space:]]chr1[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/.bed
grep "[[:space:]]chr2[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/.bed
grep "[[:space:]]chr3[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/.bed
・・・・・
                                    .bedの前消えた
アンダーバーが意味を持ってしまったんだろうな。 [hoge@boss5 rn4]# perl all_chr.20100322.plrat_uniq_gene.gbout
grep "[[:space:]]chr1[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/chr1.start_gene.bed
grep "[[:space:]]chr2[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/chr2.start_gene.bed
grep "[[:space:]]chr3[[:space:]]" rat_uniq_gene.gb| perl gb2bed.pl- > out/chr3.start_gene.bed
・・・・・・ もちろん $chr.start_gene.bed でok。
[hoge@boss5 rn4]# perl all_chr.20100322.plrat_uniq_gene.gbout
[hoge@boss5 rn4]# cd out/
[hoge@boss5 out]# ls
chr1.start_gene.bed chr15.start_gene.bed chr20.start_gene.bed chr8.start_gene.bed
chr10.start_gene.bed chr16.start_gene.bed chr3.start_gene.bed chr9.start_gene.bed
chr11.start_gene.bed chr17.start_gene.bed chr4.start_gene.bed chrX.start_gene.bed
chr12.start_gene.bed chr18.start_gene.bed chr5.start_gene.bed
chr13.start_gene.bed chr19.start_gene.bed chr6.start_gene.bed
chr14.start_gene.bed chr2.start_gene.bed chr7.start_gene.bed
[hoge@boss5 out]# head chr10.start_gene.bed
chr10 49 49 BC087103
chr10 72759 72759 BC089943
chr10 635566 635566 BC158662
chr10 666715 666715 BC166736
chr10 762521 762521 X16262
chr10 826803 826803 AB012133
chr10 1492181 1492181 BC166485
chr10 1531025 1531025 BC091175
chr10 3577406 3577406 BC081991
chr10 4249196 4249196 BC079092 outなんかあとでわからなくなるから、ディレクトリ名前付け替え
[hoge@boss5 rn4]# mkdir rn4start_bed
[hoge@boss5 rn4]# mv /home/hoge/gb_file/rn4/out/*.bed /home/hoge/gb_file/rn4/rn4start_bed/
[hoge@boss5 rn4]# cd rn4start_bed/
[hoge@boss5 rn4start_bed]# ls
chr1.start_gene.bed chr15.start_gene.bed chr20.start_gene.bed chr8.start_gene.bed
chr10.start_gene.bed chr16.start_gene.bed chr3.start_gene.bed chr9.start_gene.bed
chr11.start_gene.bed chr17.start_gene.bed chr4.start_gene.bed chrX.start_gene.bed
chr12.start_gene.bed chr18.start_gene.bed chr5.start_gene.bed
chr13.start_gene.bed chr19.start_gene.bed chr6.start_gene.bed
chr14.start_gene.bed chr2.start_gene.bed chr7.start_gene.bed [hoge@boss5 rn4]# cp all_chr.20100322.plall_chr.gb2bed.pl
名前もこれで。 #!/usr/bin/perl
#
# $| = 1; @chrarr = ('1', '2', '3', '4', '5',
'6', '7', '8', '9', '10',
'11', '12', '13', '14', '15',
'16', '17', '18', '19', '20',
'X'); if ($#ARGV < 0) {
print STDERR "\n";
print STDERR "Usage: all_chr.gb2bed.plxxx.gboutdir\n\n";
exit;
} else {
$gbfile = $ARGV[0];
$outdir = $ARGV[1]; if (-d $outdir) {
print STDERR "\n";
print STDERR "ERROR: outdir, $outdir, alredy exists!!\n\n";
exit;
} else {
system("mkdir $outdir");
}
} $gb2bed = "gb2bed.pl";
$space = "\[\[\:space\:\]\]"; foreach $i (0..$#chrarr) {
$chr = "chr$chrarr[$i]";
$start_genefile = "$outdir/$chr.start_gene.bed";
system("\grep \"$space$chr$space\" $gbfile | perl $gb2bed - > $start_genefile");
}