Group hash keys contained in a range of numbers
I have a set of cumulative attempts at several different approaches ( approaches 1 to 3
) to determine positions in the genome:
source chromosome1 bp1 chromosome2 bp2
attempt1 2L 5890205 2L 5890720
attempt2 2L 5890205 2L 5890721
attempt1 2L 22220720 2L 22255744
attempt1 3L 15568694 3L 15568866
attempt3 3R 14006279 3R 14008254
attempt1 3R 14006281 3R 14008253
attempt2 3R 14006282 3R 14008254
attempt3 3R 14006286 3R 14008254
attempt1 3R 32060908 3R 32061196
attempt1 3R 32066206 3R 32068392
attempt3 3R 32066206 3R 32068392
attempt2 3R 32066207 3R 32068393
attempt2 X 4574312 X 4576608
attempt1 X 4574313 X 4576607
attempt3 X 4574313 X 4576608
I want to find and group the positions that each attempt has identified, allowing a bit of room for error. For example, I would like to classify the first two lines ...
source chromosome1 bp1 chromosome2 bp2
attempt1 2L 5890205 2L 5890720
attempt2 2L 5890205 2L 5890721
... as one event ( event 1
), which was identified by two different attempts ( attempt1
and attempt2
). I want to classify instances like single events only on different attempts:
- Agree on position
bp1
+/-
5 (i.e. in the window5890200..5890210
) - Define the same
chromosome1
andchromosome2
(2L
) - Agree on position
bp2
+/-
5 (i.e. in the window5890715..5890725
)
I tried to achieve this by using each chromosome and bp as a separate key in the hash
my %SVs;
my $header;
# Make hash
while(<$in>){
chomp;
if ($. == 1){
$header = $_;
next;
}
my ($source, $chromosome1, $bp1, $chromosome2, $bp2) = split;
push @{$SVs{$chromosome1}{$bp1}{$chromosome2}{$bp2}}, $_;
}
}
... and then using a sliding window approach around each bp1 and bp2 value for each line:
my %events;
for my $chr1 ( sort keys %SVs ){
for my $bp1 ( sort { $a <=> $b } keys $SVs{$chr1} ){
my $w1_start = ( $bp1 - 5 );
my $w1_end = ( $bp1 + 5 );
my $window1 = "$w1_start-$w1_end";
for my $chr2 ( sort keys $SVs{$chr1}{$bp1} ){
for my $bp2 ( sort { $a <=> $b } keys $SVs{$chr1}{$bp1}{$chr2} ){
my $w2_start = ( $bp2 - 5 );
my $w2_end = ( $bp2 + 5 );
my $window2 = "$w2_start-$w2_end";
for ( $w1_start .. $w1_end ){
if ($bp1 == $_){
push @{$events{$chr1}{$window1}}, @{$SVs{$chr1}{$bp1}{$chr2}{$bp2}};
}
}
for ( $w2_start .. $w2_end ){
if ($bp2 == $_){
push @{$events{$chr2}{$window2}}, @{$SVs{$chr1}{$bp1}{$chr2}{$bp2}};
}
}
}
}
}
}
print Dumper \%events;
This achieves part of what I want, but I cannot figure out how to get the desired output:
event source chromosome1 bp1 chromosome2 bp2
1 attempt1 2L 5890205 2L 5890720
1 attempt2 2L 5890205 2L 5890721
2 attempt1 2L 22220720 2L 22255744
3 attempt1 3L 15568694 3L 15568866
4 attempt3 3R 14006279 3R 14008254
4 attempt1 3R 14006281 3R 14008253
4 attempt2 3R 14006282 3R 14008254
4 attempt3 3R 14006286 3R 14008254
5 attempt1 3R 32060908 3R 32061196
6 attempt1 3R 32066206 3R 32068392
6 attempt3 3R 32066206 3R 32068392
6 attempt2 3R 32066207 3R 32068393
7 attempt2 X 4574312 X 4576608
7 attempt1 X 4574313 X 4576607
7 attempt3 X 4574313 X 4576608
+3
source to share
1 answer
Below each equivalence class is defined by the last entry added to the equivalence class (based on my understanding of your comment above):
#!/usr/bin/env perl
use strict;
use warnings;
run(\*DATA);
sub run {
my $fh = shift;
my @header = split ' ', scalar <$fh>;
my @events = ([ get_next_event($fh, \@header)]);
while (my $event = get_next_event($fh, \@header)) {
# change the -1 in the second subscript to 0
# if you want to always compare to the first
# event added to the equivalence class
if (same_event($events[-1][-1], $event, 5)) {
push @{ $events[-1] }, $event;
next;
}
push @events, [ $event ];
}
print join("\t", event => @header), "\n";
for my $i (1 .. @events) {
for my $ev (@{ $events[$i - 1] }) {
print join("\t", $i, @{$ev}{@header}), "\n";
}
}
}
sub get_next_event {
my $fh = shift;
my $header = shift;
return unless defined(my $line = <$fh>);
return unless $line =~ /\S/;
my %event;
@event{ @$header } = split ' ', $line;
return \%event;
}
sub same_event {
my ($x, $y, $threshold) = @_;
return if $x->{chromosome1} ne $y->{chromosome1};
return if abs($x->{bp1} - $y->{bp1}) > $threshold;
return if abs($x->{bp2} - $y->{bp2}) > $threshold;
return 1;
}
__DATA__
source chromosome1 bp1 chromosome2 bp2
attempt1 2L 5890205 2L 5890720
attempt2 2L 5890205 2L 5890721
attempt1 2L 22220720 2L 22255744
attempt1 3L 15568694 3L 15568866
attempt3 3R 14006279 3R 14008254
attempt1 3R 14006281 3R 14008253
attempt2 3R 14006282 3R 14008254
attempt3 3R 14006286 3R 14008254
attempt1 3R 32060908 3R 32061196
attempt1 3R 32066206 3R 32068392
attempt3 3R 32066206 3R 32068392
attempt2 3R 32066207 3R 32068393
attempt2 X 4574312 X 4576608
attempt1 X 4574313 X 4576607
attempt3 X 4574313 X 4576608
Output:
event source chromosome1 bp1 chromosome2 bp2 1 attempt1 2L 5890205 2L 5890720 1 attempt2 2L 5890205 2L 5890721 2 attempt1 2L 22220720 2L 22255744 3 attempt1 3L 15568694 3L 15568866 4 attempt3 3R 14006279 3R 14008254 4 attempt1 3R 14006281 3R 14008253 4 attempt2 3R 14006282 3R 14008254 4 attempt3 3R 14006286 3R 14008254 5 attempt1 3R 32060908 3R 32061196 6 attempt1 3R 32066206 3R 32068392 6 attempt3 3R 32066206 3R 32068392 6 attempt2 3R 32066207 3R 32068393 7 attempt2 X 4574312 X 4576608 7 attempt1 X 4574313 X 4576607 7 attempt3 X 4574313 X 4576608
+4
source to share