#!/usr/bin/perl -w my $progname = $0; $progname =~ s!^.*/!!; my $usage = "$progname: cut out everything in one GFF file that overlaps with a second\n"; $usage .= "\n"; $usage .= "Usage: $progname \n"; $usage .= "\n"; $usage .= "Example (returns all chromosome regions that are NOT in genes):\n"; $usage .= " $progname CHROMOSOMES.gff GENES.gff\n"; $usage .= "\n"; die $usage unless @ARGV == 2; my ($file1, $file2) = @ARGV; my ($byname1, $byname2) = map (load_GFF($_), $file1, $file2); # GFF field 3 = start # GFF field 4 = end for my $seqname (sort keys %$byname1) { next unless exists $byname2->{$seqname}; my @gff2 = sort { $a->[3] <=> $b->[3] || $a->[4] <=> $b->[4] } @{$byname2->{$seqname}}; for my $gff1 (@{$byname1->{$seqname}}) { my $start = $gff1->[3]; for my $gff2 (@gff2) { next if $gff2->[4] < $gff1->[3]; last if $gff2->[3] > $gff1->[4]; if ($gff2->[3] > $start) { print join ("\t", @{$gff1}[0..2], $start, $gff2->[3] - 1, @{$gff1}[5..8]), "\n"; } $start = $gff2->[4] + 1; } if ($start <= $gff1->[4]) { print join ("\t", @{$gff1}[0..2], $start, @{$gff1}[4..8]), "\n"; } } } sub load_GFF { my ($file) = @_; my %byname; local *FILE; open FILE, "<$file"; while (my $gff = ) { chomp $gff; my @f = split /\t/, $gff, 9; push @{$byname{$f[0]}}, \@f; } return \%byname; }