#!/usr/bin/perl # tandynear.perl =item notes tandy, count blast self/near match locations per predictor answer this: if predictor has a near/far exon match, does it have also have another prediction at that location? this should be added to tandemgenes; tandystats.pm ? =item usage perl tandynear.perl [ options ] scaffold_6308/dmoj_caf060210_exons.nr.blatf8 # presumes fasta dmoj_caf060210_exons.nr opts: -[no]group; -skipids=listfile; -keeponly=IDgroup -debug=0,1,2; -[no]ignorepartialexon; -neardist=8000; -nearsteps=3; -minalign=0.5; -mineval=1e-05; dgbook% perl $td/tandynear2.perl -debug=1 dpulex1_exons.nr.blatf8 # exon_fasta_ids ref=scaffold_4 ngroup=3; ngenes=1676; ntrans=1788; nexons=11860; naltexons=12795 # poor ids=206008; skipped ids=0 Tandy exon match types per predictor group Group Total Same Near8k Near16k Near48k Far Inside ... NCBI_GNO_ 18450 8211 1768 1012 601 6271 587 NCBI_GNO_ freq 1.000 0.215 0.123 0.073 0.764 0.071 NCBI_GNO_ found 1.000 0.411 0.453 0.383 0.144 0.148 NCBI_GNO_ anyfnd 1.000 0.450 0.491 0.404 0.208 0.496 NCBI_GNO_ fnd/any 1.000 0.914 0.922 0.947 0.694 0.299 compare to tandy outputs, esp. tandemgenes.pl -act find, that processes same data and marks as eq same, near, far, .. per alt-exon. Reconcile different results. cat dpulex1_tandy6j.find | env grp=NCBI_GNO perl -ne\ 'BEGIN{$gr=$ENV{grp}} while(m/($gr\w+)\.\d+\:[\d+\-]+\=([\d-]+)/g){($d,$e)=($1,$2); print"$e\n";}'\ | sort | uniq -c 2366 -2 ; 1339 0 ; 7134 1 Group Total Same Near8k Near16k Near48k Far Inside NCBI_GNO_ -- 7134 2366 ... 1339 -- NCBI_GNO_ -- 1.0 0.33 0.19 =item further on perl $td/tandynear2.perl -refgr EISE_CGW \ -keep 'GK_NCBI_GNO|GK_EISE_CGW|dwil_GLEANR' \ -dump=match -bygene -igpart \ $em/dwil1/scaf*/*_exons.nr.blatf8 > & $em/dwil1/dwil1-tandynear2ipdump2.txt & =item more output * check dump gene id stats for quality of dupl. matches $td/tandynear2.perl -neardist=15000 -dump=matchstats -minali=0.1 \ -igpartial -nogroup -debug=1 \ scaf*/*NCBI_GNO-pexons.fa.mbout.gz > & ! dwil1-tandynear41mb2gno.txt & =item mark TE repeats ** apply to .blat/.blast output also ** gzcat $em/dwil1/scaffold*/*tandy6j.gff.gz | \ perl $td/repeatfilter.perl -act mark -r $dre/dwil_caf1_BLRtx_pilerTE.gff3.gz \ -gff stdin > dwil_caf060213_exons_tandy6j.gff set ggroup=dgri_GLEANR set ggroup=GH_NCBI_GNO GK_NCBI_GNO dwil_GLEANR GG_NCBI_GNO dere_GLEANR GM_NCBI_GNO dsec_GLEANR GA_NCBI_GNO dpse_GLEANR GE_NCBI_GNO dyak_GLEANR GF_NCBI_GNO dana_GLEANR CG_NCBI_GNO FBtr : dmel gzcat scaff*/$ggroup-pexons.fa.mbout.gz \ | $td/overlapfilter.perl -act mark -mark=terepeat -overlap d*_pilerTE.gff3.gz -itype blast -in stdin \ | $td/overlapfilter.perl -act mark -mark=dmelhsp -overlap d*-dmelhsp.gff.gz -itype blast -in stdin \ > $ggroup-pexons.fa.mbout gzcat scaff*/$ggroup-pexons.fa.gz > $ggroup-pexons.fa ; gzip *-pexons.fa* & $td/tandynear2.perl -neardist=15000 -dump=matchstats -minali=0.1 -igpartial -nogroup -debug=1 \ -overmarkfilter 'terepeat|notdmelhsp' \ $ggroup-pexons.fa.mbout.gz > & $ggroup-pexons.tandynear41mb.txt & # daphe; no pilerTE repeats; use mark=prot9hsp -overlap dpulex-prot9-hsp.gff.gz NCBI_GNO Dappu =item try with blat; megablast weak matches suspect (dmel high) GG_NCBI_GNO dere_GLEANR gzcat scaff*/dere_caf060210_exons.nr.blatf8.gz \ | grep $ggroup \ | $td/overlapfilter.perl -act mark -mark=terepeat -overlap d*_pilerTE.gff3.gz -itype blast -in stdin \ | $td/overlapfilter.perl -act mark -mark=dmelhsp -overlap d*-dmelhsp.gff.gz -itype blast -in stdin \ > $ggroup-pexons.fa.blat $td/tandynear2.perl -neardist=15000 -dump=matchstats -minali=0.1 -igpartial -nogroup -debug=1 \ -overmarkfilter 'terepeat|notdmelhsp' \ $ggroup-pexons.fa.blat.gz > & $ggroup-pexons.tandynear41blat.txt & =cut use strict; use warnings; use Getopt::Long; use constant SAMEBASE => 10; # for _sameloc, slop allowed in loca == locb our $BINSIZE = 1000 ; #was# 5000; # test with high identity: MINALIGN>0.9 our $MINEVAL = 1e-5; our $MINALIGN = 0.5; our $NEARDIST = 8000; # was 15k our $NEARSTEPS= 3; our $EQSAME=1; our $EQINSIDE=3; our $EQOVERLAP=2; #? ignore our $EQNEAR=-2; our $EQFAR=-999; ## replace w/ steps # my $EQNEAR1=-3; my $EQNEAR2=-4; my $EQFAR=-999; my $FILTER_GENEGROUP_OVERLAPS=1; # NOT USED NOW: does this want a cli option? my $overmarkfilter=""; #use constant COUNTKEYS => 0; my $COUNTKEYS = 0; #my @eqclass= ($EQSAME, $EQNEAR, $EQNEAR1, $EQNEAR2, $EQFAR, $EQINSIDE); # same, near, far #my @eqnames= ("Same","Near8k","Near15k","Near45k","Far","Inside"); # same, near, far my @eqclass= ($EQSAME, $EQINSIDE, $EQNEAR, $EQFAR); # same, near, far my @eqnames= ("Same","Inside","Near","Far"); # same, near, far my $dumpgenes=0; my $igpartialexon=0; my $groupok=1; my $showcount= 1; my $debug=0; #use constant BY_EXON => 0; #use constant BY_GENE => ! BY_EXON; my $BY_GENE= 0; # == count_by_geneid vs by_exon my (%hits, %ggroup, %altids, $gffin, %genespan, %genespanall, @fasta, $findoverlap, $skipids, $refgroup, $keeponly); my %didtypeloc=(); my $keep1exonscore=0; #? instead read both blast.tables and query.nr fasta from cmdline w/o switches # separate by ? suffix my $optok= GetOptions( # "fasta|queryfasta=s",\@fasta, "skipids=s", \$skipids, "keeponly=s", \$keeponly, # regex for id-group "MINALIGN=s", \$MINALIGN, "MINEVAL=s", \$MINEVAL, "NEARDIST=i", \$NEARDIST, "NEARSTEPS=i", \$NEARSTEPS, "overmarkfilter=s", \$overmarkfilter, "groupok!", \$groupok, "BYGENE!", \$BY_GENE, "COUNTKEYS!", \$COUNTKEYS, "refgroup=s", \$refgroup, "dumpgenes=s", \$dumpgenes, "showcount!", \$showcount, "findoverlap!", \$findoverlap, # uses MINALIGN "gffin!", \$gffin, "oneexon!", \$keep1exonscore, "igpartialexon!", \$igpartialexon, "debug=i", \$debug, ); die "usage: perl tandynear.perl [ options ] \ scaffold_6308/dmoj_caf060210_exons.nr.blatf8 # presumes fasta dmoj_caf060210_exons.nr opts: -[no]group; -skipids=[+]listfile; -keeponly=IDgroup -debug=0,1,2; -[no]ignorepartialexon; -[no]findoverlap; -neardist=$NEARDIST; -nearsteps=$NEARSTEPS; -minalign=$MINALIGN; -mineval=$MINEVAL; -overmarkfilter=terepeat|!dmelhsp .. remove/keep matches with overlap filter marks -dumpgenes=match|any|found + location|stats " unless($optok && @ARGV); #my @nears= map{ my $kb= int((1 + $NEARDIST * $_)/1000); "Near".$kb."k"; } (1..$NEARSTEPS); @eqclass= ($EQSAME, $EQINSIDE, (map{$EQNEAR * $_}(1..$NEARSTEPS)) , $EQFAR); # same, near, far @eqnames= ("Same","Inside", (map{ my $kb= int((1 + $NEARDIST * $_)/1000); "Near".$kb."k"; } (1..$NEARSTEPS)) ,"Far"); # same, near, far my %skipids=(); $skipids ||=""; my $keepids = ($skipids =~ s/^[\+\^]//) ? 1 : 0; if(-f $skipids) { open(F,$skipids); while(){chomp; s/\s.*$//; $skipids{$_}++ if(/^\w/);} close(F); } $keeponly="($keeponly)" if($keeponly && $keeponly =~ m/\|/ && $keeponly !~ m/\(/); # $FILTER_GENEGROUP_OVERLAPS= 0 unless($groupok); my $pctover= $MINALIGN; my $nskip=0; my $npoor=0; my $nexon1=0; my(%predfound, %predolap, %predcount, %predgroup, %predalign, %predmarks, %predrev, %didxeq); ## results over all scaffolds my(%predfoundover); %hits=(); # want this outside scaffold loop? # try 1 scaff at a time; problems otherwise? for( my $iarg= $#ARGV; $iarg >= 0; $iarg--) { my $outblat= $ARGV[$iarg]; warn "# processing[$iarg] $outblat\n" if $debug; %hits=(); # want this outside scaffold loop? %altids=(); %genespan=(); %didtypeloc=(); my %genealtids= (); my %refgrouplocs=(); my $fasta= $outblat; $fasta =~ s/\.gz$//; $fasta =~ s/\.\w+$//; exon_fasta_ids($fasta, \%altids, \%genespan, \%genealtids); map { $genespanall{$_} = $genespan{$_}; } keys %genespan; # save for dump; or dont erase genespan # is this genespan filter too restrictive? shouldnt chop blat matches by this # but matching gene ids .. i.e. all Gno genes overlapping in part GW genes # or GW gene region, or one of their exons if($refgroup) { foreach my $gid (sort keys %genespan) { my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); next unless($ggroup =~ m/$refgroup/); my ($tb,$te,$gref)= @{$genespan{$gid}}; $tb -= 2 * $BINSIZE; $te += 2 * $BINSIZE; # add some nearby slop my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { push @{$refgrouplocs{$ib}}, $gid; } } } my $ok=0; if(0) { # use tandy's sort, fix methods? is sort needed here? my $fixlocperl= '@v=split"\t"; ($b,$e)=@v[8,9]; if($b>$e){ @v[8,9]=($e,$b); $v[0].="-";} print join("\t",@v);'; my $sortblast8= 'sort -k2,2 -k9,9n -k10,10nr -k12,12nr'; # genome.ref,b,e ; bitscore $ok= open(BLAT, "cat $outblat | perl -ne '$fixlocperl' | $sortblast8 |") or die "sort $outblat"; } elsif($outblat=~/\.gz$/) { $ok= open(BLAT,"gunzip -c $outblat |") or die "open $outblat"; } else { $ok= open(BLAT,$outblat) or die "open $outblat"; } my ($overmarkdrop,$overmarkkeep)=("",""); if($overmarkfilter) { $overmarkdrop= $overmarkfilter; while($overmarkdrop =~ s/\|?(\!|not)(\w+)//){$overmarkkeep .= $2."\|";} $overmarkkeep=~s/\|$//; $overmarkdrop=~s/\|$//; } while() { next unless(/^\w/); chomp; my @v=split"\t"; # my($qid1,$ref,$pctid,$alen,$tb,$te,$eval,$bits)=@v[0,1,2,3,8,9,10,11]; # add gff input option? my($qid1,$ref,$pctid,$alen,$tb,$te,$tor,$eval,$bits,$attr); if($gffin) { ($ref,$tb,$te,$eval,$tor,$attr) = @v[0,3,4,5,6,8]; #($qid1) = $attr =~ m/(?:ID|Parent)=([^;]+)/; if($attr =~ m/ID=([^;]+)/) {$qid1=$1; } elsif($attr =~ m/Parent=([^;]+)/) { $qid1=$1; } $eval=0 if($eval eq "."); $tor="" if($tor eq "."); $pctid= $alen= $bits= 0; } else { ($qid1,$ref,$pctid,$alen,$tb,$te,$eval,$bits) = @v[0,1,2,3,8,9,10,11]; } ## add overlapfilter marks my $overmarks= "" . join ",",@v[12..$#v]; # my $trev= ($tor eq "-") ? -1 : 0; if($tb>$te) { $tor="-"; ($tb,$te)= ($te,$tb); } if($qid1 =~ s/-$//) { $tor="-"; } $npoor++ and next if($eval > $MINEVAL); $npoor++ and next if($overmarks and $overmarkdrop and $overmarks =~ m/$overmarkdrop/); $npoor++ and next if($overmarkkeep and $overmarks !~ m/$overmarkkeep/); # filter by palign here also # FIXME: MINALIGN is removing valid exon matches broken into HSPs # e.g. 1-exon gene matches nearby over 90% of size, but in 7 HSPs # need to look at all close exon HSPs, qb,qe range as per prothsp mapping my($qid11,$qb1,$qe1,$qo1)=split(/[:-]/,$qid1); # exon id,start,stop if($MINALIGN>0 and $alen>0 and $qe1>0) { my $elen = 1+abs($qe1-$qb1); my $palign= $alen / $elen; $npoor++ and next if($palign < $MINALIGN); } my $notinrefgroup= 0; if($refgroup) { ## && $ggroup =~ m/$refgroup/ my @sabins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); my @samelocs=(); foreach my $ib (@sabins) { push @samelocs, @ { $refgrouplocs{$ib} } if($refgrouplocs{$ib}); } my $found=0; my $olap=0; foreach my $sgid (@samelocs) { $genespan{$sgid} or next; my ($gb,$ge,$gref)= @{$genespan{$sgid}}; if( $gref eq $ref && _isoverlap($gb,$ge,$tb,$te) ) { # add _isnear() test? $olap= 1; last; } } # defer to each altid# next unless($olap||$found); $notinrefgroup=1 unless($olap); } my $altids= $altids{$qid1} || []; ## do we really want to count same match so many times? ## once per ggroup?, but need to know near/far matches along w/ same per ggroup %didxeq=(); # this now filters per ggroup,match,eq type my %didqid; # any dups in altids? foreach my $qid_be ($qid1, @$altids) { next if($keeponly and $qid_be !~ m/$keeponly/); my($qid,$qb,$qe,$qor)=split(/[:-]/,$qid_be); # exon id,start,stop $nskip++ and next unless(defined $qe); if($qb>$qe) { ($qb,$qe)=($qe,$qb); $qor="-"; } my($gid) = $qid =~ m/^(\w+)/; # using ID.exnum syntax my($ggroup)= ($groupok) ? $qid =~ m/^(\D+)/ : ("all"); # BAD for non-grouping-prefixed IDs ## need to match gene-group of each altgn to refgroup; stay if matched if($refgroup && $notinrefgroup) { my @altgn = sort keys %{ $genealtids{$gid} }; foreach my $agn (@altgn) { if($agn =~ m/$refgroup/) { $notinrefgroup= 0; last; } } next if($notinrefgroup); } # invert this also: keep only those in %skipids if(%skipids){ if($keepids) { $nskip++ and next unless($skipids{$gid} || $skipids{$qid}); } else { $nskip++ and next if($skipids{$gid} || $skipids{$qid}); # gene/transcript id here; or both? (exonid) } } # save exon match locs; after read all, recount if near/far has other prediction at sameloc my ($gb,$ge,$gref,$gid2,$gnexons,$gcdslen); if(scalar(%genespan)) { $nskip++ and next unless ($genespan{$gid}); ($gb,$ge,$gref,$gid2,$gnexons,$gcdslen) = $genespan{$gid} ? @{$genespan{$gid}} : (0) x 10; # genespan{gid} == [$gb,$ge,$gref,$gid2,$gnexons,$gcdslen] } else { ($gb,$ge,$gref,$gid2,$gnexons,$gcdslen)= ($qb,$qe,$ref,$gid,1,$qe-$qb+1); } my $eq= -1; my $isrev= 0; if( _sameloc($ref,$qb,$qe,0, $ref,$tb,$te,0) ) { $eq= $EQSAME; } elsif ($ge>0 && $gref eq $ref && _isoverlap($gb,$ge,$tb,$te)) { # gene span $eq= $EQINSIDE; } else { my $dist= _mindistance2($qb,$qe,$tb,$te); # _isnear2($qb,$qe,$tb,$te); $isrev= 0; if($qor && $tor) { $isrev = ($qor ne $tor) ? 1 : 0; } $eq= $EQFAR; for (my $iq=1; $iq<=$NEARSTEPS; $iq++) { if($dist <= $iq * $NEARDIST) { $eq= $iq * $EQNEAR; last; } } } #?? check for overabundant: repeats? keep only 1 score per exonid/eq class? #?? next if ($keep1exonscore and $didxeq{$qid_be.$eq}++); ## per-match line, check for ggroup.eq $nexon1++ and next if($keep1exonscore and $didxeq{$ggroup.$eq}++); if($debug>2) { (my $did=$qid) =~ s/\.\d+//; my $deq=($eq>0)?1:($eq<=-3)?0:-2; print "$did\t$deq\n"; } my $loc= [$tb,$te,$qid_be,$ref,$eq, $overmarks, $isrev]; ## tb,te genome loc here, not qb,qe my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE)); foreach my $ib (@bins) { push @{$hits{$ggroup}{$eq}{$ib}}, $loc; } $ggroup{$ggroup}{$eq}++; #? can we count geneids $gid instead of exons here? $ggroup{$ggroup}{total}++; } } close(BLAT); collectGroups(); ## %ggroup, %hits > (%predfound, %predolap, %predcount) } warn "# poor ids=$npoor; skipped ids=$nskip; only1-exon-eq=$nexon1\n" if($debug && $nexon1+$npoor+$nskip>0); exit if($debug>2); # print_byexon() if BY_EXON; print_bygene(); # if BY_GENE; dump_genes() if $dumpgenes; #-------------------------------------------------------------------------- # only for near classes ... found, anyfound? sub dump_genes { my @ggroup= sort keys %ggroup; my $manygroup= @ggroup>1; my @eqnear= map{$EQNEAR * $_}(1..$NEARSTEPS); my @nanear= @eqnames[2..$NEARSTEPS+1]; # want also ids from predcount, predolap ? my $matchclass; my $predlist; my $doloc = ($dumpgenes =~ /loc/); my $dostat = ($dumpgenes =~ /stat/); if($dumpgenes =~ /match/) { $matchclass= "matched"; # predcount $predlist= \%predcount; } elsif($dumpgenes =~ /any/) { $matchclass= "any_found"; $predlist= \%predfound; # and predolap, below } else { $matchclass= "found"; # predfound ? include predolap/anyfound here? $predlist= \%predfound; } my @tablehead= qw(group gid loc near found anyfnd nxfound nexon align cdslen marks); # max1found print "\n"; print "# id quality table for near-duplicate genes, class=$matchclass\n"; print join("\t",@tablehead),"\n" if($dostat); foreach my $ggroup (@ggroup) { my $eqi=0; ## urk, eqnear steps are including many of same genes; separate out dupls? foreach my $eq (@eqnear) { my @idlist= sort keys %{$predlist->{$ggroup}{$eq}}; if($matchclass =~ /any/){ my @idother= keys %{$predolap{$ggroup}{$eq}}; my %idboth= map{ $_,1} @idlist, @idother; @idlist= sort keys %idboth; } my %genestat; my $ngene= @idlist; my $neara= $nanear[$eqi]; # near_${eq} print "# group: $ggroup; $neara=$ngene\n"; ## if matched, mark ones that are found foreach my $id (@idlist) { my $mark=""; if($predfound{$ggroup}{$eq}{$id}) { # $matchclass =~ /match/ and $mark="\tfound"; } elsif($predolap{$ggroup}{$eq}{$id}) { # $matchclass =~ /match/ and $mark="\tanyfnd"; } my $loc=""; if($doloc && $genespanall{$id}) { my ($gb,$ge,$gref)= @{$genespanall{$id}}; # genespan{gid} == [$gb,$ge,$gref,$gid2,$gnexons,$gcdslen] $loc="\t$gref\t$gb\t$ge" if($gref); } ## for not BY_GENE, collect exons/gene, show counts of matched/found exons, total gene exons ## ? align total / gene cdslen ## add marks to exons.fa for dmel/etc-homology, terepeat, ... as per tandy.gff if($dostat) { # if dogenestat; assumes not BY_GENE ? my ($gid,$xn,$xb,$xe)= split(/[\.:-]/,$id); # geneid,exnum,start,stop $xn ||= 1; my ($gb,$ge,$gref,$gid2,$gnexons,$gcdslen)= @{$genespanall{$gid}}; my $palign= $predalign{$ggroup}{$eq}{$id} || 0; # sum for all hsps/eq/id my $predmarks= $predmarks{$ggroup}{$eq}{$id} || 0; # sum for all hsps/eq/id $genestat{$gid}{loc } = "$gref:$gb-$ge"; # gene stats $genestat{$gid}{nexon} = $gnexons; # gene stats $genestat{$gid}{cdslen} = $gcdslen; $genestat{$gid}{align} += $palign; # align/xn so many-dup-1exon doesn't swamp ? $genestat{$gid}{exonsfound}{$xn}++; my $xa= $genestat{$gid}{exonsalign}{$xn} || 0; $genestat{$gid}{exonsalign}{$xn} = $palign if($xa < $palign); $genestat{$gid}{found} ++ if($mark =~ /found/); $genestat{$gid}{anyfnd} ++ if($mark =~ /anyfnd/); if($predmarks) { my @mk= split",",$predmarks; my %mk=(); map{ my($k,$v)=split"=",$_; $mk{$k} += $v||1; } @mk; $predmarks= join ",", map{ $_ ."=". $mk{$_}} sort keys %mk; $genestat{$gid}{marks} = $predmarks ; } } print "$id$mark$loc\n" unless($dostat); } if($dostat && $ngene>0) { foreach my $gid (sort keys %genestat) { my($loc, $nexon,$cdslen,$align,$found,$anyfnd, $marks)= map{ $genestat{$gid}{$_}||0 } qw(loc nexon cdslen align found anyfnd marks); # my $nxfound= scalar(keys %{$genestat{$gid}{exonsfound}} ); # ^ add min/max nx count value? my $gsref= $genestat{$gid}{exonsfound} ; my @nxfound= sort{$gsref->{$b}<=>$gsref->{$a}} keys %$gsref; my $nxfound= @nxfound; my $ti=($nxfound>1) ? 1 : 0; # my $txfound= $gsref->{$nxfound[$ti]}; my $xalign=0; map{ $xalign+=$_; } values %{$genestat{$gid}{exonsalign}}; my($gidg,$gidn)= $gid =~ m/^(\D*)(.*)/; print join("\t",$gidg,$gidn,$loc,$neara, $found,$anyfnd,$nxfound,$nexon,$xalign,$cdslen,$marks),"\n"; # $txfound, } } $eqi++; } print "\n"; } } sub print_bygene { my @ggroup= sort keys %ggroup; my $manygroup= @ggroup>1; my $ttype = ($BY_GENE) ? "gene" : "exon"; $ttype .= ($COUNTKEYS) ? " ids" : " matches"; my $cc = $dumpgenes ? "#t " : ""; my @flags=(); push(@flags,reference_group => $refgroup) if($refgroup); push(@flags,count_one_exon => $keep1exonscore) if($keep1exonscore); push(@flags,found_is_overlap => $findoverlap) if($findoverlap); push(@flags,ignore_partial_match => $igpartialexon) if($igpartialexon); push(@flags,keep_only => $keeponly) if($keeponly); push(@flags,overlap_filter => $overmarkfilter) if($overmarkfilter); push(@flags,by_gene => $BY_GENE); push(@flags,count_type => (($BY_GENE)? "gene" : "exon")); push(@flags,count_by => (($COUNTKEYS)? "ids" : "matches")); push(@flags,min_align => $MINALIGN); push(@flags,min_e_value => $MINEVAL); my %flags=@flags; my $flags= join ",", map{"$_=".$flags{$_}} sort keys %flags; print $cc,"Tandy count of $ttype per predictor group\n"; print $cc,"Options: $flags\n"; print $cc,join("\t","Group ", "Stat",@eqnames),"\n"; foreach my $ggroup (@ggroup) { #?? maybe counts here should be sum of values of hash, not count of keys my ($pc,$pv,$tsame); if($COUNTKEYS) { $tsame = scalar( keys %{$predcount{$ggroup}{$EQSAME} } ) || 1; } else { $tsame=0; map{ $tsame += $_ } values %{$predcount{$ggroup}{$EQSAME}}; } my @counts= map{ if($COUNTKEYS) { $pc= scalar( keys %{$predcount{$ggroup}{$_}} ) || 0; } else { $pc=0; map{ $pc += $_ } values %{$predcount{$ggroup}{$_}}; } $pc; } @eqclass; my @freqs= map{ sprintf("%.3f", $_ / $tsame ); } @counts; my $sgroup= sprintf("%-15s", substr($ggroup,0,14)); # my @isrev= map{ # if($COUNTKEYS) { # $pc= scalar( keys %{$predrev{$ggroup}{$_}} ) || 0; # } else { # $pc=0; map{ $pc += $_ } values %{$predrev{$ggroup}{$_}}; # } # $pc; # } @eqclass; my @foundfreq= map{ if($COUNTKEYS) { $pc= scalar( keys %{$predcount{$ggroup}{$_}} ) || 1; $pv= scalar( keys %{$predfound{$ggroup}{$_}} ) ; } else { $pc=0; map{ $pc += $_ } values %{$predcount{$ggroup}{$_}}; $pv=0; map{ $pv += $_ } values %{$predfound{$ggroup}{$_}}; } $pc = ($showcount || $pc<1) ? 1 : $pc; my $fmt= ($showcount)? "%d" : "%.3f"; sprintf($fmt, $pv / $pc); } @eqclass; my @anyfnd= map{ if($COUNTKEYS) { $pc= scalar( keys %{$predcount{$ggroup}{$_}} ) || 1; $pv= scalar( keys %{$predolap{$ggroup}{$_}} ) ; } else { $pc=0; map{ $pc += $_ } values %{$predcount{$ggroup}{$_}}; $pv=0; map{ $pv += $_ } values %{$predolap{$ggroup}{$_}}; } $pc = ($showcount || $pc<1) ? 1 : $pc; my $fmt= ($showcount)? "%d" : "%.3f"; sprintf($fmt, $pv / $pc); } @eqclass; my @foundofany= map{ if($COUNTKEYS) { $pc= scalar( keys %{$predolap{$ggroup}{$_}} ) || 1; $pv= scalar( keys %{$predfound{$ggroup}{$_}} ) ; } else { $pc=0; map{ $pc += $_ } values %{$predolap{$ggroup}{$_}}; $pc ||=1; $pv=0; map{ $pv += $_ } values %{$predfound{$ggroup}{$_}}; } $pc ||=1; sprintf("%.3f", $pv / $pc); } @eqclass; print $cc,join("\t",$sgroup,"count",@counts),"\n"; # $total print $cc,join("\t",$sgroup,"freq",@freqs),"\n"; # print $cc,join("\t",$sgroup,"rever",@isrev),"\n"; # not useful print $cc,join("\t",$sgroup,"found",@foundfreq),"\n"; print $cc,join("\t",$sgroup,"anyfnd",@anyfnd),"\n" if $manygroup; print $cc,join("\t",$sgroup,"fnd/any",@foundofany),"\n" if $manygroup; print "\n"; } } # sub print_byexon { # # my @ggroup= sort keys %ggroup; # my $manygroup= @ggroup>1; # print "Tandy exon match types per predictor group\n"; # print join("\t","Group ", "Stat",@eqnames),"\n"; # foreach my $ggroup (@ggroup) { # my $total = $ggroup{$ggroup}{total}; # my $tsame = $ggroup{$ggroup}{$EQSAME} || 1; # # # my @counts= map{ $ggroup{$ggroup}{$_} || 0 } @eqclass; #? or $predcount{$ggroup}{$_} # my @counts= map{ $predcount{$ggroup}{$_} || 0 } @eqclass; # my @freqs= map{ sprintf("%.3f", $_ / $tsame ); } @counts; # my $sgroup= sprintf("%-15s",$ggroup); # # my @foundfreq= map{ # my $pc = ($showcount)? 1: ($predcount{$ggroup}{$_} || 1); # my $fmt= ($showcount)? "%d" : "%.3f"; # sprintf($fmt, $predfound{$ggroup}{$_} / $pc); # } @eqclass; # # my @anyfnd= map{ # my $pc = ($showcount)? 1: ($predcount{$ggroup}{$_} || 1); # my $fmt= ($showcount)? "%d" : "%.3f"; # sprintf($fmt, $predolap{$ggroup}{$_} / $pc); # } @eqclass; # # my @foundofany= map{ # my $pany= $predolap{$ggroup}{$_} || 1; # sprintf("%.3f", $predfound{$ggroup}{$_} / $pany); # } @eqclass; # # my @ftotal= map{ # $predcount{$ggroup}{$_} || 1; # } @eqclass; # # # if manygroup, compute/print foundfreq / anyfound freq # # print join("\t",$sgroup,"count",@counts),"\n"; # $total # print join("\t",$sgroup,"freq",@freqs),"\n"; # print join("\t",$sgroup,"found",@foundfreq),"\n"; # print join("\t",$sgroup,"anyfnd",@anyfnd),"\n" if $manygroup; # print join("\t",$sgroup,"fnd/any",@foundofany),"\n" if $manygroup; # print "\n"; ## if($sgroup ne $lastgroup); $lastgroup= $sgroup; # # anyfound almost same as found # # overlap almost same as found # # print join("\t",$sgroup,"ftotal",@ftotal),"\n"; # same as @counts # } # # } #-------------------- sub collectGroups { my @ggroup= sort keys %ggroup; my $manygroup= @ggroup>1; foreach my $ggroup (@ggroup) { # if(BY_EXON) { # $predcount{$ggroup}{$EQSAME} ||= 1; # $predolap{$ggroup}{$EQSAME} ||= 1; # $predfound{$ggroup}{$EQSAME} ||= 1; # } my @eqlist= @eqclass; ## NOT NOW: shift @eqlist; # drop EQSAME ? foreach my $eq (@eqlist) { # if(BY_EXON) { # $predcount{$ggroup}{$eq} ||= 0; # $predfound{$ggroup}{$eq} ||= 0; # $predolap{$ggroup}{$eq} ||= 0; # } my (@duplocs,%diddup); my @bins= sort keys %{ $hits{$ggroup}{$eq} }; foreach my $ib (@bins) { # lots o these; all in scaffold here ?? want that? push @duplocs, @ { $hits{$ggroup}{$eq}{$ib} }; } foreach my $dp (@duplocs) { # match exons from dupl class and same? my($dqb,$dqe,$dqid_be,$dref,$deq, $overmarks, $disrev)= @$dp; #^^ == [$tb,$te,$qid_be,$ref,$eq, $overmarks]; ##next if($diddup{$dp}++); next if($diddup{$dqb.$dqe.$dqid_be.$dref.$deq}++); next unless($deq == $eq); # redundant my $id= $dqid_be; $id =~ s/\W.*$// if($BY_GENE); $predalign{$ggroup}{$eq}{$id} += 1 + abs($dqe-$dqb); #? push list? $predmarks{$ggroup}{$eq}{$id} .= $overmarks."," if($overmarks); if ($eq == $EQSAME) { $predfound{$ggroup}{$EQSAME}{$id}++; # if $found; # always found same $predolap{$ggroup}{$EQSAME}{$id}++; # if ($found || $otherfound); $predcount{$ggroup}{$EQSAME}{$id}++; # $predrev{$ggroup}{$EQSAME}{$id}++ if($disrev); #??? next; } my @sabins= (int($dqb/$BINSIZE) .. int($dqe/$BINSIZE)); my @samelocs=(); foreach my $ib (@sabins) { push @samelocs, @ { $hits{$ggroup}{$EQSAME}{$ib} } if($hits{$ggroup}{$EQSAME}{$ib}); } my $found=0; my $olap=0; my $isrev= 0; foreach my $sp (@samelocs) { my($sqb,$sqe,$sqid0,$sref,$seq,$sxmark,$sisrev)= @$sp; next unless($sref eq $dref && $seq == $EQSAME); if( _sameloc($sref,$sqb,$sqe,0, $dref,$dqb,$dqe,0) ) { $found= 1; last; } elsif( _isoverlap_pct($sqb,$sqe,$dqb,$dqe) ) { $olap= 1; } # $isrev=1 if($sisrev != $disrev); #?? what here } $isrev= $disrev; my $otherfound=0; if($manygroup) { my @anysame=(); foreach my $gg (@ggroup) { next if($gg eq $ggroup); # dont need to retest foreach my $ib (@sabins) { push @anysame, @ { $hits{$gg}{$EQSAME}{$ib} } if($hits{$gg}{$EQSAME}{$ib}); } } foreach my $sp (@anysame) { my($sqb,$sqe,$sqid0,$sref,$seq)= @$sp; next unless($sref eq $dref && $seq == $EQSAME); if( _sameloc($sref,$sqb,$sqe,0, $dref,$dqb,$dqe,0) ) { $otherfound= 1; last; } } } $found ||= $olap if($findoverlap); $predfound{$ggroup}{$eq}{$id}++ if $found; ## $predfoundover{$ggroup}{$eq}{$id}++ if ($olap || $found); # replace predfound option? # $predrev{$ggroup}{$eq}{$id}++ if $isrev; $predolap{$ggroup}{$eq}{$id}++ if ($found || $otherfound); $predcount{$ggroup}{$eq}{$id}++; } } } } sub min { return ($_[1] < $_[0] ? $_[1] : $_[0]); } sub max { return ($_[1] > $_[0] ? $_[1] : $_[0]); } sub _isoverlap { my($gb,$ge, $qb,$qe)= @_; return ($gb <= $qe && $ge >= $qb) ? 1 : 0; } sub _isoverlap_pct { my($tb,$te, $lb,$le)= @_; 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; my $pover= $maxo/$leno; $over = 0 if $pover < $pctover; } return $over; } sub _troverlap { my($gb,$ge,$gref, $trlist)= @_; foreach my $tr (@$trlist) { #[$gb,$ge,$gref,$gid]= $tr; my $overlap = $gref eq $$tr[2] && _isoverlap($gb,$ge, $$tr[0],$$tr[1]); return $tr if $overlap; } } sub _mina { my $a= shift; my $b= shift; return $a unless(defined $b); return _mina( (abs($a) < abs($b) ? $a : $b), @_); } sub _isnear { my($gb,$ge, $qb,$qe)= @_; #? change to test not midpoint but smallest dist ? Has no sig. effect # my $dmin= _mina($gb - $qe, $ge - $qb, $gb - $qb, $ge - $qe ); # return (abs($dmin) < $NEARDIST) ? 1 : 0; my $gm= int(($gb+$ge)/2); my $qm= int(($qb+$qe)/2); return (abs($gm - $qm) < $NEARDIST) ? 1 : 0; } sub _isnear2 { my($gb,$ge, $qb,$qe)= @_; my $gm= int(($gb+$ge)/2); my $qm= int(($qb+$qe)/2); return abs($gm - $qm); } sub _mindistance2 { my($qb,$qe, $tb,$te)= @_; my $bd= abs($tb - $qe); # g above q my $ed= abs($qb - $te); # g below q ; dont need tb - qb, te - qe test return ($ed < $bd) ? $ed : $bd; } # 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 # # if ($tb <= $qe && $te >= $qb) { return (0,$rev); } # or -? # # #? what of _isinside ? # 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 _sameloc { my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_; # $aor ||= 0; $bor ||= 0; # $ab ||=0; $ae ||=0; $bb ||=0; $be ||=0; #? need return ((abs($ab-$bb) <= SAMEBASE) && (abs($ae-$be) <= SAMEBASE)) ? 1 : 0 ; } sub exon_id_crossref { my($altids, $genealt)= @_; ## prnote "exon_id_crossref for $queryfa\n" if $debug; # 1. crossref the exonids from alt lists my %xaltid=(); my $nxin= 0; foreach my $xid (sort keys %{$altids}) { $xaltid{$xid}{$xid}=1; $nxin++; my($gid) = $xid =~ m/^(\w+)/; $genealt->{$gid}{$gid}= 1; my $alist= $altids->{$xid} || []; ## or next; foreach my $axid (@$alist) { $xaltid{$xid}{$axid}= 1; $xaltid{$axid}{$xid}= 1; $xaltid{$axid}{$axid}= 1; # ensure self match # ditto for gene ids my($agid) = $axid =~ m/^(\w+)/; $genealt->{$gid}{$agid}= $genealt->{$agid}{$gid}= $genealt->{$agid}{$agid}= 1; } } # 1b. update alt lists with all xrefs; must keep xloc as part of xid my $nxout= 0; foreach my $xid (sort keys %xaltid) { my @blist = sort keys %{$xaltid{$xid}}; $nxout++; if(@blist) { my $alist = $altids->{$xid}||[]; my %ablist= map{ $_,1 } (@$alist, @blist); my @ablist= sort keys %ablist; $altids->{$xid}= \@ablist; } } return $altids; } sub exon_fasta_ids { my($queryfa, $altids, $genespan, $genealts)= @_; my(%gloc, %didexongloc, %alttr, %isalttr); my($nexons, $naltexons, $ntrans, $ngenes, $atref, $atstrand, $isalttr)=(0)x10; my ($ref,$strand,$loc)=(0,0,"",""); # my $newaltids= {}; $altids ||= {}; $genespan ||= {}; $genealts ||= {}; $queryfa.=".gz" unless( -f $queryfa); warn "#exon_fasta_ids Error: $queryfa not found\n" and return unless( -f $queryfa); my $ok=0; if($queryfa =~ /.gz$/) { $ok= open(GREP, "gunzip -c $queryfa | grep '^>' |") ; } else { $ok= open(GREP, "grep '^>' $queryfa|"); } return ($altids, $genespan, $genealts) unless($ok); while(){ my($xid)=m/>(\S+)/ or next; my($aid); if(m/altids=(\S+)/){$aid=$1;} ($ref,$strand,$loc)=(0,0,"",""); if(m/loc=([^;\s]+)/) { ($ref,$loc,$strand)= split(/:/,$1); ($atref, $atstrand)= ($ref,$strand); } if($keeponly && $xid !~ m/$keeponly/) { if($aid && $aid =~ m/^($keeponly\w+\.\d+:\d+-\d+)/) { $xid= $&; ## $1; } else { next; #? } } my($gid,$xn,$xb,$xe) = $xid =~ m/^(\w+)\.(\d+):(\d+)-(\d+)/; unless($gid && $xe) { ($gid,$xb,$xe) = $xid =~ m/^(\w+):(\d+)-(\d+)/; $xn=1; } warn "bad xid: $gid.$xn:$xb-$xe = $xid\n" unless($gid && $xe); # for no .xn exon num -- older bug ?? # GI_EISE_CEX_11198402:577674-577795 #??? unless($didexongloc{$xid}) { } $gloc{$gid} or $gloc{$gid}=[]; push @{$gloc{$gid}}, [$xn,$ref,$xb,$xe,$strand]; $didexongloc{$xid}++; $nexons++; if($aid) { $aid =~ s/pm,.*$// if ($igpartialexon); # no partials here ?? # my @altids= split(",",$aid); # my @altids= grep{ $xid ne $_ } split(",", $aid); #? exclude $xid here or not # ## check that each has ID.exnum:locb-loce syntax ? some dont # ## m/^(\w+)\.(\d+):(\d+)-(\d+)/; my @altids= grep{ # $xid ne $_ and if($keeponly && $_ !~ /$keeponly/){ $_="skip"; } (m/^\w+\.\d+:\d+\-\d+/) } split(",", $aid); #? exclude $xid here or not $altids->{$xid}= \@altids; } } close(GREP); ## urk; this is done for each input scaffold on same $altids ref; ## at the least, this is inefficient; avoid it ## FIXME: need to crossref all altids, otherwise get biased distribution of groups ## ?? do we need to update %gloc{$agid} with new crossref ids? exon_id_crossref($altids, $genealts); my $missinstrand=0; foreach my $xid (sort keys %{$altids}) { # NOT for all ref/scaffold, just input one next if($didexongloc{$xid}); my $alist= $altids->{$xid} or next; foreach my $ad (@$alist) { $naltexons++; my($agid,$axn,$axb,$axe) = $ad =~ m/^(\w+)\.(\d+):(\d+)-(\d+)/; unless($agid && $axe) { ($agid,$axb,$axe) = $xid =~ m/^(\w+):(\d+)-(\d+)/; $axn=1; } next unless($agid); # where this error? # # double check overlap alt-trs ? # # for GeneWise and like that call many transcripts/exons at same location # # if FILTER_GENEGROUP_OVERLAPS # my $isdup= $alttr{$agid} || # filter_genegroup_overlaps( $ref, $agid, $axn, $axb, $axe); # if($isdup) { # my($ftype)= ($groupok) ? $agid =~ m/^(\D+)/ : ("all"); # $isalttr++; $isalttr{$ftype}++; $alttr{$agid}++; # # delete $gloc{$gid}; #?? # #? delete $altids->{$xid}; #? this is right? alt-tr is alt-tr everywhere? # # next; # } $gloc{$agid} or $gloc{$agid}=[]; push @{$gloc{$agid}}, [$axn,$atref,$axb,$axe,$missinstrand]; $didexongloc{$ad}++; } } # # now copy this set into main altids # foreach my $xid (sort keys %{$newaltids}) { # NOT for all ref/scaffold, just input one # $altids->{$xid}= $newaltids->{$xid}; # } ## not good enough yet; input IDs are transcript ids ## need to reloop and look for tr overlaps to call gene span ## but this is merging alternate predictor groups into 1 gene, do we want that? my $nerr=0; # my @genelist; my %genelist; # per ggroup my $fref="unknown"; foreach my $gid (sort keys %gloc) { my($gb,$ge,$gref,$gstrand)=(-1,0,0,0); my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); my ($gnexons, $gcdslen)=(0,0); foreach my $ex (@{$gloc{$gid}}) { my($xn,$xref,$xb,$xe,$xstrand)= @$ex; unless($xe) { $nerr++; next;} # error ?? $gb=$xb if($gb==-1 || $xb<$gb); $ge=$xe if($xe>$ge); $gref=$xref; $gstrand= $xstrand; $gcdslen += 1+abs($xe-$xb); $gnexons++; } # ... merge genespans of overlapping transcripts my ($gene,$trid); my $genelist= $genelist{$ggroup} || []; if($gene= _troverlap( $gb,$ge,$gref, $genelist) ) { $gb= min($gb, $$gene[0]); $ge= max($ge, $$gene[1]); $$gene[0]= $gb; $$gene[1]= $ge; # $trid= $$gene[3]; } $ntrans++; $fref= $gref; my $trloc= [$gb,$ge,$gref,$gid,$gnexons,$gcdslen]; # add nexons, cdslen for gene-align test $genespan->{$gid}= $trloc; # $genespan->{$trid}= $trloc if($trid); # redundant w/ above push( @{$genelist{$ggroup}}, $trloc) unless($gene); } my $ngroup= scalar(keys %genelist); map{ $ngenes += @ {$genelist{$_}} } keys %genelist; warn "# exon_fasta_ids ref=$fref ngroup=$ngroup; ngenes=$ngenes; ntrans=$ntrans; nexons=$nexons; naltexons=$naltexons\n" if $debug>1; warn "# genespan errs=$nerr\n" if($debug && $nerr>0); if($debug && $isalttr>0) { my $atrs= join ",", map{ $_."=".$isalttr{$_};} sort keys %isalttr; warn "# alt-tr genes skipped=$atrs\n" ; } return ($altids, $genespan, $genealts); } sub filter_genegroup_overlaps { my($ref, $gid,$xn,$ab,$ae)= @_; # global: %didtypeloc; reset per scaffold return 0 unless($FILTER_GENEGROUP_OVERLAPS); $ref ||="thisscaffold"; my $loc= "$ab-$ae"; my($ftype)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); my $isolap=0; my @bins= (int($ab/$BINSIZE) .. int($ae/$BINSIZE)); LBIN: foreach my $ib (@bins) { $didtypeloc{$ftype.$ref}{$ib} or next; my @locs= @{$didtypeloc{$ftype.$ref}{$ib}}; foreach my $lc (@locs) { my ($lb,$le)= split "-",$lc; if(_isoverlap($ab,$ae,$lb,$le)) { $isolap=1; last LBIN; } } } return 1 if ($isolap); foreach my $ib (@bins) { push @{$didtypeloc{$ftype.$ref}{$ib}}, $loc; } return 0; } __END__ =item data prep gzcat scaff*/*NCBI*fa.mbout.gz | \ $td/overlapfilter.perl -act mark -mark=terepeat -overlap dgri_caf1_BLRtx_pilerTE.gff3.gz -itype blast -in stdin | \ $td/overlapfilter.perl -act mark -mark=dmelhsp -overlap dgri-dmelhsp.gff.gz -itype blast -in stdin \ > ! GH_NCBI_GNO-pexons.fa.mbout =item R stats tn1 <- read.table("dspp/dere1-tandynear41mb2gno.txt", header=T, strip.white=T) tn1.title<-"dere/gnomon/tandynear4" tn2 <- read.table("dspp/dere1-tandynear41mb2glr.txt", header=T, strip.white=T) tn2.title<-"dere/gleanr/tandynear4" tn3 <- read.table("dspp/dere1-tandynear41mb2cgw.txt", header=T, strip.white=T) tn3.title<-"dere/genewise/tandynear4" tw1 <- read.table("dspp/dwil1-tandynear41mb2gno.txt", header=T, strip.white=T) tw1.title<-"dwil/gnomon/tandynear4" tw2 <- read.table("dspp/dwil1-tandynear41mb2glr.txt", header=T, strip.white=T) tw2.title<-"dwil/gleanr/tandynear4" tw3 <- read.table("dspp/dwil1-tandynear41mb2cgw.txt", header=T, strip.white=T) tw3.title<-"dwil/genewise/tandynear4" paste(title); sum(tn$found>0); sum(tn$found==0); mean(tn$nxfound/tn$nexon); mean(tn$align/tn$cdslen) [1] "dere/gnomon/tandynear4" [1] 358 [1] 1291 [1] 0.6399089 [1] 0.7336216 #............... > tn1[1:10,] group gid found anyfnd nxfound nexon align cdslen 1 GK_NCBI_GNO_ 32000882 0 0 1 2 122 123 2 GK_NCBI_GNO_ 32001406 0 0 1 2 462 3090 3 GK_NCBI_GNO_ 32001668 2 0 2 2 822 822 4 GK_NCBI_GNO_ 32001916 0 0 2 5 128 2427 5 GK_NCBI_GNO_ 32003641 0 0 1 5 165 2100 6 GK_NCBI_GNO_ 32003903 0 0 1 8 423 2279 7 GK_NCBI_GNO_ 32003957 0 0 1 1 140 813 8 GK_NCBI_GNO_ 32004481 0 0 1 1 3515 789 9 GK_NCBI_GNO_ 32005005 0 0 1 1 3812 780 10 GK_NCBI_GNO_ 32005074 0 0 2 2 258 1878 =item more results w/ more near steps dghome2% perl $td/tandynear2.perl -debug=1 -neardist=5000 -nearstep=10 scaffold_?/dpulex1_exons.nr.blatf8 # exon_fasta_ids ref=scaffold_9 ngroup=3; ngenes=1049; ntrans=1126; nexons=7177; naltexons=7322 # exon_fasta_ids ref=scaffold_8 ngroup=3; ngenes=950; ntrans=1068; nexons=8122; naltexons=8595 # exon_fasta_ids ref=scaffold_7 ngroup=3; ngenes=924; ntrans=1041; nexons=7699; naltexons=7817 # exon_fasta_ids ref=scaffold_6 ngroup=3; ngenes=1461; ntrans=1540; nexons=8768; naltexons=11204 # exon_fasta_ids ref=scaffold_5 ngroup=3; ngenes=1034; ntrans=1157; nexons=8768; naltexons=9437 # exon_fasta_ids ref=scaffold_4 ngroup=3; ngenes=1676; ntrans=1788; nexons=11860; naltexons=12795 # exon_fasta_ids ref=scaffold_3 ngroup=3; ngenes=1642; ntrans=1792; nexons=12490; naltexons=12894 # exon_fasta_ids ref=scaffold_2 ngroup=3; ngenes=1665; ntrans=1841; nexons=13482; naltexons=14649 # exon_fasta_ids ref=scaffold_1 ngroup=3; ngenes=1707; ntrans=1892; nexons=13926; naltexons=14743 # poor ids=561038; skipped ids=0 Tandy exon match types per predictor group Group Stat Same Inside Near5k Near10k Near15k Near20k Near25k Near30k Near35k Near40k Near45k Near50k Far DP_DGIL_SNO_ count 169393 7887 6962 9597 5540 4000 3130 707 868 1457 892 802 232570 DP_DGIL_SNO_ freq 1.000 0.047 0.041 0.057 0.033 0.024 0.018 0.004 0.005 0.009 0.005 0.005 1.373 DP_DGIL_SNO_ found 1.000 0.141 0.345 0.473 0.502 0.597 0.748 0.505 0.612 0.720 0.663 0.640 0.511 DP_DGIL_SNO_ anyfnd 1.000 0.249 0.391 0.516 0.528 0.655 0.758 0.560 0.624 0.731 0.692 0.680 0.568 DP_DGIL_SNO_ fnd/any 1.000 0.565 0.883 0.918 0.950 0.910 0.987 0.902 0.980 0.985 0.958 0.941 0.900 Dappu count 62433 3838 1529 2563 1462 1474 1069 143 138 534 336 242 68783 Dappu freq 1.000 0.061 0.024 0.041 0.023 0.024 0.017 0.002 0.002 0.009 0.005 0.004 1.102 Dappu found 1.000 0.109 0.181 0.348 0.373 0.596 0.674 0.175 0.087 0.581 0.750 0.628 0.601 Dappu anyfnd 1.000 0.545 0.405 0.484 0.523 0.676 0.749 0.476 0.384 0.753 0.777 0.694 0.727 Dappu fnd/any 1.000 0.200 0.445 0.718 0.712 0.881 0.900 0.368 0.226 0.771 0.966 0.905 0.828 NCBI_GNO_ count 62938 5417 3554 3843 2202 1976 1389 347 334 720 449 416 103190 NCBI_GNO_ freq 1.000 0.086 0.056 0.061 0.035 0.031 0.022 0.006 0.005 0.011 0.007 0.007 1.640 NCBI_GNO_ found 1.000 0.171 0.351 0.395 0.454 0.570 0.703 0.452 0.644 0.739 0.537 0.659 0.584 NCBI_GNO_ anyfnd 1.000 0.477 0.402 0.479 0.506 0.642 0.763 0.519 0.689 0.747 0.601 0.712 0.650 NCBI_GNO_ fnd/any 1.000 0.359 0.875 0.824 0.897 0.888 0.921 0.872 0.935 0.989 0.893 0.926 0.898 melon.% perl $td/tandynear2.perl -neardist=5000 -nearstep=10 -debug=1 $em/cele2/II/*_exons.nr.blatf8 # exon_fasta_ids ref=II ngroup=3; ngenes=10673; ntrans=10847; nexons=42398; naltexons=49843 # poor ids=199200; skipped ids=0 Tandy exon match types per predictor group Group Stat Same Inside Near5k Near10k Near15k Near20k Near25k Near30k Near35k Near40k Near45k Near50k Far Coding_transcript count 48567 1006 3423 2116 875 766 305 216 131 390 421 220 7956 Coding_transcript freq 1.000 0.021 0.070 0.044 0.018 0.016 0.006 0.004 0.003 0.008 0.009 0.005 0.164 Coding_transcript found 1.000 0.232 0.348 0.386 0.450 0.730 0.390 0.449 0.221 0.508 0.492 0.486 0.361 Coding_transcript anyfnd 1.000 0.377 0.436 0.599 0.518 0.752 0.443 0.648 0.603 0.621 0.561 0.568 0.458 Coding_transcript fnd/any 1.000 0.615 0.798 0.644 0.870 0.970 0.881 0.693 0.367 0.818 0.877 0.856 0.789 Genefinder count 45226 3751 4147 2830 1137 912 404 280 219 478 436 230 22232 Genefinder freq 1.000 0.083 0.092 0.063 0.025 0.020 0.009 0.006 0.005 0.011 0.010 0.005 0.492 Genefinder found 1.000 0.209 0.334 0.595 0.421 0.572 0.329 0.632 0.434 0.590 0.424 0.478 0.308 Genefinder anyfnd 1.000 0.390 0.378 0.650 0.484 0.647 0.438 0.646 0.530 0.621 0.541 0.526 0.348 Genefinder fnd/any 1.000 0.537 0.884 0.915 0.871 0.885 0.751 0.978 0.819 0.949 0.784 0.909 0.886 twinscan count 41029 5000 3139 2753 1161 928 340 236 192 419 322 202 12441 twinscan freq 1.000 0.122 0.077 0.067 0.028 0.023 0.008 0.006 0.005 0.010 0.008 0.005 0.303 twinscan found 1.000 0.159 0.371 0.585 0.391 0.600 0.397 0.725 0.547 0.613 0.565 0.470 0.471 twinscan anyfnd 1.000 0.285 0.431 0.644 0.465 0.629 0.471 0.729 0.562 0.613 0.615 0.564 0.569 twinscan fnd/any 1.000 0.558 0.860 0.909 0.841 0.954 0.844 0.994 0.972 1.000 0.919 0.833 0.829 #........ =cut