20110406作業履歴 2つのgenomegraphのdataを平均するプログラム

201104062つのgenomegraphのdataを平均するプログラム


##------------------------------------##
#
# averaging 2 genomegraph files
#
##------------------------------------##

usage:

perl all_2files2genomegraph.pl genomegraphfile_1 genomegraph_2 > out.txt


[hoge@boss5 log2]# perl all_2files2genomegraph.pl /home/hoge/aCGH/log2/293/wo_repeatprobe/wo_100k_dist_probe/discart_overlap_new_start/chr_293_100k_dist/chr_293_100k_ave_genomegraph/allchr_293_100k_ave.genomegraph /home/hoge/aCGH/log2/317/wo_repeatprobe/wo_100k_dist_probe/discart_overlap_new_start/chr_317_100k_dist/chr_317_100k_ave_genomegraph/allchr_317_100k_ave.genomegraph > /home/hoge/aCGH/log2/293_317/all_293_317.genomegraph


またgenomegraphにuploadしてみる。
geneと一緒に乗るといいんだけど、2つしか乗らないんだね。
自分で書いた方がいいかなやっぱり。
ますますtreat群とcontrolとの区別がイマイチに

318-292
293-317








#!/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_2files2genomegraph.pl allchr_292_100k_ave.genomegraph allchr_293_100k_ave.genomegraph\n\n";
exit;
} else {
$avefile1 = $ARGV[0];
$avefile2 = $ARGV[1];
}

#-----------------
# read avefiles
#-----------------

undef(@avechr1);
undef(@avepos1);
undef(@aveval_1);
undef(@avechr2);
undef(@avepos2);
undef(@aveval2);
open(AVEFILE1, "<$avefile1") || die "Can't open $avefile1";
while (<AVEFILE1>) {
chomp;
@col = split(/\s+|\t/, $_);
push(@avechr1, $col[0]);
push(@avepos1, $col[1]);
push(@aveval_1, $col[2]);
}
close(AVEFILE1);


open(AVEFILE2, "<$avefile2") || die "Can't open $avefile2";
while (<AVEFILE2>) {
chomp;
@col2 = split(/\s+|\t/, $_);
push(@avechr2, $col2[0]);
push(@avepos2, $col2[1]);
push(@aveval2, $col2[2]);
}
close(AVEFILE2);

foreach $i (0..$#chrarr) {
$chr = "chr$chrarr[$i]";
foreach $h (0..$#avepos1) {
if (($chr eq $avechr1[$h]) && ($chr eq $avechr2[$h]) && ($avepos1[$h] == $avepos2[$h])) {
$ave_aveval = ($aveval_1[$h] + $aveval2[$h])/2;
print "$chr\t$avepos1[$h]\t$ave_aveval\n";
}
}
}

20110405作業履歴プローブをaveしてgenomegraphに乗せるバッチファイル

20110405作業履歴100kdistで取ったプローブをaveしてgenomegraphに乗せるバッチファイル ##------------------------------------##
#
# batch file for ave2genomegraph
#
##------------------------------------##
usage: perl /home/hoge/aCGH/log2/all_chr.ave2genomegraph.plallchr_292_100k_dist.txt chr_292_10
0k_ave_genomegraph

ARGV[0] は log2file
ARGV[1] は outdir 出来あがったディレクトリに移動して、

cat chr*_292_100k_dist.txt > allchr_292_100k_dist.txt upload ucsc genomegraph





2種しか座標乗せられないのが残念だ






all_chr.ave2genomegraph.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.ave2genomegraph.plchr*_nnn_nnnk_dist.txt chr_nnn_nnnk_ave_genomegraph\n\n";
exit;
} else {
$log2file = $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");
}
} $microarray = (split(/\_/, $log2file))[1];
$dist = (split(/\_/, $log2file))[2];
$ave2genomegraph = "/home/hoge/aCGH/log2/ave2genomegraph.pl";
$space = "\[\[\:space\:\]\]"; foreach $i (0..$#chrarr) {
$chr = "chr$chrarr[$i]";
$ave_genomegraph = "$outdir/$chr\_$microarray\_$dist\_ave.genomegraph";
system("cut -f 5-9 $log2file | \grep \"$space$chr$space\" - | perl $ave2genomegraph - > $ave_genomegraph");
}