#!/usr/bin/perl # overlapbestgene.perl =item overlapbestgene special case of overlap detection: -overlap = tandy tandem gene 'models' (match), maybe also HSP/exons -input = gene predictions (many predictors) output = closest gene prediction(s?) matching tandy models =item from overlapfilter perl overlapfilter -act keep|drop|mark -mark=terepeat -overlaps terepeats.gff -itype=gff|blast -input stdin|tandy.gff > tandyfilt.gff apply identified/predicted TE repeat GFF to tandy or blast results generalized to filter or mark GFF for any GFF overlap data (protein & EST matches, repeats) gzcat *exons_tandy6j.gff.gz | perl $td/overlapfilter.perl \ -act mark -mark=proteinhsp -over *prot9-hsp.gff.gz -in stdin > dere_exons_tandy6jm.gff =cut use strict; use warnings; use Getopt::Long; use constant SAMEBASE => 0; # for _sameloc, slop allowed in loca == locb use constant { ACT_DROP=>1, ACT_KEEP=>2, ACT_MARK=>3, ACT_MARK_WITH_ID=>4 }; use constant { kOVERLAP=>1, kSAMELOC=>2, kNEARLOC=>3, kBESTGENE=>4, kBESTEXONS=>5, }; use constant BESTEXON_SCORE => 2; our $BINSIZE = 1000 ; #was# 5000; our $NEARDIST = 500; # was 15k; needs to be < BINSIZE our $debug=1; my ($overlaps,$overlaplist,%geneexons,$markidtype,$input,$itype,$action,$actid,$ok,$mark); my ($overtype,$typeover,$sametypes,$pctover,)= (kBESTGENE,"",1,0); my $exontypes="CDS,exon,match_part,HSP"; my $mrnatypes="mRNA,match"; # fixme ?? need this my $correct_partial_hsp=1; my $optok= GetOptions( "overlaps=s", \$overlaps, "typeover=s", \$typeover, "exontypes=s", \$exontypes, "mrnatypes=s", \$mrnatypes, "pctover=i", \$pctover, "NEARDIST=i", \$NEARDIST, "BINSIZE=i", \$BINSIZE, "mark=s", \$mark, "input=s", \$input, "itype=s", \$itype, "midtype=s", \$markidtype, # return not ID= but other attribute or score "action=s", \$action, "debug!", \$debug, ); die " usage: perl overlapbestgene -act keep|mark|markid -mark=bestgene -overlaps tandymatch.gff -typeover bestgene|bestexons -exontypes $exontypes -itype=gff|blast -input stdin|tandy.gff > tandyfilt.gff " unless($optok and $action and (-f $overlaps or -f $input)); # -neardist=$NEARDIST (for near type, base distance) # -midtype=ID|score|Name|... (for markid, attribute or score to mark, ID default) # -pctover=0 (for overlap type, min % overlap) ## for matching IDs b/n updates; need also match type: gene/mRNA/exon/CDS/... $mark ||= "bestgene"; $itype ||= "gff"; $actid= ($action =~ /keep/) ? ACT_KEEP : ($action =~ /mark/) ? ACT_MARK : ACT_DROP; $actid= ACT_MARK_WITH_ID if($actid == ACT_MARK && $action =~ /id/i); $overtype= ($typeover =~ /bestexon/) ? kBESTEXONS : ($typeover =~ /best/) ? kBESTGENE : ($typeover =~ /same/) ? kSAMELOC : ($typeover =~ /near/) ? kNEARLOC : ($typeover =~ /over/) ? kOVERLAP : kBESTGENE; $sametypes= ($typeover =~ /feat/) ? 1 : 0; #?? $pctover= $pctover/100.0 if($pctover); $BINSIZE= int($NEARDIST*2) if($overtype == kNEARLOC and $NEARDIST > $BINSIZE); %geneexons=(); # for kBESTEXONS $exontypes =~ s/[,; ]+/|/g; $mrnatypes =~ s/[,; ]+/|/g; # allow -overlap gff == stdin if -input is a file? ##local *R; my $ovh; ##= *OVR; #$ok = ($overlaps =~ /.gz$/) ? open(R,"gunzip -c $overlaps |") : open(R,$overlaps); if($overlaps =~ /.gz$/) { $ok= open(OVR,"gunzip -c $overlaps |"); $ovh= *OVR; } elsif($overlaps =~ /^(stdin|-)/) { $ovh= *STDIN; $ok=1; } else { $ok= open(OVR,$overlaps); $ovh= *OVR; } die "bad -overlaps=$overlaps" unless($ok); $overlaplist= collect_overlaps($ovh); close($ovh); my $inh= *STDIN; $ok = ($input =~ /.gz$/) ? open($inh,"gunzip -c $input |") : ($input =~ /^(stdin|-)/) ? $inh= *STDIN : open($inh,$input); die "bad -input=$input" unless($ok); my $nr=0; if($itype =~ /blast/i) { $nr= filter_blast($inh); } else { $nr= filter_gff($inh); } warn"#overlaps found=$nr\n" if $debug; #.................. sub filter_gff { my($inh)= @_; my $nr=0; my @inlines; while(<$inh>){ unless(/^\w/){ next if(/^(#n |$)/); print and next; } my($ref,$src,$typ,$tb,$te,$tp,$to, @gffmore)= split "\t"; next unless(@gffmore); # for kBESTEXONS, require input sorted by loc? i.e. genes > exons my $gid=""; if($gffmore[-1] =~ m/\bID=([^;]+)/) { $gid=$1; } elsif($gffmore[-1] =~ m/\bParent=([^;]+)/) { $gid=$1; } if( BESTEXON_SCORE == 1 and $overtype == kBESTEXONS and $typ =~ /^($exontypes)$/) { next unless($gid); # push exons into gene/match overlaps .. need %geneexons hash my $rloc= [$tb,$te,$ref,$gid,$to,$typ]; # change to string; save mem push( @{$geneexons{$gid}}, $rloc); # use same geneexons? next; } elsif( BESTEXON_SCORE == 1 and $overtype == kBESTEXONS) { push(@inlines, $_) if($typ =~ /^($mrnatypes)$/); next; } elsif( BESTEXON_SCORE > 1 and $overtype == kBESTEXONS and $typ =~ /^($exontypes)$/) { next; } my $isover= ($overtype == kBESTEXONS) ? bestexons($ref,$tb,$te,$to,$typ) : bestgene($ref,$tb,$te,$to,$typ); if($isover) { # now $isover == ID of match $nr++; next if($actid == ACT_DROP); if($actid == ACT_MARK) { s/$/;$mark=1/; } elsif($actid == ACT_MARK_WITH_ID) { s/$/;$mark=$isover/ ; } } else { next if($actid == ACT_KEEP); } print; } if(BESTEXON_SCORE == 1 and $overtype == kBESTEXONS) { while($_ = shift(@inlines)) { my($ref,$src,$typ,$tb,$te,$tp,$to, @gffmore)= split"\t"; my $gid=""; if($gffmore[-1] =~ m/\bID=([^;]+)/) { $gid=$1; } elsif($gffmore[-1] =~ m/\bParent=([^;]+)/) { $gid=$1; } my $texons= $geneexons{$gid}; my $isover= bestexons($ref,$tb,$te,$to,$typ,$gid,$texons); if($isover) { # now $isover == ID of match $nr++; next if($actid == ACT_DROP); if($actid == ACT_MARK) { s/$/;$mark=1/; } elsif($actid == ACT_MARK_WITH_ID) { s/$/;$mark=$isover/ ; } } else { next if($actid == ACT_KEEP); } print; } } return $nr; } sub filter_blast # ncbi format=8,9 blast table { my($inh)= @_; my $nr=0; while(<$inh>){ unless(/^\w/){ print and next; } my($qid,$ref,$pid,$align,$xa,$xb,$qb,$qe,$tb,$te,@bmore)= split"\t"; # fix blast loc swap for orient my $to='+'; ($tb,$te,$to)= ($te,$tb,'-') if($tb>$te); my $typ="HSP"; my $isover= ($overtype == kBESTGENE) ? bestgene($ref,$tb,$te,$to,$typ) : ($overtype == kSAMELOC) ? sameloc($ref,$tb,$te,$to,$typ) #($ref,$tb,$te,$to,$typ) : ($overtype == kNEARLOC) ? nearloc($ref,$tb,$te) : overlaps($ref,$tb,$te); if($isover) { $nr++; next if($actid == ACT_DROP); if($actid == ACT_MARK) { s/$/\t$mark=1/; } elsif($actid == ACT_MARK_WITH_ID) { s/$/\t$mark=$isover/;} } else { next if($actid == ACT_KEEP); } print; } return $nr; } my $warns=0; # handle array of input sub _min { return ($_[1] < $_[0]) ? $_[1] : $_[0]; } sub _max { return ($_[1] > $_[0]) ? $_[1] : $_[0]; } sub overlaps { my($ref,$tb,$te)= @_; my @lid; return 0 unless($overlaplist->{$ref}); # warn join",",("overlap",$ref,$tb,$te),"\n" if $debug and $warns++<10; my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaplist->{$ref}{$ib} or next; my @locs= @{$overlaplist->{$ref}{$ib}}; foreach my $rloc (@locs) { my ($lb,$le,$lid)= @{$rloc}[0,1,3]; my $over= ($tb <= $le && $te >= $lb) ? 1 : 0; # add option to screen out trival 5% UTR overlaps .. pctoverlap=i if($over and $pctover) { # my $maxo= _max( $le - $tb, $te - $lb); # wrong my ($bb,$be)= ( _max($tb,$lb), _min($te,$le) ); my $maxo= abs($be - $bb); my $leno= _min( abs($le - $lb), abs($te - $tb)) || 1; $over = 0 if $maxo/$leno < $pctover; } push @lid, $lid if($over); # ^^ collect *all* overlap ids as list to return } } if(@lid) { my %lid= map{$_,1}@lid; return join",", sort keys %lid; } return 0; } sub nearloc { my($ref,$tb,$te)= @_; my @lid; return 0 unless($overlaplist->{$ref}); # warn join",",("nearloc",$ref,$tb,$te),"\n" if $debug and $warns++<10; my $tm= int(($tb+$te)/2); # change to min-distance from ends my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaplist->{$ref}{$ib} or next; my @locs= @{$overlaplist->{$ref}{$ib}}; foreach my $rloc (@locs) { my ($lb,$le,$lid)= @{$rloc}[0,1,3]; # my $lm= int(($lb+$le)/2); # push @lid, $lid if (abs($lm - $tm) < $NEARDIST); my $over= ($tb <= $le && $te >= $lb) ? 1 : 0; if($over and $pctover) { my ($bb,$be)= ( _max($tb,$lb), _min($te,$le) ); my $maxo= abs($be - $bb); my $leno= _min( abs($le - $lb), abs($te - $tb)) || 1; $over = 0 if $maxo/$leno < $pctover; } next if($over); # skip any overlaps ?? pctover? my $bd= abs($tb - $le); my $ed= abs($lb - $te); my $mind= ($ed < $bd) ? $ed : $bd; push @lid, $lid if ($mind < $NEARDIST); } } if(@lid) { my %lid= map{$_,1}@lid; return join",", sort keys %lid; } return 0; } # bestexons: improve bestgene using exons (CDS, match_part/HSP) # assume input=match,HSP for tandems; overlap=CDS-exons only w/ parent IDs # ** tandy HSPs are partial exon matches (of dupl region), need some slack # in measuring best overlap CDS with such partial match: CDS/exon extending beyond should # match fully? sub bestexons { my($ref,$tb,$te,$to,$typ,$geneid,$texons)= @_; my @lid; return 0 unless($overlaplist->{$ref}); $to ||= '+'; ($tb,$te,$to)= ($te,$tb,'-') if($tb>$te); my %genescore; # my $texons= $texons{$geneid}; #?? my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaplist->{$ref}{$ib} or next; my @locs= @{$overlaplist->{$ref}{$ib}}; # note here overlaplist == genes, use id to get geneexons foreach my $rloc (@locs) { my ($lb,$le,$lid,$lo,$ltyp)= @{$rloc}[0,1,3,4,5]; next unless($to eq $lo); next if ($genescore{$lid}); my $ovexons= $geneexons{$lid} or next; # need to align target-exons and overexons and count overlaps, score=non-overlaps # see tdprotnear2.perl : _same_exons if (BESTEXON_SCORE == 1) { # NOTE THIS SCORE high is best $genescore{$lid}= _same_exons($texons, $ovexons); # note order in is important } elsif(BESTEXON_SCORE > 1) { $genescore{$lid}= _same_exon2( $tb,$te, $ovexons); # note order in is important } } } # NOTE THIS SCORE high is best my @bestid = sort{$genescore{$b} <=> $genescore{$a}} keys %genescore; # smallest score; but want all near same return 0 unless(@bestid); my $sc0= $genescore{$bestid[0]}; # any range 0..100000 my $sccut= $sc0 - 50; if(0 and @bestid>1) { # score range should be option my $scmax= $genescore{$bestid[-1]}; $sccut = 0.10 * abs($sc0 - $scmax); $sccut = $sc0 - _min( 100, $sccut); } foreach my $id (@bestid) { my $sc=$genescore{$id}; last if ($sc < $sccut); push(@lid,$id); } return join ",",@lid; } sub _same_exon2 { my( $tb,$te, $bex)= @_; ## here $tb,te == tandy gene match span, $bex == gene predicted exons ## also reduce score by distance outside of tb,te; i.e. skipover exon my $nbx= @$bex; my $n= $nbx; my $lena= abs($te-$tb); my $lenb= 0; for(my $i=0; $i<$nbx; $i++) { $lenb += abs($bex->[$i]->[1] - $bex->[$i]->[0]); } my $nalign= 0; for(my $i=0; $i<$n; $i++) { my($bb,$be)= ($bex->[$i]->[0], $bex->[$i]->[1]); my $na=0; if($be < $tb) { # outside $na= $be - $tb; # neg } elsif($bb > $te) { # outside $na= $te - $bb; # neg } else { $na= _max(0, _min($te, $be) - _max($tb, $bb)); } $nalign += $na; $lena -= $na; $lenb -= $na; } return $nalign - ($lena + $lenb); } sub _same_exons { my( $aex, $bex)= @_; # aex, bex are exon list [$b,$e] # exons are sorted by loc ** REVERSE sorted when -strand # note aex, bex order important here; use only aex length as score base # NO, need to reduce score where blen >> overlap # but need to reduce score where blen << alen?, i.e. b is subset of a, # overlap==blen << alen should score below overlap==blen==alen # ?count 2x each overlap base? 2xolap - blen - alen = 0 # * reduce alen,blen for each overlap base olap - blen_left - alen_left > 0 for good olap # ** tandy HSPs (aex) are partial exon matches (of dupl region), need some slack # in measuring best overlap CDS (bex) with such partial match: # CDS/exon extending beyond should score fully? # ** tandy HSPs not quite good enough to match apparent tandem genes? # try scoring tandy gene match region: length of gene-cds inside - gene-cds outside my $nax= @$aex; my $nbx= @$bex; my $n= _min($nax,$nbx); my $lena=0; for(my $i=0; $i<$nax; $i++) { $lena += abs($aex->[$i]->[1] - $aex->[$i]->[0]); } my $lenb=0; for(my $i=0; $i<$nbx; $i++) { $lenb += abs($bex->[$i]->[1] - $bex->[$i]->[0]); } return 0 - ($lena + $lenb) if($n<1); my @xs=(); ## a:------ i index is bad for case of offset genes ## b: ------ ## BUT this can include Skipovers, e.g. ## a: ---.. ...--| :a ## b1: -----| -----| :b2 ## look for alignment of exons first my($ia,$ib)= (0,0); my($ab,$ae)= ($aex->[$ia]->[0], $aex->[$ia]->[1]); my($bb,$be)= ($bex->[$ib]->[0], $bex->[$ib]->[1]); ## REV orient fix my $rev=0; if($n>1) { $rev=($ab > $aex->[1]->[0]) ? 1:0; } if ( ($rev ? $ae<$bb : $be<$ab) ) { for( ; $ib<$nbx; $ib++) { ($bb,$be)= ($bex->[$ib]->[0], $bex->[$ib]->[1]); last if( ($rev ? $ae>$bb : $be>$ab) ); } } elsif( ($rev ? $be<$ab : $ae<$bb) ) { for( ; $ia<$nax; $ia++) { ($ab,$ae)= ($aex->[$ia]->[0], $aex->[$ia]->[1]); last if( ($rev ? $be>$ab : $ae>$bb) ); } } my $nalign= 0; $n= _min($nax - $ia, $nbx - $ib); for(my $i=0; $i<$n; $i++) { # note these may be reverse sorted, if gene is rev ($ab,$ae)= ($aex->[$i+$ia]->[0], $aex->[$i+$ia]->[1]); ($bb,$be)= ($bex->[$i+$ib]->[0], $bex->[$i+$ib]->[1]); # my $na= _min(_max( 0, $be - $ab), _max(0, $ae - $bb)); # $na= _min($na, _min($ae-$ab, $be-$bb)); my $na= _max( 0, _min( $be, $ae) - _max( $bb, $ab) ); ## fixme for blen > alen if( $correct_partial_hsp and $ab >= $bb and $ae <= $be) { $na= abs($be-$bb); } $nalign += $na; $lena -= $na; $lenb -= $na; } return $nalign - ($lena + $lenb); ## return score == overlapped bases - non-overlapped bases ( bad < 0 > good) } # bestgenes: want best overlap, without excess gene span (skipov, dublup errors) # should use CDS-exons/match_part overlaps to count best # want instead of call per predict gene, to call per overlap match location # but find/return best matching predict gene sub bestgene { my($ref,$tb,$te,$to,$typ)= @_; my @lid; return 0 unless($overlaplist->{$ref}); $to ||= '+'; ($tb,$te,$to)= ($te,$tb,'-') if($tb>$te); my %genescore; my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaplist->{$ref}{$ib} or next; my @locs= @{$overlaplist->{$ref}{$ib}}; foreach my $rloc (@locs) { my ($lb,$le,$lid,$lo,$ltyp)= @{$rloc}[0,1,3,4,5]; next unless($to eq $lo); my $score= abs($tb-$lb) + abs($te-$le) ; # best = min(score) $genescore{$lid}= $score; #? } } my @bestid = sort{$genescore{$a} <=> $genescore{$b}} keys %genescore; # smallest score; but want all near same return 0 unless(@bestid); my $sc0= $genescore{$bestid[0]}; # any range 0..100000 my $sccut= $sc0 + 50; if(@bestid>1) { # score range should be option my $scmax= $genescore{$bestid[-1]}; $sccut = 0.10 * abs($scmax - $sc0); $sccut = $sc0 + _min( 100, $sccut); } ## fixme near scores foreach my $id (@bestid) { my $sc=$genescore{$id}; last if ($sc > $sccut); push(@lid,$id);} return join ",",@lid; } sub sameloc { my($ref,$tb,$te,$to,$typ)= @_; my @lid; return 0 unless($overlaplist->{$ref}); # warn join",",("sameloc",$ref,$tb,$te,$to,$typ),"\n" if $debug and $warns++<10; $to ||= '+'; ($tb,$te,$to)= ($te,$tb,'-') if($tb>$te); my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaplist->{$ref}{$ib} or next; my @locs= @{$overlaplist->{$ref}{$ib}}; foreach my $rloc (@locs) { my ($lb,$le,$lid,$lo,$ltyp)= @{$rloc}[0,1,3,4,5]; # [$tb,$te,$ref,$gid,$to,$typ] # also need match ft types?, orient, next if($sametypes and not($typ eq $ltyp and $to eq $lo)); #? push @lid, $lid if(abs($tb-$lb) <= SAMEBASE and abs($te-$le) <= SAMEBASE); } } if(@lid) { my %lid= map{$_,1}@lid; return join",", sort keys %lid; } return 0; } sub collect_overlaps { my($gff)= @_; my %overlaps=(); my $nr=0; while(<$gff>){ next unless(/^\w/); chomp; my($ref,$src,$typ,$tb,$te,$tp,$to,@gffmore)= split"\t"; $nr++; my($gid); if($markidtype) { if($markidtype =~ /^score/i) { $gid=$tp; } elsif($gffmore[-1] =~ m/\b$markidtype=([^;]+)/) { $gid=$1; } } else { if($gffmore[-1] =~ m/\bID=([^;]+)/) { $gid=$1; } elsif($gffmore[-1] =~ m/\bParent=([^;]+)/) { $gid=$1; } } unless(defined $gid) { $gid = "N".$nr; } my $rloc= [$tb,$te,$ref,$gid,$to,$typ]; # change to string; save mem if($overtype == kBESTEXONS and $typ =~ /^($exontypes)$/) { # push exons into gene/match overlaps .. need %geneexons hash push( @{$geneexons{$gid}}, $rloc); next; } elsif($overtype == kBESTEXONS and $typ !~ /^($mrnatypes)$/) { next; } my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { push( @{$overlaps{$ref}{$ib}}, $rloc); } } warn"#collect_overlaps=$nr\n" if $debug; return \%overlaps; } # sub _isoverlap { # my($gb,$ge, $qb,$qe)= @_; # return ($gb <= $qe && $ge >= $qb) ? 1 : 0; # }