#!/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 [ -minalign 0.9 -nogroup ] \ scaffold_6308/dmoj_caf060210_exons.nr scaffold_6308/dmoj_caf060210_exons.nr.blatf8 =cut use strict; use warnings; use Getopt::Long; use constant SAMEBASE => 10; # for _sameloc, slop allowed in loca == locb our $NEARDIST = 8000; # was 15k our $BINSIZE = 1000 ; #was# 5000; # test with high identity: MINALIGN>0.9 our $MINEVAL = 1e-5; our $MINALIGN = 0.5; my $EQSAME=1; my $EQINSIDE=3; my $EQOVERLAP= 2; #? ignore my $EQNEAR=-2; my $EQNEAR1=-3; my $EQNEAR2=-4; my $EQFAR=-9; # add more near steps: $NEARDIST, 2x, 3x, > 3x my @eqclass= ($EQSAME, $EQNEAR, $EQNEAR1, $EQNEAR2, $EQFAR, $EQINSIDE); # same, near, far my @eqnames= ("Same","Near8k","Near15k","Near45k","Far","Inside"); # same, near, far my $groupok=1; my $debug=0; my (%hits, %ggroup, %altids, %genespan, @fasta, $skipids); #? 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, "MINALIGN=s", \$MINALIGN, "MINEVAL=s", \$MINEVAL, "NEARDIST=s", \$NEARDIST, "groupok!", \$groupok, "debug!", \$debug, ); die "usage: perl tandynear.perl [ -minalign 0.9 -nogroup ] \ scaffold_6308/dmoj_caf060210_exons.nr scaffold_6308/dmoj_caf060210_exons.nr.blatf8 " unless($optok); for( my $i=$#ARGV; $i>=0; $i--) { my $arg= $ARGV[$i]; if($arg =~ m/\.(nr|fa)$/) { push(@fasta,$arg); splice(@ARGV,$i,1); } } #warn "# fasta=@fasta\n# bltab=@ARGV\n\n"; # exit; my @nears= map{ my $kb= int((1+$_)/1000); "Near".$kb."k"; } ($NEARDIST, 2*$NEARDIST, 6*$NEARDIST); @eqnames= ("Same",@nears,"Far","Inside"); # same, near, far my %skipids=(); if($skipids) { open(F,$skipids); while(){chomp; $skipids{$_}++;} close(F); } foreach my $fasta (@fasta) { exon_fasta_ids($fasta, \%altids, \%genespan); } my $nskip=0; while(<>) { chomp; my @v=split"\t"; my($qid1,$ref,$pctid,$alen,$tb,$te,$eval,$bits)=@v[0,1,2,3,8,9,10,11]; my $tstrand=1; if($tb>$te) { $tstrand=-1; ($tb,$te)= ($te,$tb); } my %didqid; # any dups in altids? my $altids= $altids{$qid1} || []; foreach my $qid0 ($qid1, @$altids) { my($qid,$qb,$qe)=split(/[:-]/,$qid0); # exon id,start,stop my($gid) = $qid =~ m/^(\w+)/; # using ID.exnum syntax my($ggroup)= ($groupok) ? $qid =~ m/^(\D+)/ : ("all"); # BAD for non-grouping-prefixed IDs $nskip++ and next if($skipids{$gid}); # gene/transcript id here; or both? (exonid) $nskip++ and next if($skipids{$qid}); # both? (exonid) my ($gb,$ge,$gref)= $genespan{$gid} ? @{$genespan{$gid}} : (0,0,0); my $elen = 1+abs($qe-$qb); my $palign= $alen / $elen; next if($palign < $MINALIGN); next if($eval > $MINEVAL); # save exon match locs; after read all, recount if near/far has other prediction at sameloc my $eq= -1; if( _sameloc($ref,$qb,$qe,0, $ref,$tb,$te,0) ) { $eq= $EQSAME; # } elsif (_isoverlap($qb,$qe,$tb,$te)) { # need _isoverlap also ? # # problem was missing $tstrand fix; now few of these; not same exons # # warn "overlap: $ref:$qid $qb-$qe <> $tb-$te \n"; # $eq= $EQOVERLAP; } elsif ($ge>0 && $gref eq $ref && _isoverlap($gb,$ge,$tb,$te)) { # gene span $eq= $EQINSIDE; # } elsif (_isnear($qb,$qe,$tb,$te)) { # always same $ref here # $eq= $EQNEAR; # } else { # far match # $eq= $EQFAR; } else { my $dist= _isnear2($qb,$qe,$tb,$te); if($dist < $NEARDIST) { $eq= $EQNEAR; } # 8k x 2 elsif($dist < 2*$NEARDIST) { $eq= $EQNEAR1; } #16k x 2 elsif($dist < 6*$NEARDIST) { $eq= $EQNEAR2; } #48k x 2 else { $eq= $EQFAR; } } my $loc= [$tb,$te,$qid0,$ref,$eq]; ## 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}++; $ggroup{$ggroup}{total}++; } } warn "# skipped ids n=$nskip\n" if($debug && $nskip>0); my @ggroup= sort keys %ggroup; my $manygroup= @ggroup>1; my(%predfound, %predolap, %predcount); foreach my $ggroup (@ggroup) { $predcount{$ggroup}{$EQSAME}= $predolap{$ggroup}{$EQSAME}= $predfound{$ggroup}{$EQSAME}= 1; my @eqlist= @eqclass; shift @eqlist; # drop EQSAME foreach my $eq (@eqlist) { $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} } if($hits{$ggroup}{$eq} and $hits{$ggroup}{$eq}{$ib}); } foreach my $dp (@duplocs) { # match exons from dupl class and same? next if($diddup{$dp}++); my($dqb,$dqe,$dqid0,$dref,$deq)= @$dp; next unless($deq == $eq); # redundant 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; foreach my $sp (@samelocs) { my($sqb,$sqe,$sqid0,$sref,$seq)= @$sp; next unless($sref eq $dref && $seq == $EQSAME); if( _sameloc($sref,$sqb,$sqe,0, $dref,$dqb,$dqe,0) ) { $found= 1; last; } elsif( _isoverlap($sqb,$sqe,$dqb,$dqe) ) { $olap= 1; } } my $otherfound=0; my @anysame=(); if($manygroup) { 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; } } } $predfound{$ggroup}{$eq}++ if $found; $predolap{$ggroup}{$eq}++ if ($found || $otherfound); $predcount{$ggroup}{$eq}++; } } } print "Tandy exon match types per predictor group\n"; print join("\t","Group ", "Total",@eqnames),"\n"; foreach my $ggroup (@ggroup) { my $total = $ggroup{$ggroup}{total}; my $tsame = $ggroup{$ggroup}{$EQSAME} || 1; my @counts= map{ $ggroup{$ggroup}{$_} || 0 } @eqclass; my @freqs= map{ sprintf("%.3f", $_ / $tsame ); } @counts; my $sgroup= sprintf("%-15s",$ggroup); my @foundfreq= map{ my $pc= $predcount{$ggroup}{$_} || 1; sprintf("%.3f", $predfound{$ggroup}{$_} / $pc); } @eqclass; my @anyfnd= map{ my $pc= $predcount{$ggroup}{$_} || 1; sprintf("%.3f", $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,$total,@counts),"\n"; 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 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 _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 _sameloc { my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_; # $aor ||= 0; $bor ||= 0; $ab ||=0; $ae ||=0; $bb ||=0; $be ||=0; #? need return # $aor eq $bor && # strand/orient (abs($ab-$bb) < SAMEBASE) && (abs($ae-$be) < SAMEBASE) ; } sub exon_id_crossref { my($altids)= @_; ## 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 $alist= $altids->{$xid} or next; foreach my $axid (@$alist) { $xaltid{$xid}{$axid}= 1; $xaltid{$axid}{$xid}= 1; $xaltid{$axid}{$axid}= 1; # ensure self match } } # 1b. update alt lists with all xrefs; must keep xloc as part of xid ## ? we are missing self/xid in altlist ?? why? 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)= @_; $altids ||= {}; $genespan ||= {}; my(%gloc,%didexongloc); my($nexons, $naltexons, $ntrans, $ngenes, $atref, $atstrand)=(0)x10; my $newaltids= {}; open(GREP, "grep '^>' $queryfa|") or return $altids; # die "grep $queryfa"; while(){ my($xid)=m/>(\S+)/ or next; my ($ref,$strand,$loc)=(0,0,"",""); if(m/loc=([^;\s]+)/) { ($ref,$loc,$strand)= split(/:/,$1); ($atref, $atstrand)= ($ref,$strand); } 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++; #>GI_DGIL_SNO_28460615.1:13347039-13347238 type=CDS:DGIL_SNO; loc=scaffold_6473:13347039-13347238:-1; altids=GI_DGIL_SNO_28460615.1:13347039-13347238,GI_NCBI_GNO_32327927.1:13347039-13347238,GI_PACH_GMP_15216073.1:13347039-13347238,GI_PACH_GMP_15216074.1:13347039-13347238,GI_PACH_GMP_15216075.1:13347039-13347238,GI_RGUI_GID_mRNA_40009118.9:13347039-13347238,GLEAN_19757.1:13347039-13347238,dmoj_GLEANR_15521.1:13347039-13347238 if(m/altids=(\S+)/) { my $aid= $1; $aid =~ s/pm,.*$//; # no partials here my @altids= grep{$xid ne $_} split(/,/, $aid); #? exclude $xid here or not $newaltids->{$xid}= \@altids; if(0) { foreach my $ad (@altids) { next if($didexongloc{$ad}); $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); #{ ($agid) = $xid =~ m/^(\w+)/; } # where this error? $gloc{$agid} or $gloc{$agid}=[]; push @{$gloc{$agid}}, [$axn,$ref,$axb,$axe,$strand]; $didexongloc{$ad}++; } } } } ## 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($newaltids); if(1) { my $missinstrand=0; foreach my $xid (sort keys %{$newaltids}) { # NOT for all ref/scaffold, just input one next if($didexongloc{$xid}); my $alist= $newaltids->{$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? $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"); 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; } # ... 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]; $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; warn "# genespan errs=$nerr\n" if($debug && $nerr>0); return ($altids,$genespan); } __END__ =item results # dpulex, all of clean genome # exon_fasta_ids ref=scaffold_999 ngroup=3; ngenes=7; ntrans=7; nexons=17; naltexons=1 # exon_fasta_ids ref=scaffold_998 ngroup=1; ngenes=1; ntrans=1; nexons=2; naltexons=0 ... # exon_fasta_ids ref=scaffold_10 ngroup=3; ngenes=1366; ntrans=1430; nexons=6564; naltexons=4191 # exon_fasta_ids ref=scaffold_1 ngroup=3; ngenes=1707; ntrans=1892; nexons=10220; naltexons=6096 Tandy exon match types per predictor group Group Total Same Near Far Inside DP_DGIL_SNO_ 219436 161009 14844 40068 3515 DP_DGIL_SNO_ freq 1.000 0.092 0.249 0.022 DP_DGIL_SNO_ found 1.000 0.345 0.272 0.253 DP_DGIL_SNO_ anyfnd 1.000 0.383 0.303 0.269 DP_DGIL_SNO_ fnd/any 1.000 0.901 0.899 0.939 Dappu 166033 134424 9550 19982 2077 Dappu freq 1.000 0.071 0.149 0.015 Dappu found 1.000 0.265 0.316 0.156 Dappu anyfnd 1.000 0.428 0.453 0.358 Dappu fnd/any 1.000 0.619 0.697 0.436 NCBI_GNO_ 236201 177386 16726 39025 3064 NCBI_GNO_ freq 1.000 0.094 0.220 0.017 NCBI_GNO_ found 1.000 0.445 0.406 0.309 NCBI_GNO_ anyfnd 1.000 0.486 0.462 0.381 NCBI_GNO_ fnd/any 1.000 0.915 0.879 0.810 # Tandy exon match types per predictor group # Group Total Same Near Far # DP_DGIL_SNO_ 219436 161009 17469 40154 # DP_DGIL_SNO_ freq 1.000 0.108 0.249 # DP_DGIL_SNO_ found 1.000 0.343 0.272 # DP_DGIL_SNO_ anyfnd 1.000 0.378 0.303 # # Dappu 166033 134424 11356 20209 # Dappu freq 1.000 0.084 0.150 # Dappu found 1.000 0.250 0.313 # Dappu anyfnd 1.000 0.419 0.452 # # NCBI_GNO_ 236201 177386 19571 39196 # NCBI_GNO_ freq 1.000 0.110 0.221 # NCBI_GNO_ found 1.000 0.426 0.405 # NCBI_GNO_ anyfnd 1.000 0.472 0.461 # dpulex, all of clean genome, >= 90% alignment perl $td/tandynear.perl -minalign 0.9 $em/daphe/scaffold_*/dpulex1_exons.nr.blatf8 $em/daphe/scaffold_*/dpulex1_exons.nr Tandy exon match types per predictor group Group Total Same Near Far DP_DGIL_SNO_ 196384 161008 10669 24517 DP_DGIL_SNO_ freq 1.000 0.066 0.152 DP_DGIL_SNO_ found 1.000 0.540 0.419 DP_DGIL_SNO_ anyfnd 1.000 0.589 0.461 Dappu 153514 134424 6249 12828 Dappu freq 1.000 0.046 0.095 Dappu found 1.000 0.439 0.473 Dappu anyfnd 1.000 0.723 0.666 NCBI_GNO_ 214521 177386 11736 25381 NCBI_GNO_ freq 1.000 0.066 0.143 NCBI_GNO_ found 1.000 0.686 0.599 NCBI_GNO_ anyfnd 1.000 0.752 0.664 # # subset of full, clean dpulex annots # perl $td/tandynear.perl $em/daphe/scaffold_?/dpulex1_exons.nr.blatf8 $em/daphe/scaffold_?/dpulex1_exons.nr # # Tandy exon match types per predictor group # Group Total Same Near Far # DP_DGIL_SNO_ 50920 33691 3167 14054 # DP_DGIL_SNO_ freq 1.000 0.094 0.417 # DP_DGIL_SNO_ found 1.000 0.404 0.310 # DP_DGIL_SNO_ anyfnd 1.000 0.439 0.338 # # Dappu 39711 30787 2114 6808 # Dappu freq 1.000 0.069 0.221 # Dappu found 1.000 0.257 0.442 # Dappu anyfnd 1.000 0.436 0.567 # # NCBI_GNO_ 56499 38035 3823 14639 # NCBI_GNO_ freq 1.000 0.101 0.385 # NCBI_GNO_ found 1.000 0.435 0.443 # NCBI_GNO_ anyfnd 1.000 0.479 0.521 # # # using MINALIGN=0.9 reduces near,far; increases found rates (to .7,.8) but # # ratio of group-found / any-found is ~same # perl $td/tandynear.perl -MINALIGN 0.9 $em/daphe/scaffold_?/dpulex1_exons.nr.blatf8 $em/daphe/scaffold_?/dpulex1_exons.nr # # Tandy exon match types per predictor group # Group Total Same Near Far # DP_DGIL_SNO_ 44267 33691 1891 8685 # DP_DGIL_SNO_ found 0.000 0.654 0.470 # DP_DGIL_SNO_ anyfnd 0.000 0.707 0.506 # # Dappu 36575 30787 1118 4670 # Dappu found 0.000 0.473 0.613 # Dappu anyfnd 0.000 0.784 0.761 # # NCBI_GNO_ 50229 38035 2189 10005 # NCBI_GNO_ found 0.000 0.734 0.614 # NCBI_GNO_ anyfnd 0.000 0.798 0.692 #.................................................... # # Dros. moj., all euchromatin, predictor subset # perl $td/tandynear.perl $em/dmoj2/scaffold_6*/dmoj_caf060210_exons.nr $em/dmoj2/scaffold_6*/dmoj_caf060210_exons.nr.blatf8 # Tandy exon match types per predictor group # Group Total Same Near Far # GI_BREN_NSC_ 79484 55327 2569 21534 # GI_BREN_NSC_ freq 1.000 0.046 0.389 # GI_BREN_NSC_ found 1.000 0.253 0.034 * 5% off tandem modelling # GI_BREN_NSC_ anyfnd 1.000 0.306 0.052 # # GI_DGIL_SNO_ 187083 115342 4594 67028 # GI_DGIL_SNO_ freq 1.000 0.040 0.581 # GI_DGIL_SNO_ found 1.000 0.283 0.054 * 5% off tandem modelling # GI_DGIL_SNO_ anyfnd 1.000 0.325 0.061 # # GI_NCBI_GNO_ 125723 113050 4260 8374 * oddly low Far count # GI_NCBI_GNO_ freq 1.000 0.038 0.074 # GI_NCBI_GNO_ found 1.000 0.404 0.108 * 3% off tandem modelling # GI_NCBI_GNO_ anyfnd 1.000 0.435 0.120 # # GLEAN_ 139399 113387 4639 21317 # GLEAN_ freq 1.000 0.041 0.188 # GLEAN_ found 1.000 0.386 0.059 * 6% off tandem modelling # GLEAN_ anyfnd 1.000 0.441 0.096 # Dros. moj., all euchrom, subset predictor, with inside gene class melon.% perl $td/tandynear.perl $em/dmoj2/scaffold_6*/dmoj_caf060210_exons.nr $em/dmoj2/scaffold_6*/dmoj_caf0602 10_exons.nr.blatf8 # exon_fasta_ids ref=scaffold_6680 ngenes=12987; ntrans=13738; nexons=36091; naltexons=31448 # exon_fasta_ids ref=scaffold_6654 ngenes=1571; ntrans=1688; nexons=4622; naltexons=4681 # exon_fasta_ids ref=scaffold_6541 ngenes=1188; ntrans=1213; nexons=2154; naltexons=696 # exon_fasta_ids ref=scaffold_6540 ngenes=18096; ntrans=19092; nexons=53273; naltexons=49477 # exon_fasta_ids ref=scaffold_6500 ngenes=14912; ntrans=15450; nexons=39866; naltexons=35724 # exon_fasta_ids ref=scaffold_6498 ngenes=1342; ntrans=1397; nexons=3647; naltexons=2163 # exon_fasta_ids ref=scaffold_6496 ngenes=14180; ntrans=14903; nexons=43049; naltexons=42059 # exon_fasta_ids ref=scaffold_6482 ngenes=1325; ntrans=1358; nexons=2941; naltexons=1824 # exon_fasta_ids ref=scaffold_6473 ngenes=7766; ntrans=8233; nexons=20059; naltexons=14977 # exon_fasta_ids ref=scaffold_6359 ngenes=2068; ntrans=2214; nexons=5830; naltexons=4577 # exon_fasta_ids ref=scaffold_6328 ngenes=2153; ntrans=2318; nexons=6239; naltexons=5576 # exon_fasta_ids ref=scaffold_6308 ngenes=1660; ntrans=1781; nexons=4385; naltexons=3728 Tandy exon match types per predictor group Group Total Same Near Far Inside GI_BREN_NSC_ 79484 55327 2045 21530 582 GI_BREN_NSC_ freq 1.000 0.037 0.389 0.011 GI_BREN_NSC_ found 1.000 0.270 0.034 0.167 GI_BREN_NSC_ anyfnd 1.000 0.325 0.052 0.210 GI_BREN_NSC_ fnd/any 1.000 0.831 0.652 0.795 GI_DGIL_SNO_ 187083 115342 3981 67027 733 GI_DGIL_SNO_ freq 1.000 0.035 0.581 0.006 GI_DGIL_SNO_ found 1.000 0.283 0.054 0.237 GI_DGIL_SNO_ anyfnd 1.000 0.327 0.061 0.258 GI_DGIL_SNO_ fnd/any 1.000 0.863 0.891 0.921 GI_NCBI_GNO_ 125723 113050 3854 8356 463 * odd low Far count GI_NCBI_GNO_ freq 1.000 0.034 0.074 0.004 GI_NCBI_GNO_ found 1.000 0.419 0.108 0.233 GI_NCBI_GNO_ anyfnd 1.000 0.447 0.120 0.283 GI_NCBI_GNO_ fnd/any 1.000 0.938 0.898 0.824 GLEAN_ 139399 113387 4102 21301 609 GLEAN_ freq 1.000 0.036 0.188 0.005 GLEAN_ found 1.000 0.402 0.059 0.241 GLEAN_ anyfnd 1.000 0.448 0.096 0.350 GLEAN_ fnd/any 1.000 0.897 0.610 0.690 # Dros. moj., subset euchromatin, all predictor perl $td/tandynear.perl $em/dmoj1/scaffold_6*/dmoj_caf060210_exons.nr $em/dmoj1/scaffold_6*/dmoj_caf060210_exons.nr.blatf8 | & more # exon_fasta_ids ref=scaffold_6680 ngroup=11; ngenes=27884; ntrans=38764; nexons=115848; naltexons=480864 # exon_fasta_ids ref=scaffold_6541 ngroup=11; ngenes=1875; ntrans=1996; nexons=3592; naltexons=3653 # exon_fasta_ids ref=scaffold_6473 ngroup=11; ngenes=15372; ntrans=20482; nexons=58949; naltexons=234972 Tandy exon match types per predictor group Group Total Same Near Far Inside GI_BATZ_CNA_ 16916 11402 275 5135 104 GI_BATZ_CNA_ freq 1.000 0.024 0.450 0.009 GI_BATZ_CNA_ found 1.000 0.225 0.025 0.183 GI_BATZ_CNA_ anyfnd 1.000 0.320 0.062 0.221 GI_BATZ_CNA_ fnd/any 1.000 0.705 0.402 0.826 GI_BREN_NSC_ 29738 21875 541 7196 126 GI_BREN_NSC_ freq 1.000 0.025 0.329 0.006 GI_BREN_NSC_ found 1.000 0.288 0.024 0.190 GI_BREN_NSC_ anyfnd 1.000 0.375 0.054 0.286 GI_BREN_NSC_ fnd/any 1.000 0.768 0.450 0.667 GI_DGIL_SNO_ 59305 36330 881 21882 212 GI_DGIL_SNO_ freq 1.000 0.024 0.602 0.006 GI_DGIL_SNO_ found 1.000 0.279 0.036 0.226 GI_DGIL_SNO_ anyfnd 1.000 0.338 0.047 0.255 GI_DGIL_SNO_ fnd/any 1.000 0.826 0.761 0.889 GI_EISE_CEX_ 88374 82910 861 4021 582 GI_EISE_CEX_ freq 1.000 0.010 0.048 0.007 GI_EISE_CEX_ found 1.000 0.261 0.019 0.301 GI_EISE_CEX_ anyfnd 1.000 0.454 0.039 0.378 GI_EISE_CEX_ fnd/any 1.000 0.575 0.490 0.795 GI_EISE_CGW_ 154467 150220 1344 1887 1016 GI_EISE_CGW_ freq 1.000 0.009 0.013 0.007 GI_EISE_CGW_ found 1.000 0.365 0.144 0.322 GI_EISE_CGW_ anyfnd 1.000 0.476 0.144 0.531 GI_EISE_CGW_ fnd/any 1.000 0.766 0.996 0.606 GI_NCBI_GNO_ 68740 63979 1259 3341 161 GI_NCBI_GNO_ freq 1.000 0.020 0.052 0.003 GI_NCBI_GNO_ found 1.000 0.442 0.032 0.435 GI_NCBI_GNO_ anyfnd 1.000 0.496 0.060 0.491 GI_NCBI_GNO_ fnd/any 1.000 0.891 0.543 0.886 GI_PACH_GMP_ 137234 135429 610 807 388 GI_PACH_GMP_ freq 1.000 0.005 0.006 0.003 GI_PACH_GMP_ found 1.000 0.018 0.073 0.000 GI_PACH_GMP_ anyfnd 1.000 0.525 0.119 0.479 GI_PACH_GMP_ fnd/any 1.000 0.034 0.615 0.000 GI_RGUI_GID_ 76460 66332 814 8778 536 GI_RGUI_GID_ freq 1.000 0.012 0.132 0.008 GI_RGUI_GID_ found 1.000 0.294 0.028 0.397 GI_RGUI_GID_ anyfnd 1.000 0.386 0.059 0.476 GI_RGUI_GID_ fnd/any 1.000 0.761 0.482 0.835 GLEAN_ 104873 91665 1831 11207 170 GLEAN_ freq 1.000 0.020 0.122 0.002 GLEAN_ found 1.000 0.430 0.031 0.041 GLEAN_ anyfnd 1.000 0.497 0.076 0.247 GLEAN_ fnd/any 1.000 0.865 0.406 0.167 TRdmoj_ 135272 130424 1541 2421 886 TRdmoj_ freq 1.000 0.012 0.019 0.007 TRdmoj_ found 1.000 0.500 0.114 0.281 TRdmoj_ anyfnd 1.000 0.539 0.139 0.468 TRdmoj_ fnd/any 1.000 0.929 0.822 0.600 dmoj_GLEANR_ 111404 98647 1999 10571 187 dmoj_GLEANR_ freq 1.000 0.020 0.107 0.002 dmoj_GLEANR_ found 1.000 0.454 0.030 0.000 dmoj_GLEANR_ anyfnd 1.000 0.509 0.076 0.299 dmoj_GLEANR_ fnd/any 1.000 0.893 0.389 0.000 #.. old.... # perl $td/tandynear.perl $em/dmoj1/scaffold_6*/dmoj_caf060210_exons.nr $em/dmoj1/scaffold_6*/dmoj_caf060210_exons.nr.blatf8 # # fasta=/bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6680/dmoj_caf060210_exons.nr /bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6541/dmoj_caf060210_exons.nr /bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6473/dmoj_caf060210_exons.nr # # bltab=/bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6473/dmoj_caf060210_exons.nr.blatf8 /bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6541/dmoj_caf060210_exons.nr.blatf8 /bio/bio-grid/mb/EVidenceModeler/dmoj1/scaffold_6680/dmoj_caf060210_exons.nr.blatf8 # # Tandy exon match types per predictor group # Group Total Same Near Far # GI_BATZ_CNA_ 16916 11402 379 5135 # GI_BATZ_CNA_ found 0.000 0.214 0.025 # GI_BATZ_CNA_ anyfnd 0.000 0.293 0.062 # # GI_BREN_NSC_ 29738 21875 664 7199 # GI_BREN_NSC_ found 0.000 0.271 0.024 # GI_BREN_NSC_ anyfnd 0.000 0.360 0.054 # # GI_DGIL_SNO_ 59305 36330 1093 21882 # GI_DGIL_SNO_ found 0.000 0.269 0.036 # GI_DGIL_SNO_ anyfnd 0.000 0.322 0.047 # # GI_EISE_CEX_ 88374 82910 1443 4021 * includes alt transcripts (?) # GI_EISE_CEX_ found 0.000 0.277 0.019 # GI_EISE_CEX_ anyfnd 0.000 0.423 0.039 # # GI_EISE_CGW_ 154467 150220 2356 1891 * includes alt transcripts (?) # GI_EISE_CGW_ found 0.000 0.347 0.143 # GI_EISE_CGW_ anyfnd 0.000 0.501 0.144 # # GI_NCBI_GNO_ 68740 63979 1409 3352 # GI_NCBI_GNO_ found 0.000 0.445 0.032 * 11% off tandem modelling (.45/.50) # GI_NCBI_GNO_ anyfnd 0.000 0.500 0.059 # # GI_PACH_GMP_ 137234 135429 948 857 * includes alt transcripts (?) # GI_PACH_GMP_ found 0.000 0.012 0.069 # GI_PACH_GMP_ anyfnd 0.000 0.512 0.137 # # GI_RGUI_GID_ 76460 66332 1232 8896 # GI_RGUI_GID_ found 0.000 0.334 0.032 # GI_RGUI_GID_ anyfnd 0.000 0.422 0.063 # # GLEAN_ 104873 91665 1998 11210 # GLEAN_ found 0.000 0.397 0.031 # GLEAN_ anyfnd 0.000 0.476 0.076 # # TRdmoj_ 135272 130424 2414 2434 * includes alt transcripts (?) # TRdmoj_ found 0.000 0.423 0.114 # TRdmoj_ anyfnd 0.000 0.516 0.138 # # dmoj_GLEANR_ 111404 98647 2179 10578 # dmoj_GLEANR_ found 0.000 0.417 0.029 # dmoj_GLEANR_ anyfnd 0.000 0.492 0.076 #.................................................... # Dros. mel. -- redo tandy w/ all CAF1 predictors # ** includes alt splice exons?; :( ids are all FBtr transcript ; cant tell gene span from this # should reframe _isnear to use gene-span, distinguish inside/outside of genespan perl $td/tandynear.perl -skip dmel.ncrna.idlist -nogroup $em/dmel4/{2,3,4,X,U}*/dmel4_exons.nr.blatf8 $em/dmel4/{2,3,4,X,U}*/dmel4_exons.nr # exon_fasta_ids ref=X ngroup=1; ngenes=2058; ntrans=3138; nexons=10009; naltexons=140 # exon_fasta_ids ref=4 ngroup=1; ngenes=89; ntrans=205; nexons=875; naltexons=15 # exon_fasta_ids ref=3R ngroup=1; ngenes=3060; ntrans=4753; nexons=15996; naltexons=174 # exon_fasta_ids ref=3L ngroup=1; ngenes=2376; ntrans=3737; nexons=12089; naltexons=128 # exon_fasta_ids ref=2R ngroup=1; ngenes=2524; ntrans=4057; nexons=13336; naltexons=2152 # exon_fasta_ids ref=2L ngroup=1; ngenes=2333; ntrans=3565; nexons=11483; naltexons=647 # genespan errs=1 Tandy exon match types per predictor group Group Total Same Near Far Inside all 75589 63175 3537 8609 268 all freq 1.000 0.056 0.136 0.004 all found 1.000 0.701 0.891 0.086 # Dros.mel, -minalign 0.9 Tandy exon match types per predictor group Group Total Same Near Far Inside all 74402 63175 2918 8201 108 all freq 1.000 0.046 0.130 0.002 all found 1.000 0.847 0.935 0.204 # Tandy exon match types per predictor group # Group Total Same Near Far Inside # all 113796 65000 24650 23553 593 # all freq 1.000 0.379 0.362 0.009 * odd high near freq; ?? artifact of something? # all found 1.000 0.951 0.948 0.567 # * non-coding exons/genes are included in dmel4 annots, some high-dupl trna/ncrna # * 2R has 145% Near > Same ** ; 2L has 20% Near; X has 7% Near; 3L,3R have 4% Near/Same; ?????? # ^^^ this is the problem; 2R has some 2200 matches to 280 ncRNA genes; # 3L has 465 matches to 100 ncRNA genes; 2L has 911 matches to 138 ncRNA #.................................................... # C. elegans, with alt transcripts parsing perl $td/tandynear.perl -nogroup $em/cele1/{I,V,X}*/cele1_exons.nr.blatf8 $em/cele1/{I,V,X}*/cele1_exons.nr # exon_fasta_ids ngenes=2743; ntrans=3788; nexons=24169; naltexons=14974 # exon_fasta_ids ngenes=4930; ntrans=6026; nexons=30065; naltexons=10418 # exon_fasta_ids ngenes=3210; ntrans=4505; nexons=23568; naltexons=17456 # exon_fasta_ids ngenes=2578; ntrans=3951; nexons=20580; naltexons=14278 # exon_fasta_ids ngenes=3401; ntrans=4675; nexons=23497; naltexons=13082 # exon_fasta_ids ngenes=2785; ntrans=4092; nexons=22443; naltexons=14234 Tandy exon match types per predictor group Group Total Same Near Far Inside all 253063 226386 10849 14161 1667 all freq 1.000 0.048 0.063 0.007 all found 1.000 0.389 0.315 0.142 # C. elegans, -minalign 0.9 Tandy exon match types per predictor group Group Total Same Near Far Inside all 240318 226386 5960 7469 503 all freq 1.000 0.026 0.033 0.002 all found 1.000 0.690 0.578 0.384 #.................................................... # Apis mell. perl $td/tandynear.perl -nogroup $em/amel4/chr*/amel4_exons.nr.blatf8 $em/amel4/chr*/amel4_exons.nr Tandy exon match types per predictor group Group Total Same Near Far all 47565 46225 719 600 all freq 1.000 0.016 0.013 all found 1.000 0.234 0.323 =item fixed bug odd: not finding unpredicted exons I know are there from map views ... ** problem was we used exon-pred location (qb,qe) not genome match loc (tb,te) correct table with genome (tb,te) locations melon.% cat $em/daphd/scaffold_4/dpulex1_exons.nr.blatf8 | perl $td/tandynear.perl -fa $em/daphd/scaffold_4/dpulex1_exons.nr Tandy exon match types per predictor group Group Total Same Near Far DP_DGIL_SNO_ 35241 16328 3930 14983 DP_DGIL_SNO_ found 0.000 0.419 0.151 DP_DGIL_SNO_ anyfnd 0.000 0.465 0.158 Dappu 9402 7144 1292 966 Dappu found 0.000 0.220 0.226 << the missing Dappu predictions for tandems Dappu anyfnd 0.000 0.490 0.310 NCBI_GNO_ 12149 7840 2022 2287 NCBI_GNO_ found 0.000 0.453 0.221 NCBI_GNO_ anyfnd 0.000 0.485 0.254 melon.% cat scaffold_6680/dmoj_caf060210_exons.nr.blatf8 | perl $td/tandynear.perl -fa scaffold_6680/dmoj_caf060210_exons.nr Tandy exon match types per predictor group Group Total Same Near Far GI_BREN_NSC_ 18007 14817 450 2740 GI_BREN_NSC_ found 0.000 0.320 0.016 GI_BREN_NSC_ anyfnd 0.000 0.391 0.038 GI_DGIL_SNO_ 35983 24270 782 10931 GI_DGIL_SNO_ found 0.000 0.263 0.039 GI_DGIL_SNO_ anyfnd 0.000 0.317 0.047 GI_NCBI_GNO_ 23024 21493 703 828 GI_NCBI_GNO_ found 0.000 0.427 0.016 GI_NCBI_GNO_ anyfnd 0.000 0.444 0.024 GLEAN_ 25734 21967 697 3070 GLEAN_ found 0.000 0.374 0.038 GLEAN_ anyfnd 0.000 0.438 0.072 =cut