2月26日の作業履歴

単なるメモ。dos2unix文字コード変換　　コマンド。。。。。。。。。。

とりあえず相変わらずだけどプローブを100k間隔で取ってくるプログラムを作る。

その前に、プローブ間隔でこの間隔から大幅に離れている場所を出そう。

そして、あとでゲノムブラウザに載せられるといいな。

プローブ間隔を出すのにstartから前のプローブのendを引いて間隔を出すプログラム。
#!/usr/bin/perl
#
# grep "[[:space:]]chr*[[:space:]]" agilent244a.txt | perl probe_dist.2010.pl -
# $start[$i+1] - $end[$i] ######################################

$| = 1;

$infile = $ARGV[0];
open(INFILE, "< $infile") || die "Can't open $infile";
while(<INFILE>) {
chomp;
push(@agiarray, $_);
push(@start, (split(/\t/, $_))[2]);
push(@end, (split(/\t/, $_))[3]);
push(@id, (split(/\t/, $_))[4]);
#print "$start $end $id\n";
}

foreach $i (0..$#end) {
if ($i !=0) {
$distance = $start[$i] - $end[$i-1];
$tempid = $id[$i];
print "$tempid\t$id[$i-1]\t$distance\n";
} else {
$distance = $start[$i] - 0;
print "$id[0]\t0\t$distance\n";
}
}

~
[minako1@minako agilent]# grep "[[:space:]]chr10[[:space:]]" agilentCgh244a.txt | perl probe_dist.2
.pl - | head
A_86_P04312176 0 181
A_86_P04312201 A_86_P04312176 8129
A_86_P04312216 A_86_P04312201 4206
A_86_P04312225 A_86_P04312216 4739
A_86_P04312228 A_86_P04312225 8540
A_86_P04312237 A_86_P04312228 3268
A_86_P04312247 A_86_P04312237 3244
A_86_P04312257 A_86_P04312247 2761
A_86_P04312278 A_86_P04312257 4771
A_86_P04312293 A_86_P04312278 4182

調べたいchrを選択してからプログラムに入れる。

probeが無い部分はどれだけか？と思うので、
こんなプログラム。

"20100219.pl" 30 lines, 658 characters written
#!/usr/bin/perl
#
#

$| = 1;

$infile = $ARGV[0];
$select_dist = $ARGV[1];
open(INFILE, "< $infile") || die "Can't open $infile";
while(<INFILE>) {
chomp;
push(@agiarray, $_);
push(@start, (split(/\t/, $_))[2]);
push(@end, (split(/\t/, $_))[3]);
push(@id, (split(/\t/, $_))[4]);
#print "$start $end $id\n";
}

foreach $i (0..$#end) {
if ($i !=0) {
$distance = $start[$i] - $end[$i-1];
$tempid = $id[$i];
print "$tempid\t$id[$i-1]\t$distance\n" if ($select_dist < $distance);
} else {
$distance = $start[$i] - 0;
print "$id[0]\t0\t$distance\n" if ($select_dist < $distance);
}
}

[minako1@minako agilent]# grep "[[:space:]]chr10[[:space:]]" agilentCgh244a.txt | perl longerthan_selectdist.pl - 100000 | head
A_86_P04317089 A_86_P04316823 166272
A_86_P04320226 A_86_P04320204 151042
A_86_P04344161 A_86_P04344114 152947
A_86_P04364781 A_86_P04364780 110002
A_86_P04415780 A_86_P04415772 323738
A_86_P04418327 A_86_P04418324 173550
A_86_P04440080 A_86_P04440075 317143
A_86_P04462526 A_86_P04462516 104600
A_86_P04468130 A_86_P04468127 124428
A_86_P04468740 A_86_P04468456 101972
[minako1@minako agilent]#
~
後で100k間隔でプローブを取る準備をしているのだ。
NNNNNなのかどうか知らないけど、chr10ではこのように間隔の広い部分が発見。

これらプローブがどこにあるのか出してくれるプログラムを書こう。
#!/usr/bin/perl
#
#

$| = 1;

$infile = $ARGV[0];
$agilent = $ARGV[1];
open(INFILE, "< $infile") || die "Can't open $infile";
while(<INFILE>) {
chomp;
push(@pre, (split(/\t/, $_))[0]);
push(@post, (split(/\t/, $_))[1]);
}
close(INFILE);

@allid = push(@pre, @post);
print "@allid\n";

open(AGILENT, "< $agilent") || die "Can't open $agilent";
while(<AGILENT>) {
chomp;
push(@agiarray, $_);
push(@chr, (split(/\t/, $_))[1]);
push(@start, (split(/\t/, $_))[2]);
push(@end, (split(/\t/, $_))[3]);
push(@id, (split(/\t/, $_))[4]);
}
close(AGILENT);

foreach $i (0..$#id) {
$tempid = $id[$i];
# print "\$tempid = $tempid\n";
foreach $h (0..$#allid) {
# print "\$i = $i\n";
$tempallid = $allid[$h];
if ($tempid eq $tempallid) {
print "\$tempid= $tempid = \$allid[$h] = $allid[$h]\n";
# print "$chr[$i]\t$start[$i]\t$end[$i]\$id[$i]\n";
}
}
}

[minako1@minako agilent]# grep "[[:space:]]chr10[[:space:]]" agilentCgh244a.txt | perl longerthan_s
electdist.pl - 100000 | perl 20100226.pl - agilentCgh244a.txt | head
34

途中経過だけど、何かわからない。
また、明日。