I have a file characterizing genomic regions that looks like this:
chrom chromStart chromEnd PGB
chr1 12874 28371 2
chr1 15765 21765 1
chr1 15795 28371 2
chr1 18759 24759 1
chr1 28370 34961 1
chr3 233278 240325 1
chr3 239279 440831 2
chr3 356365 362365 1
Basically PGB describes the category of the genomic region characterised by its chromosome number (chrom), start (chromStart) and end (chromEnd) coordinates.
I wish to collapse the overlapping regions such that overlapping regions of PGB = 1 and 2 are in a new category, PGB = 3. Output being:
chrom chromStart chromEnd PGB
chr1 12874 15764 2
chr1 15765 24759 3
chr1 24760 28369 2
chr1 28370 28371 3
chr1 28372 34961 1
chr3 233278 239278 1
chr3 239279 240325 3
chr3 240326 356364 2
chr3 356365 440831 3
Basically I wish to obtain an output file which reports unique regions. There are a two criteria.
First, if PGB (column 4) is identical between rows, merge range. eg.
chrom chromStart chromEnd PGB
chr1 1 10 1
chr1 5 15 1
output
chrom chromStart chromEnd PGB
chr1 1 15 1
Second, if PGB is different between rows, chr (column 1) is identical, and the ranges overlap (col2 and 3), report overlapping range as PGB = 3 as well as the ranges unique to their individual categories.
eg.
chrom chromStart chromEnd PGB
chr1 30 100 1
chr1 50 150 2
output
chrom chromStart chromEnd PGB
chr1 30 49 1
chr1 50 100 3
chr1 101 150 2
I hope that illustrates the problem better.
I've created a script that I believe accomplishes this goal.
use List::Util qw(max);
use strict;
use warnings;
# Skip Header row
<DATA>;
# With only 60k rows, going to just slurp all the data
my %chroms;
while (<DATA>) {
chomp;
my ($chrom, $start, $end, $pgb) = split ' ';
# Basic Data Validation.
warn "chromStart is NaN, '$start'" if $start =~ /\D/;
warn "chromEnd is NaN, '$end'" if $end =~ /\D/;
warn "range not ordered, '$start' to '$end'" if $start > $end;
warn "PGB only can be 1 or 2, '$pgb'" if ($pgb ne '1' && $pgb ne '2');
push @{$chroms{$chrom}{$pgb}}, [$start, $end];
}
print "chrom chromStart chromEnd PGB\n";
# Process each Chrom
for my $chrom (sort keys %chroms) {
for my $pgb (sort keys %{$chroms{$chrom}}) {
my $ranges = $chroms{$chrom}{$pgb};
# Sort Ranges
@$ranges = sort {$a->[0] <=> $b->[0] || $a->[1] <=> $b->[1]} @$ranges;
# Combine overlapping and continguous ranges.
# - Note because we're dealing with integer ranges, 1-4 & 5-8 are contiguous
for my $i (0..$#$ranges-1) {
if ($ranges->[$i][1] >= $ranges->[$i+1][0] - 1) {
$ranges->[$i+1][0] = $ranges->[$i][0];
$ranges->[$i+1][1] = max($ranges->[$i][1], $ranges->[$i+1][1]);
$ranges->[$i] = undef;
}
}
@$ranges = grep {$_} @$ranges;
}
# Create pgb=3 for overlaps.
# - Save old ranges into aliases, and then start fresh
my $pgb1array = $chroms{$chrom}{1};
my $pgb2array = $chroms{$chrom}{2};
my @ranges;
# Always working on the first range in each array, until one of the arrays is empty
while (@$pgb1array && @$pgb2array) {
# Aliases to first element
my $pgb1 = $pgb1array->[0];
my $pgb2 = $pgb2array->[0];
# PGB1 < PGB2
if ($pgb1->[1] < $pgb2->[0]) {
push @ranges, [@{shift @$pgb1array}, 1]
# PGB2 < PGB1
} elsif ($pgb2->[1] < $pgb1->[0]) {
push @ranges, [@{shift @$pgb2array}, 2]
# There's overlap for all rest
} else {
# PGB1start < PGB2start
if ($pgb1->[0] < $pgb2->[0]) {
push @ranges, [$pgb1->[0], $pgb2->[0]-1, 1];
$pgb1->[0] = $pgb2->[0];
# PGB2start < PGB1start
} elsif ($pgb2->[0] < $pgb1->[0]) {
push @ranges, [$pgb2->[0], $pgb1->[0]-1, 2];
$pgb2->[0] = $pgb1->[0];
}
# (Starts are equal now)
# PGB1end < PGB2end
if ($pgb1->[1] < $pgb2->[1]) {
$pgb2->[0] = $pgb1->[1] + 1;
push @ranges, [@{shift @$pgb1array}, 3];
# PGB2end < PGB1end
} elsif ($pgb2->[1] < $pgb1->[1]) {
$pgb1->[0] = $pgb2->[1] + 1;
push @ranges, [@{shift @$pgb2array}, 3];
# PGB2end = PGB1end
} else {
push @ranges, [@$pgb1, 3];
shift @$pgb1array;
shift @$pgb2array;
}
}
}
# Append whichever is left over
push @ranges, map {$_->[2] = 1; $_} @$pgb1array;
push @ranges, map {$_->[2] = 2; $_} @$pgb2array;
printf "%-8s %-11s %-11s %s\n", $chrom, @$_ for (@ranges);
}
1;
__DATA__
chrom chromStart chromEnd PGB
chr1 12871 12873 2
chr1 12874 28371 2
chr1 15765 21765 1
chr1 15795 28371 2
chr1 18759 24759 1
chr1 28370 34961 1
chr3 233278 240325 1
chr3 239279 440831 2
chr3 356365 362365 1
And the result:
chrom chromStart chromEnd PGB
chr1 12871 15764 2
chr1 15765 24759 3
chr1 24760 28369 2
chr1 28370 28371 3
chr1 28372 34961 1
chr3 233278 239278 1
chr3 239279 240325 3
chr3 240326 356364 2
chr3 356365 362365 3
chr3 362366 440831 2
Note: I created this mostly as an exercise for myself, but if you use it in any way, please let me know.