#!/usr/bin/perl # tblastnear.perl =item usage tblastnear.perl: special case of overlapfilter.perl merge with blastnear distance output perl tblastnear -overlaps genes.gff -input tblastn.blexon.gff > tblastnear.table * ? want genes-cds.gff overlaps for anything (remove confusion w/ inside genes) =item input tblastn gff input tblastn.blexon.gff is now result of blast9protstats2.pl: each HSP has {geneID}_(G|S|o)\d+ syntax for multiple gene matches zmore dpse_pg1f.blexonf.gz ##gff-version 3 2 GA_BREN_NSC_ HSP 780 866 60.5 - . Parent=GA_BREN_NSC_50093814_G1;T arget=GA_BREN_NSC_50093814:55-83;align=29 2 dpse_GLEANR_ HSP 1512 1772 174 - . Parent=dpse_GLEANR_4729_G1;Targe t=dpse_GLEANR_4729:1-87;align=87 2 dpse_GLEAN_ HSP 1512 1772 174 - . Parent=dpse_GLEAN_05195_G1;Targe t=dpse_GLEAN_05195:1-87;align=87 2 GA_BREN_NSC_ HSP 1608 1772 111 - . Parent=GA_BREN_NSC_50093814_G1;T arget=GA_BREN_NSC_50093814:1-55;align=55 2 GA_NCBI_GNO_ HSP 1608 1772 111 - . Parent=GA_NCBI_GNO_32303796_G1;T arget=GA_NCBI_GNO_32303796:1-55;align=55 =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 }; use constant SHOWREV => 1; my $BINSIZE = 1000 ; # NO, NOT HERE (in overlapfilter tho) #must be >= NEARDIST #was# 1000; my $NEARDIST = 15000; # common with tandy exon choice my $PCT_ISINSIDE = 0.40; # 0.40; ?? our $debug=1; my ($overlaps,$overlaplist,$markidtype,$input,$itype,$action,$actid,$ok,$mark); my ($overtype,$typeover,$sametypes,$pctover,)= (kOVERLAP,"",1,0); my (%generef, %genebe, %galtids, %genepairs, %modelerrs, %distclass, %nearsteps, %revsteps, %nearpair, %nearrev); my $mrnatype= $ENV{mrnatype} || "mRNA"; my $ngene= $ENV{ngene} || 0; my $eval_max= 1e-10; my $bitscore_min=0; my $align=0.1; # align=0.5 is too restrictive: Dmoj/Est6 tandems fail my $ALLPAIRSCORE=1; my $MIXGROUPS= 0; my $groupok= 1; my $remove2ndgenes= 0; my $hspadd= 0; #?? my $skipgeneoverlaps= 0; my $clustersize= 0; my $clustertop=0; my @clusize; $itype="gff"; $typeover="overlap"; # complex of all: same, near, far ... $action="keep"; # what?? $overlaps= $input=""; my $optok= GetOptions( "input=s", \$input, "overlaps=s", \$overlaps, "eval=s",\$eval_max, "bitscore=s",\$bitscore_min, "align=s",\$align, # not for blast.GFF?? "pctover=i", \$pctover, "NEARDIST=i", \$NEARDIST, "BINSIZE=i", \$BINSIZE, "itype=s", \$itype, "hspadd!", \$hspadd, "remove2ndgenes!", \$remove2ndgenes, "skipgeneoverlaps!", \$skipgeneoverlaps, "debug!", \$debug, #?? MIXGROUPS , groupok # "typeover=s", \$typeover, # "mark=s", \$mark, # "midtype=s", \$markidtype, # return not ID= but other attribute or score # "action=s", \$action, ); die "usage: perl tblastnear -overlaps genes.gff -itype gff -input tblastn[.gff] > tblastnear.table options: -neardist=$NEARDIST (base distance) -itype=gff|blast -bitscore=$bitscore_min : score filter -pctover=0 (for overlap type, min % overlap) -hspadd : combine gene-hsps/location for quality filter (tblastn2gff sorted by gene) -skipgeneoverlaps : ignore matches to genes.gff locations (new only) -remove2ndgenes : remove 2ndary gene match to hsp per group,location (not desired) " unless($optok and (-f $overlaps or $input));# and $action # usage: perl overlapfilter -act keep|drop|mark|markid -mark=terepeat # -overlaps terepeats.gff -typeover overlap|sameloc|samefeat|near # -pctover=0 (for overlap type, min % overlap) # -neardist=$NEARDIST (for near type, base distance) # -midtype=ID|score|Name|... (for markid, attribute or score to mark, ID default) # -itype=gff|blast -input stdin|tandy.gff > tandyfilt.gff $mark ||= "genematch"; $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 =~ /same/) ? kSAMELOC : ($typeover =~ /near/) ? kNEARLOC : kOVERLAP; # $sametypes= ($typeover =~ /feat/) ? 1 : 0; #?? $pctover= $pctover/100.0 if($pctover); # NOT HERE# $BINSIZE= int($NEARDIST*2) if($overtype == kNEARLOC and $NEARDIST > $BINSIZE); #............ my ($noloc,$nskip, $nchain)=(0) x 9; my %noloc=(); # by group $overlaplist={}; if($overlaps) { my $ovh=undef; 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); } else { die "bad -overlaps=undef" if($skipgeneoverlaps); } #............ 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= tblast2near($inh); } elsif($itype =~ /gff/i) { $nr= tblastgff2near($inh); } else { die "cannot read $itype input format"; } close($inh); # unless STDIN ? # warn"#overlaps found=$nr\n" if $debug; OUTPUT: outdist(); #.................. sub getId { my($attrib)= @_; if(!$attrib) { return ""; } elsif($attrib =~ m/ID=([^;\s]+)/) { return $1; } elsif($attrib =~ m/Parent=([^;\s]+)/) { return $1; } return ""; } ## add tblastn table input direct? w/o blast2gff filtering sub tblast2near { my($inh)= @_; my $nover=0; my %grouphsp=(); # filter out multiple hsp-loc/predictor, keep multiple over predictors ## see blast9protstats2.pl ## assume tblastn with query == geneid, source = ref chromosome/scaffold id while(<$inh>){ unless(/^\w/){ next; } # next if(/^(#n |$)/); print and chomp; my @v= split "\t"; my ($qid, $tref, $pctident, $alignment_length, $mismatches, $gap_openings, $q_start, $q_end, $s_start, $s_end, $prob, $bit_score ) = @v; next unless($bit_score); # got also some error logs mixed in # warn,die if all input not in this table size ## hsp2near does this test; move it?? if($bitscore_min) { $nskip++ and next if($bit_score < $bitscore_min); } elsif($eval_max) { $nskip++ and next if($prob > $eval_max); } my($s_strand,$q_strand)= ('+','+'); if ($s_start > $s_end) { $s_strand='-'; ($s_start,$s_end)= ($s_end,$s_start); } if ($q_start > $q_end) { $q_strand='-'; ($q_start,$q_end)= ($q_end,$q_start); } # if($hspadd) { } #?? ## need to change tiloc.qid to keep hsp2near from matching only self ## but need to keep separate hsp-exons of same gene match under 1 id .. as per G,S,o of blast9protstats2 #my $tkey= "$qid:$q_start-$q_end"; # need more : sid.s_start or s_binloc ?? my $tbin= int(($s_start+$s_end)/ (2 * $BINSIZE)); my $tid= "$qid:$q_start-$q_end.$tref:$tbin"; my($tsrc, $ttyp)= ("tblastn","HSP"); #? my $tiloc= join("\t",$tref,$tsrc,$ttyp,$s_start,$s_end,$bit_score,$s_strand,$tid,"blast_hsp"); my $ok= hsp2near( $qid, \%grouphsp, $tiloc, $tref,$s_start,$s_end,$bit_score, ); } return $nover; } sub tblastgff2near { my($inh)= @_; my $nover=0; my %grouphsp=(); # filter out multiple hsp-loc/predictor, keep multiple over predictors my $lasthspid=""; my $lastref; my $genelines=""; my $ngenehsp=0; while(<$inh>){ unless(/^\w/){ next; } # next if(/^(#n |$)/); print and # this is tblastn match location, matching gene id for HSP my($tref,$tsrc,$ttyp,$tb,$te,$tp,$to, $tx,$gffmore)= split"\t"; # warn,die if all input not in this table size my ($iid, $iidhsp); $iidhsp= $iid= getId($gffmore); next unless($iid); ## ^^ ?? collect all HSP/gene at these gene-loc first ## measure n_exons, n_alignsize for quality of match filter ?? ## but need to a. sort by $tsrc, $iid over gene-region 1st and separate diff genes ## unsorted tblastn2gff input should be grouped by query gene: use that to collect all gene-hsps? if($hspadd) { my $inline= $_; if($iidhsp eq $lasthspid) { # and $tref eq $lastref # dont need $genelines .= $_; $ngenehsp++; next; } elsif($ngenehsp>0) { my($kidhsp, $kid, $kref,$ksrc,$ktyp,$kb,$ke,$kp,$ko, @kmore); $kb=-1; $ke=0; $kp=0; my @genelines= split "\n", $genelines; foreach my $gl (@genelines) { my($jref,$jsrc,$jtyp,$jb,$je,$jp,$jo,@jmore)= split"\t",$gl; ($kref,$ksrc,$ktyp,$ko)= ($jref,$jsrc,$jtyp,$jo); $kb= ($kb<0) ? $jb : _min($kb,$jb); $ke= _max($ke,$je); $kp += $jp; @kmore= @jmore; # need only one for iid } $kidhsp= $kid= getId($kmore[-1]); my $kiloc= join("\t",$kref,$ksrc,$ktyp,$kb,$ke,$kp,$ko,$kidhsp,"blast_hsp"); my $ok= hsp2near( $kid, \%grouphsp, $kiloc, $kref,$kb,$ke,$kp, ); #now regen new gff input line tref ... } else { # ?? not here } # .. urk, need to stop here to collect any more hsps $lasthspid= $iidhsp; $genelines= $inline; $ngenehsp=1; $lastref=$tref; } else { my $tiloc= join("\t",$tref,$tsrc,$ttyp,$tb,$te,$tp,$to,$iidhsp,"blast_hsp"); my $ok= hsp2near( $iid, \%grouphsp, $tiloc, $tref,$tb,$te,$tp, ); } } # while in # validmodelerrs(\%genepairloc,\%moderrtemp); return $nover; } sub hsp2near { my($iid, $grouphsp, $tiloc, $tref, $tb, $te, $tp)= @_; my $nover=0; $nskip++, return 0 if($bitscore_min and $tp < $bitscore_min); #? what if score == eval? # what if we are collecting gene-hsps? ## input tblastn IDs have '_(G1|G2|Sn|on)' appended: GA_NCBI_GNO_32410721_G1 $iid =~ s/_(G|S|o)\d+$//; my($igroup)= ($groupok) ? $iid =~ m/^(\D+)/ : ("all"); ## ^ input $tsrc should be igroup; use that? # have multiple predictor hsp per location; keep all or 1 only??? # probably don't want to remove, these are the dupl. gene matches we want to count if($remove2ndgenes) { my $tbin= int(($tb+$te)/1000); # want only hsp-sized bin here return 0 if( $grouphsp->{$igroup}{"$tref.$tbin"}++ ); } my $iloc= $genebe{$iid}; unless($iloc){ $noloc{$igroup}++; return 0; } ; # require all tblast hsp to have valid gene gff my ($iref,$isrc,$ityp,$ib,$ie,$ip,$io,$iid2)= split "\t", $iloc; # 8-col gff line # dont confuse the tref,tb,te genome location with $iid gene's genome loc # my $tiloc= join("\t",$tref,$tsrc,$ttyp,$tb,$te,$tp,$to,$iidhsp,"blast_hsp"); ## step thru J genes at same locus; FIXME: add case of no match gene ?? ## use overlaplist to skip this hsp if it overlaps predicted gene (same group) if($skipgeneoverlaps) { my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $bini (@bins) { $overlaplist->{$tref}{$bini} or next; foreach my $jloc ( split( "\n", $overlaplist->{$tref}{$bini} ) ) { my ($jref,$jsrc,$jtyp,$jb,$je,$jp,$jo,$jid,$jflags)= split "\t", $jloc; next unless($jid); #? error my($jgroup)= ($groupok) ? $jid =~ m/^(\D+)/ : ("all"); next unless( $jgroup eq $igroup); my $isover_hsp= ($jref eq $tref && $tb <= $je && $te >= $jb) ? 1 : 0; if ($isover_hsp) { $nskip++; return 0; } } } } #NOT# my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); ## double dang; this my ib clashes above iid ib #NOT# foreach my $bini (@bins) # was ib if(1) # no bin loop here ... { my @locs=(); # ($tiloc); # .. test case with NO overlaps ; just $tiloc (tblastn hsp) ## unless($skipgeneoverlaps) # not here; skip below if no @locs ##? maybe yes, exclude these overlap genes: want to test only tblast hsp vs iigene? ## wide tb/binsize is finding near genes with NO relation to this hsp, from overlaplist ## want ONLY genes that overlap with this tiloc hsp, that shows similarity ## dont need @bins loop here? # unless($skipgeneoverlaps) { # @locs= split "\n", $overlaplist->{$tref}{$bini} if $overlaplist->{$tref}{$bini}; # } # foreach my $jloc ($tiloc, @locs) my $jloc= $tiloc; { my ($jref,$jsrc,$jtyp,$jb,$je,$jp,$jo,$jid,$jflags)= split "\t", $jloc; # 8-col gff line next unless($jid); #?? error $jflags||=""; my($jgroup)= ($groupok) ? $jid =~ m/^(\D+)/ : ("all"); next unless( $MIXGROUPS || $jgroup eq $igroup); my $grp= $igroup; # "$igroup.$jgroup" #?? # if($skipgeneoverlaps and $jflags !~ /blast_hsp/) { # my $isover_hsp= ($jref eq $tref && $tb <= $je && $te >= $jb) ? 1 : 0; # next if($isover_hsp); # } # ?? this is bad; need to add some location to hash? # want to count gene IDs, not hsps/exons; # but? want to allow multiple gene I x gene J matches (at diff locations)? # switch to $iidhsp for tiloc, so can do muliple matches to iid gene next if($genepairs{$iid}{$jid} ); # not $dupgeneok and () $genepairs{$iid}{$jid}++; # store distance instead of count for pairs? if( $iref ne $jref) { $distclass{$grp}{'scaf'}{$iid}++; if($ALLPAIRSCORE) { $nearpair{"$iid.$jid"}= 399; } # = UNLINKSCORE } else { my ($dist,$rev)= _mindistance( $ib,$ie,$io, $jb,$je,$jo ); # returns dist:1 == _isinside ; dist:2 == small overlap $dist=0 if($dist < 2); # for now; mostly subset gene calls of same larger model if ($dist<=0) { $nover++; $distclass{$grp}{'same'}{$iid}++; if($ALLPAIRSCORE) { $nearpair{"$iid.$jid"}= 199; } # = OVERLAPSCORE; skip this one?? } else { if($dist>$NEARDIST) { $distclass{$grp}{'far'}{$iid}++; $nearsteps{$grp}{999}{$iid}++; if($ALLPAIRSCORE) { $nearpair{"$iid.$jid"}= 299; } # = FARSCORE if($rev) { $revsteps{$grp}{999}{$iid}++; } } else { $distclass{$grp}{'near'}{$iid}++; my $distep= 1 + int( $dist / 1000); $distep= 0 if($dist < 30); my $pairscore= $distep; $nearpair{"$iid.$jid"}= $pairscore; $nearsteps{$grp}{$distep}{$iid}++; if ($rev) { $revsteps{$grp}{$distep}{$iid}++; } } } } } } return 1; # means ok, added } sub outopts { my @ggroup= sort keys %distclass; my $manygroup= @ggroup>1; my $cc = "# "; ## $dumpgenes ? "#t " : ""; my @flags=(); # push(@flags,pct_align => $pctalign); # push(@flags,pct_bitscore => $pctbits); push(@flags,max_e_value => $eval_max); push(@flags,min_bitscore => $bitscore_min) if $bitscore_min; push(@flags,pct_overlap => $pctover) if($pctover); ## push(@flags,pct_inside => $PCT_ISINSIDE) if($PCT_ISINSIDE); push(@flags,remove_2ndgenes => $remove2ndgenes) if($remove2ndgenes); push(@flags,skip_geneoverlaps => $skipgeneoverlaps) if($skipgeneoverlaps); push(@flags,hsp_add => $hspadd) if($hspadd); # push(@flags,pairfilter => $pairfilter) if($pairfilter); # push(@flags,keepgffattr => $keepgffattr) if($keepgffattr); push(@flags,mixed_groups => $MIXGROUPS) if($MIXGROUPS); my %flags=@flags; my $flags= join ",", map{"$_=".$flags{$_}} sort keys %flags; print $cc,"Tandy-blastnear counts of predictor groups\n"; # $ttype per print $cc,"Options: $flags\n"; print $cc,"Data : in=$input, over=$overlaps\n"; # print $cc,join("\t","Group ", "Stat ",@eqnames),"\n"; } sub outdist { my($gscaf,$gfar,$gnear,$gsame, $gpairs)=(0)x9; my($pscaf,$pfar,$pnear)= (0)x9; outopts(); my @grp= sort keys %distclass; @grp= sort keys %noloc unless(@grp); foreach my $grp (@grp) { $gscaf= scalar(keys %{ $distclass{$grp}{'scaf'} } ); $gfar = scalar(keys %{ $distclass{$grp}{'far'} } ); $gnear= scalar(keys %{ $distclass{$grp}{'near'} } ); $gsame= scalar(keys %{ $distclass{$grp}{'same'} } ); my $modelerrs= ""; # # join ",", map{ my $c= $modelerrs{$grp}{$_}; "$_:$c" } sort keys %{ $modelerrs{$grp} }; ## need $grp here : $ngene || my $ngenetotal = scalar( grep /^$grp/, keys %genebe ); # count all mRNA/prot IDs; fixme for full AA subsets $gpairs=scalar( grep /^$grp/, keys %genepairs); my $noloc= $noloc{$grp} || 0; print "source=$grp; " if $grp; # was $source print "ngenes=$ngenetotal; npaired=$gpairs; ngroup=$nchain; noloc=$noloc; skip=$nskip; "; print "clusters=$clustersize..$clustertop; " if($clustersize); print "clusters=all; " unless($clustersize); print "neardist=$NEARDIST\n"; print "modelerrs=$modelerrs\n" if($modelerrs); print " gscaf=$gscaf; gfar=$gfar; gnear=$gnear; gsame=$gsame\n"; # want range of near steps: 1000, .., 5000, 10000, 15000 .. 50000 , far, scaf foreach my $kb (sort {$a<=>$b} keys %{$nearsteps{$grp}}) { my $cnear =scalar(keys %{$nearsteps{$grp}{$kb}}); print " kb$kb=$cnear," } print "\n"; if(SHOWREV) { foreach my $kb (sort {$a<=>$b} keys %{$revsteps{$grp}}) { my $crev =scalar(keys %{$revsteps{$grp}{$kb}}); print " rv$kb=$crev," } print "\n"; foreach my $kb (sort {$a<=>$b} keys %{$revsteps{$grp}}) { my $cnear =scalar(keys %{$nearsteps{$grp}{$kb}}) || 1; my $crev =scalar(keys %{$revsteps{$grp}{$kb}}); my $rp= sprintf("%.3f",$crev/$cnear); print " revp$kb=$rp," # or $crev? both } print "\n"; } print "\n"; } #... end $grp ## FIXME # if($shownear) { print " near-genes:\n"; # change to all dupl pairs: ALLPAIRSCORE == showall # map{ # my $v= $nearpair{$_}; my $pr=1; # print $_,"\t$v\n" if($pr); # } sort keys %nearpair; # } } #................. my $warns=0; sub _min { return ($_[1] < $_[0]) ? $_[1] : $_[0]; } sub _max { return ($_[1] > $_[0]) ? $_[1] : $_[0]; } sub _isinside { #assume overlap? check for g-span >> q-span ## Q is inside T my($tb,$te,$to, $qb,$qe,$qo, $checkswap)= @_; my $tlen= abs($te - $tb) || 1; my $qlen= abs($qe - $qb) || 1; if($checkswap and $qlen > $tlen) { #swap ($qlen,$tlen)=($tlen,$qlen); ($tb,$te,$to, $qb,$qe,$qo)= ($qb,$qe,$qo, $tb,$te,$to); } if($qb > $tb and $qe < $te) { return 1 if( $qlen/$tlen < $PCT_ISINSIDE ); # 0.33; #? use pctover here? } return 0; } sub _mindistance { my($tb,$te,$to, $qb,$qe,$qo)= @_; my $rev= ($to and $qo and ($to eq "-" or $qo eq "-") and ($to ne $qo)) ? 1 : 0; # reversed not counted as overlap my $over= ($tb <= $qe && $te >= $qb) ? 1 : 0; if($over) { my $isin= _isinside($tb,$te,$to, $qb,$qe,$qo, 1); return (1,$rev) if($isin); # special case } if($over and $pctover) { my ($bb,$be)= ( _max($tb,$qb), _min($te,$qe) ); my $maxo= abs($be - $bb); my $leno= _min( abs($qe - $qb), abs($te - $tb)) || 1; $over = 0 if $maxo/$leno < $pctover; return (2,$rev) unless($over); # special case 2; treat different from _isinside? } if($over) { return (0,$rev); } else { # assume gb,ge, qb,qe are ordered my $bd= abs($tb - $qe); # g above q my $ed= abs($qb - $te); # g below q ; dont need gb - qb, ge - qe test return ($ed < $bd) ? ($ed,$rev) : ($bd,$rev); } } sub collect_overlaps { my($gff)= @_; my %overlaps=(); my $nr=0; while(<$gff>){ next unless(/^\w/); my($ref,$src,$typ,$tb,$te,$tp,$to,$tx,$gffmore)= split"\t"; my $gid= getId($gffmore); next unless($te>0 and $gid); ## must have id my $rloc= join("\t", $ref,$src,$typ,$tb,$te,$tp,$to,$gid); if($typ =~ /$mrnatype/) { $genebe{$gid}= $rloc; } #? save gene and cds in overlaps, or just one?? my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { $overlaps{$ref}{$ib} .= $rloc."\n"; # one string better for mem use } $nr++; } warn"#collect_overlaps=$nr\n" if $debug; return \%overlaps; }