- 论坛徽章:
- 0
|
本帖最后由 davexu2 于 2014-03-27 16:08 编辑
输入文件1格式，tab键分隔，共10列。其他列均不为空：
chr1 955597
chr1 984302
chr1 1431165
chr1 1558792
chr1 1562080
chr1 3414105
chr1 3416449
chr1 6162054
输入文件2格式，tab键分隔，共15列，其他列不为空：
1 2 chr1 0 10000
2 2 chr1 10000 12820
3 2 chr1 12820 12828
文件2中，4,5列是区间。文件1第2列中，任何一个数字都会在文件2的某个区间中。这个区间的数字，临行的首尾是重复的。
文件1的第一列和文件2的第三列还需要先匹配，才能开始找区间。
理论上来说，文件1中所有行都执行匹配，所以不会执行遍历哈希那一步。但实际上，有大约1半的输出是遍历哈希得到的。不知道bug在何处，求大神指导。
脚本如下：
#!/usr/bin/perl -w
use strict;
use FindBin qw($Bin $Script);
my ($line,@tmps,%hash);
open IN,"$ARGV[0]";
while($line=<IN>){
chomp $line;
@tmps=split /\t/,$line;
$hash{$tmps[0]}{$tmps[1]}=$line;
}
close IN;
open OUT,">$ARGV[2]";
open BZ2,"bzip2 -dc $ARGV[1]|";
while($line=<BZ2>){
next if ($line=~ /^#/);
chomp $line;
@tmps=split /\t/,$line;
next unless $hash{$tmps[2]};
my $res=&mini($hash{$tmps[2]});
next if !defined$res;
print OUT "$res\t$tmps[3]\t$tmps[4]\t".($res>$tmps[3])."\t".($tmps[4]>$res)."\n";
if( exists $hash{$tmps[2]}{$tmps[4]}){
print "$hash{$tmps[2]}{$tmps[4]}\t$tmps[2]\t$tmps[3]\t$tmps[4]\t$tmps[5]\t$tmps[6]\n";
delete $hash{$tmps[2]}{$tmps[4]};
}
elsif(($res>$tmps[3])and($tmps[4]>$res)){
print "$hash{$tmps[2]}{$res}\t$tmps[2]\t$tmps[3]\t$tmps[4]\t$tmps[5]\t$tmps[6]\n";
delete $hash{$tmps[2]}{$res};
redo;
}
}
close BZ2;
#下面的语句理论上不会执行
foreach my $mkey(sort keys %hash){
foreach my $skey(sort keys %{$hash{$mkey}}){
print "$hash{$mkey}{$skey}\n";
}
}
#上面理论上不会执行
sub mini{
my @par=@_;
my $res;
my @tmp=sort {$a<=>$b} keys %{$par[0]};
$res=$tmp[0];
$res;
} |
|
免费注册
发表于 2014-03-27 16:05
