#!/usr/bin/perl # tdprotnear2.perl use strict; use warnings; use Getopt::Long; =item inputs 1. source gene id list $dp-protnear15k.id1 or $dp-protnear15k.ids (paired) : this one mixes predictors 2. source gene gff (add CDS gff?): $dp-$predictor-genes.note.gff.gz 3. target predictor gene gff (add CDS gff??) 4. tandy marked tandem exons: need for tandems missed by predictors -- maybe all needed here is subset of tandy HSPs with no source gene from input predictors -- need to extract tandy info of exon set: .. best source gene set? =item outputs -- counts -- any marked gene gff? =item version 2 * need to use input of pred-exon search (tandy) to find all case 4 errs (TDMISSED) * change to match gene-exons not gene-span ** * change to count below gene model error classes, instead of inside/outside gene span * try to process, count all predictor classes, rather than one reference class (any/all can be reference and target/test classes) =item notes Protein near gene sets: for each predictor w/ duplicates, near-group, find overlap of genes, CDS-exons (?) among all/subset other predictors. Count these: 1. substantial agreement (same and dupl. genes, mostly same exons) 2. same-gene overlap/no-overlap (missed source gene?) 3. extended-overlap (other has one gene/exon overlap for source and dupl genes) 3a. all dupl. exons matched; 3b. dupl. exons skipped 4. tandem-missed (other found same source gene, but no model for tandem) 5. ? foreach source predictor gene id - get source gene GFF (and CDS GFF?) == SG1 - get all source tandem genes (foreach nearstep or for nearregion?) == @SDn - overlapfilter target predictor (gene and/or CDS) GFF == OG1, @ODn - count: SG1 x OG1 [1,2], SDn x OG1 [3 a,b], SDn x ODn [4], ? SG1 x ODn [5?] Start with limited protnear geneid set as above: $dp-protnear15k.id1 gene model error classes: - 2. SCMISSED: no gm at source gene (can also be tandem, ~ 4 TDMISSED) - 3a. DOUBLED : gm contains 2+ gene-exon sets, double size of tandem dupl; protein match has part matches to other full genes at low/high parts -- common among ab-initio predictors - 3b. SKIPOVER: gm has 1 gene-exon set, but inside of gene span are 1+ other ~complete gene-exon sets -- common among gene-mapping modellers - 4. TDMISSED: no gm where other predictors find tandem gene (not quite same as skipover) differs from 2 (no models) in that 1 gm in tandem region is found, not others -- includes cases where all gene predictors miss duplicate exon set, need pred-exon search to find. =item case 1 : one source gene, target gene set set pdid=/bio/bio-grid/dpulex/prots/tanprot/dros/idlists set aadros=/bio/bio-grid/dpulex/prots/tanprot/dros set dp=dgri set predict=gnomon set pother=genewise foreach dp ($dspp) echo ========[ $dp ]============================ set othergff=$aadros/$dp*$pother*.note.genes.gz set skipids=$em/${dp}1/$dp*genewise.altgn2 $td/tdprotnear.perl -debug=1 -pctover 75 \ -tandemid $aadros/idlists/$dp*$predict*.protnear15k.ids \ -source $aadros/$dp*$predict*.note.genes.gz \ -other $aadros/$dp*$pother*.note.genes.gz \ -skipid $skipids end =item case 2: all predictors versus one source gene set set aadros=/bio/bio-grid/dpulex/prots/tanprot/dros cd $aadros set dp=dere $td/tdprotnear.perl -debug=1 -pctover 75 \ -tandemid blastnear/${dp}_aapred6.blastnear15 \ -refgroup dere_GLEANR_ -source ${dp}_dmelhsp_note_all.genes.gz \ -other ${dp}_dmelhsp_note_all.genes.gz =item version 2 test: ** must have -refgroup and == gene ID - digits set dp=dere $td/tdprotnear2.perl -debug=1 -pctover 75 \ -tandemid blastnear/${dp}_aapred6.blastnear21e \ -refgroup GG_NCBI_GNO_ \ -gff $em/${dp}1/${dp}*predict.gff.gz \ ## ${dp}_dmelhsp_note_all.genes.gz << no good; need allpredmodels.gff w/ CDS exons =item problems * alt-tr versus nearby duplicates: is that still a problem? - gnomon has some alt-tr, more for near-dmel - genewise has lots (40k mrna/13k genes for dvir?, 20k for others) * gnomon has utrs that expand gene size: use CDS overlap instead? - problem for Inside/Outside/Overlap classifying based on gene span =cut # use constant SAMEBASE => 10; # for _sameloc, slop allowed in loca == locb # try wider match? use constant SAMEBASE => 27; # for _sameloc, slop allowed in loca == locb use constant BINSIZE => 500; # 1000 ; #was# 5000; our $MINEVAL = 1e-5; our $MINBITSCORE = 0; our $MINALIGN = 0.1; # 0.5; our $NEARDIST = 15000; # was 15k our $NEARSTEPS= 3; our $EQSAME=1; our $EQOVERLAP=2; our $EQINSIDE=3; our $EQOUTSIDE=4; our $EQPARTOVER=5; our $EQNEAR=-2; our $EQFAR=-999; our $EQNONE=-1; # these are new gene model error classes: 3a,3b,2,4 below our $EQDOUBLED=6; our $EQSKIPOVER=7; our $EQSCMISSED=8; our $EQTDMISSED=9; our $EQALTTR=10; my $groupok=1; my $pctover=20; my $pctinside=50; my $mrnatype= $ENV{mrnatype} || "mRNA"; my $exontype= $ENV{exontype} || "CDS"; my $useexons=1; my ($clustersize, $clustertop, $debug, $dumpgenes, $nskip, $npoor, ) = (0) x 20; my ($ok, $tandemidpairs, $blastptable, $sourcegff, $othergff, @clusize, @neargroup, %ggroup, %genespan, %exons, $skipids, $refgroup, $keeponly); my $skipchr='^U'; # for dpse only ?? better to fix input data my $skipgrp=""; my $optok= GetOptions( "tandemidpairs=s", \$tandemidpairs, "gff|sourcegff=s", \$sourcegff, #? dont need to distinguish source/other gff? # "othergff=s", \$othergff, # same as source gff now; stdin ok "refgroup=s", \$refgroup, # use this against gene IDs to separate source/other genes # "blastptable=s", \$blastptable, # drop this input; use blastnear for dup gene pairs # "eval=s",\$MINEVAL, "bitscore=s",\$MINBITSCORE, "align=s",\$MINALIGN, "skipids=s", \$skipids, # == +keepids; for TE genes, etc. "mrnatype=s", \$mrnatype,"exontype=s", \$exontype, # "keeponly=s", \$keeponly, # regex for id-group # "NEARDIST=i", \$NEARDIST, "NEARSTEPS=i", \$NEARSTEPS, "neargroup=i", \@neargroup, "pctoverlap=i", \$pctover, "pctinside=i", \$pctinside, "clustersize=i", \@clusize, "dumpgenes=s", \$dumpgenes, "skipchr=s", \$skipchr, "skipgrp=s", \$skipgrp, "debug=i", \$debug, ); die "usage: perl tdprotnear.perl [ options ] -tandemidpairs dsec-protnear15k.ids : as ID1.ID2 list of tandem genes -gff allpredictors.gff.gz|stdin : contains tandemid gene group and others -refgroup NCBI_GNO : reference prediction group for source x other tests opts: -skipids=idlist ; -pctoverlap=75 ; -pctinside=50 ; -dumpgenes=any|found -neargroup=9000 -neargroup=15000 : limit range of tandem distance -mrnatype=mRNA -exontype=CDS : what gff types to collect -skipchr='^U' : skip gff reference/chromosome pattern " unless($optok && ($tandemidpairs) && $sourcegff); # or $blastptable && $othergff # OR -blastptable dsec-pred6.blastp.gz : blastp table of matches # -eval 1e-5 or -bitscore 100 # -other allpredictors.gff.gz : other gene preds to compare to # needs input fix for cluster size among all pairs, not just protnear15k.ids .. # -cluster=10 -cluster=80 : min,max dupl source genes # ?? input blastp table and find pairs from that? $pctover= $pctover/100.0 if($pctover); $pctinside= ($pctinside) ? $pctinside/100.0 : 0.5; ## must have if(@clusize) { $clustersize= $clusize[0]; $clustertop= $clusize[1] || 0; } if($skipgrp) { $skipgrp =~ s/,/\|/g; } #? space also if($skipchr) { $skipchr =~ s/,/\|/g; } # my @eqclass= ($EQSAME, $EQINSIDE, $EQOUTSIDE, $EQOVERLAP, $EQPARTOVER, $EQNONE); # no Near, Far yet # ## (map{$EQNEAR * $_}(1..$NEARSTEPS)) , $EQFAR); # same, near, far # my @eqnames= ("Same", "Inside", "Outside", "Overlap", "Partov", "None"); ## (map{ my $kb= int((1 + $NEARDIST * $_)/1000); "Near".$kb."k"; } (1..$NEARSTEPS)) ## ,"Far"); # same, near, far # is None == SCMISSED and/or TDMISSED ? my @eqclass= ($EQSAME, $EQSCMISSED, $EQDOUBLED, $EQSKIPOVER, $EQALTTR, $EQOVERLAP, $EQNONE, $EQTDMISSED ); my @eqnames= ("Same", "NoSrc", "Dubl", "Skipov", "Alttr", "Ovlap", "None", "NoDup"); #, "NoDup" ## new gene model error classes: # - DOUBLED: gm contains 2+ gene-exon sets, double size of tandem dupl; # protein match has part matches to other full genes at low/high parts # - SKIPOVER: gm has 1 gene-exon set, but inside of gene span are 1+ other # ~complete gene-exon sets # - TDMISSED: no gm where other predictors find tandem gene (not quite same as skipover) my %eqnames; for (my $i=0; $i<@eqclass; $i++) { $eqnames{ $eqclass[$i] }= $eqnames[$i]; } # want tandem id pairs here instead ? see below my %skipids=(); my $keepids = 0; if($skipids) { $keepids = ($skipids =~ s/^[\+\^]//) ? 1 : 0; open(F,$skipids); while(){chomp; s/\s.*$//; $skipids{$_}++ if(/^\w/);} close(F); } #........ source duplicate pairs IN ............ my $nidpairs=0; my %pairs=(); # source gene pairs; combine with opairs # my %opairs=(); # other predictor pairs my %xpairs=(); # cross-predictor prot matches .. for err values my %altpairs=(); # detected as alt-tr; removed from pairs inputTandemidpairs($tandemidpairs); # sets %pairs, %xpairs # find real refgroup from pairs match to option unless( $refgroup and $pairs{$refgroup} ) { my @somegroup= sort keys %pairs; my ($rgroup) = grep /$refgroup/, @somegroup; $refgroup= ($rgroup) ? $rgroup : $somegroup[0]; } # drop this one # inputBlastpTable($blastptable, \%pairs ) if ($blastptable); #........ end source duplicate pairs IN ............ #........ source gff IN ............ %genespan=(); %ggroup=(); # this is main result: other-genes x ref genes overlap types my ($nsourcin, $allgrouplocs) = gffin( $sourcegff ); warn "# $refgroup: idpairs=$nidpairs; ngenes=$nsourcin; poor_ids=$npoor; skipped_ids=$nskip; \n" if($debug); my @allgroup= sort keys %genespan; ##my $refgrouplocs= {}; # makeRefgroupLocs($refgroup); ##my $allgrouplocs= makeAllgroupLocs(\@allgroup); #? do in gffin() and use to filter alt-tr? print_options( $refgroup, \@allgroup); my @sgene= sort keys %{$pairs{$refgroup}}; # test only tandem genes from id pairs list ## require refgroup genes be in genespan @sgene= grep { $genespan{$refgroup}{$_} } @sgene; my $testtan=0; foreach my $ogroup (@allgroup) { next if($ogroup eq $refgroup); # or test it?? $testtan= 0; # results in %ggroup; need only 1 group now not all %ggroup=(); # this is main result: other-genes x ref genes overlap types foreach my $sgene (@sgene) { my @standems= sort keys %{$pairs{$refgroup}{$sgene}}; my ($nosgene, $nmatch,$ndupl)= test_srctan_exists( $refgroup, $sgene, \@standems, $allgrouplocs, $ogroup); test_tandems( $refgroup, $sgene, \@standems, $allgrouplocs, $ogroup) unless($nosgene); } print_group( $ogroup, \%ggroup, $refgroup, \@sgene); ## results now ... } #.......... main subs ...................... sub makeAllgroupLocs { # my $allgrouplocs= makeAllgroupLocs(\@allgroup); my($allgroup)= @_; my %locs=(); my($ngene,$ngroup)=(0,0); foreach my $ggroup (@$allgroup) { $ngroup++; foreach my $gid (sort keys %{$genespan{$ggroup}}) { my ($tref,$tsource,$ttype,$tb,$te,$tor,)= split "\t",$genespan{$ggroup}{$gid}; my @bins= (int($tb/BINSIZE) .. int($te/BINSIZE)); foreach my $ib (@bins) { push @{$locs{$ggroup}{$tref}{$ib}}, $gid; } $ngene++; } } warn "# allgrouplocs ngroup=$ngroup, ngene=$ngene\n" if $debug>1; return \%locs; } # per predictor patches # PACH_GMP: only CDS, exon > construct mRNA; parent gene == Name=CGnnnn-Px; drop -Px my %pachsav=(); sub gffin_PACH_GMP { my($v)=@_; my @morev=(); if($v->[2] eq "CDS") { my($id)= $v->[8] =~ m/Parent=(\w+)/; my($gn)= $v->[8] =~ m/Name=(\w+)/; my($tb,$te)= ($v->[3],$v->[4]); if($pachsav{gn} && $pachsav{gn} ne $gn) { $pachsav{tab}->[8]= "ID=".$pachsav{id}.";Parent=".$pachsav{gn}; $pachsav{tab}->[2]= "mRNA"; # $mrnatype @morev=( $pachsav{tab}, $v ); # return last gene + curr CDS } $pachsav{tab}= [ @$v ]; $pachsav{gn}= $gn; $pachsav{id}= $id; $pachsav{tb}=$tb unless($pachsav{$id}{tb} && $tb > $pachsav{tb}); $pachsav{te}=$te unless($pachsav{$id}{te} && $te < $pachsav{te}); } return @morev; } # BREN_NSC: gene > mRNA sub gffin_BREN_NSC { my($v)=@_; $v->[2]= "mRNA" if($v->[2] eq "gene"); # $mrnatype ?? } =item gffin( sourcegff) =cut use constant CHECK_ALTTR => 1; sub gffin { my($sourcegff)= @_; my $inh= *STDIN; my $ok=0; if($sourcegff=~/^(stdin|\-)$/i) { $inh= *STDIN; $ok=1; } elsif($sourcegff=~/\.gz$/) { $ok= open(IN,"gunzip -c $sourcegff |"); $inh= *IN; } else { $ok=open(IN,$sourcegff); $inh= *IN; } die "open $sourcegff" unless($ok); my %didgene=(); my %locbins=(); my $nsourcin=0; my $nexon=0; my $nalttr=0; while(<$inh>){ next unless(/^\w/); chomp; my @v=split "\t"; next if($skipchr and $v[0]=~m/$skipchr/); # dang Dpse "U" vs "Unknown_group" ... next if($skipgrp and $v[1]=~m/$skipgrp/); # OR/must? use $ggroup below from gid my @allv= (\@v); if($v[1] eq "BREN_NSC") { gffin_BREN_NSC(\@v); } elsif($v[1] eq "PACH_GMP") { my @morev= gffin_PACH_GMP(\@v); @allv= @morev if(@morev); } # need to collect CDS & add new mRNA @v foreach my $vref (@allv) { @v= @$vref; my($tref,$fsource,$ftype,$tb,$te,$eval,$tor,$attr) = @v[0,1,2,3,4,5,6,8]; my ($gid,$pid); if($attr =~ m/ID=([^;]+)/) { $gid=$1; } if($attr =~ m/Parent=([^;]+)/) { $pid=$1; $gid= $pid unless($gid);} next unless($gid); # count? my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); $tor="" if(!$tor || $tor eq "."); if(%skipids){ if($keepids) { $nskip++ and next unless($skipids{$gid}); } else { $nskip++ and next if($skipids{$gid}); # gene/transcript id here; or both? (exonid) } } #? add ggroup here, as hash key? # ** need to weed out alt-tr somewhere; here? ncbi Parent gene id on mRNA will help that one if($ftype =~ m/$mrnatype/) { my $isalttr= 0; my @bins= (int($tb/BINSIZE) .. int($te/BINSIZE)); if(CHECK_ALTTR) { if($pid and $pid ne $gid) { # remove alttr; keep 1st found? is it longest? # also check pairs; do we have last, this mrna id? if($didgene{$pid}) { my @mrnaid= split ",", $didgene{$pid}; map { $isalttr++ if $pairs{$ggroup}{$_}; } @mrnaid; if($isalttr && $pairs{$ggroup}{$gid}) { $altpairs{$ggroup}{$gid}= delete $pairs{$ggroup}{$gid}; } } $nalttr++ and next if $isalttr; $didgene{$pid} .= $gid.","; } } ##? want this or not for other group check? if(CHECK_ALTTR) { foreach my $ib (@bins) { next unless( $locbins{$ggroup}{$tref}{$ib} ); foreach my $ag ( @{$locbins{$ggroup}{$tref}{$ib}}) { next unless ($pairs{$ggroup}{$ag}); my($aref,$asource,$atype,$ab,$ae,$aor)= split( "\t", $genespan{$ggroup}{$ag}); if(_isoverlap($ab,$ae,$tb,$te)) { $isalttr++; last; } } last if($isalttr); } if($isalttr && $pairs{$ggroup}{$gid}) { $altpairs{$ggroup}{$gid}= delete $pairs{$ggroup}{$gid}; } $nalttr++ and next if $isalttr; } $nsourcin++; $genespan{$ggroup}{$gid} = join("\t",$tref,$fsource,$ftype,$tb,$te,$tor,$gid)."\n"; foreach my $ib (@bins) { push @{$locbins{$ggroup}{$tref}{$ib}}, $gid; } } elsif($ftype =~ m/$exontype/) { # $useexons and $exons{$gid} .= "$tb\t$te\n"; $nexon++; } } # @allv } close($inh); warn "# gffin genes=$nsourcin, exons=$nexon, alt_tr=$nalttr\n" if $debug; return ($nsourcin,\%locbins); } =item testgene( $refgroup, $refgrouplocs, $ogroup, $ogene) for each target gene, classify duplicate status to ref duplicate gene set: same, source missed, tandems (of ref-source gene) missed, tandem/source doubled, tandem/source skipover global in: refgrouplocs{binned}, genespan{refid} global out: sets ggroup{group}{tid}{refid} = eqvalue =item test_tandems( $refgroup, $sgene, \@standems, $allgrouplocs, $ogroup); foreach sgene in pairs of source genes: 1. sgene tests : any other gene matches? (same, overlap) == SCMISSED -- is ogene overlap a DOUBLED gene? a SKIPOVER : need all sgene,tandems foreach tandems of sgene 2. @tgene tests: any other gene match (same, overlap == TDMISSED) 3. DOUBLED, SKIPOVER tests across sgene, tgenes ## ** Need to change this to start from ref gene list, class matches to all other genes =cut sub get_genelocs { my($ggroup, $gene)= @_; return () unless ($genespan{$ggroup}{$gene}); my($gref,$gsource,$gtype,$gb,$ge,$gor)= split( "\t", $genespan{$ggroup}{$gene}); my @gexons= split( "\n", $exons{$gene}||""); my @gexonv= map{ my @xbe=split "\t"; \@xbe; } @gexons; if($gor eq "-") { # got missing or vals? @gexonv= sort{ $b->[0] <=> $a->[0] } @gexonv; } else { @gexonv= sort{ $a->[0] <=> $b->[0] } @gexonv; } return ($gref,$gsource,$gtype,$gb,$ge,$gor,\@gexonv); } sub test_srctan_exists { my( $refgroup, $sgene, $standems, $allgrouplocs, $ogroup)= @_; my $ndupl=0; my $nmatch=0; my $nosgene= 0; foreach my $gn ($sgene, @$standems) { my($gref,$gsource,$gtype,$gb,$ge,$gor,$gexonv)= get_genelocs($refgroup,$gn); return (-1,0,0) unless($ge or $ndupl>0); next unless($ge); my @sabins= (int($gb/BINSIZE) .. int($ge/BINSIZE)); my $matchgn= 0; BINL: foreach my $ib (@sabins) { next unless $allgrouplocs->{$ogroup}{$gref}{$ib}; foreach my $ogene (@{$allgrouplocs->{$ogroup}{$gref}{$ib} }) { my($tref,$tsource,$ttype,$tb,$te,$tor,$otexonv)= get_genelocs($ogroup,$ogene); if( _isoverlap($gb,$ge,$tb,$te)) { $matchgn++; last BINL; } } } $ndupl++; $nmatch++ if($matchgn); if($ndupl == 1 && $nmatch==0) { $nosgene=1; $ggroup{$ogroup}{source}{$sgene}{none}= $EQSCMISSED; } } if($ndupl > 1 && $nosgene==0 && $nmatch < $ndupl) { $ggroup{$ogroup}{source}{$sgene}{duplmissed}= $EQTDMISSED; } return ($nosgene, $nmatch, $ndupl); } sub test_tandems { my( $refgroup, $sgene, $standems, $allgrouplocs, $ogroup)= @_; my $debugtest= 0; my($gref,$gsource,$gtype,$gb,$ge,$gor,$gexonv)= get_genelocs($refgroup,$sgene); return -1 unless($ge and @$gexonv); # require exons my %tdgenes=(); foreach my $tgene (@$standems) { my @tgenev= get_genelocs($refgroup,$tgene); $tdgenes{$tgene}= \@tgenev; # if(@tgenev); # save even if empty; see below } my %oids=(); my @sabins= (int($gb/BINSIZE) .. int($ge/BINSIZE)); foreach my $ib (@sabins) { next unless $allgrouplocs->{$ogroup}{$gref}{$ib}; map{ $oids{$_}=1; } @{$allgrouplocs->{$ogroup}{$gref}{$ib} }; } my $ntandem= @$standems; my $noids= scalar(keys %oids); $debugtest= $debug>1 && $testtan++ < 200 && $testtan > 100; # && (($testtan % 10) == 4); warn "# test_tandems[$testtan] ref=$refgroup, sgene=$sgene, ntandem=$ntandem, nother=$noids\n" if $debugtest; my %didloc=(); my $olap=0; my $eq= $EQNONE; my @ogene= (sort keys %oids); # ## want to distinguish EQSCMISSED and EQTDMISSED : each gene is both source and tandem here # ## so src missed contains tandem missed: separate case of all missed in family vs 1+ found # ## see test_srctan_exists # if(@ogene == 0) { # $eq= $EQSCMISSED; # $EQTDMISSED, ; # $ggroup{$ogroup}{source}{$sgene}{none}= $eq; # #? $ggroup{$ogroup}{other}{none}{$sgene}= $eq; # } foreach my $ogene (@ogene) { next unless $genespan{$ogroup}{$ogene}; # error; count? next if ($didloc{$sgene.$ogene}++); my($tref,$tsource,$ttype,$tb,$te,$tor,$otexonv)= get_genelocs($ogroup,$ogene); next unless($otexonv and @$otexonv); # require exons # next unless($gref eq $tref); # dont need; above select only for gref == tref next if($gor && $tor && $gor ne $tor); next unless(_isoverlap($gb,$ge,$tb,$te)); my @tdover= (); # this is only for cases doubled, skipover below? foreach my $tgene (@$standems) { my($sref,$ssource,$stype,$sb,$se,$sor)= @ { $tdgenes{$tgene} }; if ($sref && $sref eq $tref && $sor eq $tor && _isoverlap($sb,$se,$tb,$te)) { push(@tdover,$tgene); } } my $oerr= 0; my $blerr= $xpairs{$ogroup}{$ogene}{$sgene} || 0; $oerr = $blerr if($blerr && $blerr != $EQSAME); $olap= 1; $eq= $EQOVERLAP; # change these inside, partover tests: use exons, count # 3a. DOUBLED : 2+ source genes contained in exon set of target # _same_exons subset? >> want also protein eq info for 3a. # .. want to see all,most of 2+ sgid exons matching otexonv # 3b. SKIPOVER : 2+ source genes inside target exon set but only 1 set used # .. want to see some exons of 2+ sgid matching @otexonv # also count 2,4 SCMISSSED, TDMISSED # .. cant count source or source-tandems missed here with iterate over other genes # full vs part overlap? with exon matching this is less useful ## some, many -3 missing exons are alt-tr like diffs; should class that way? my($issame, $samex, $ngex, $ntex)= _same_exons( $gexonv, $otexonv, $gor); my $nsamex= (ref $samex) ? scalar(@$samex) : 0; # my $issame= _same_exons( $gexonv, $otexonv, $gor); # issame : 1 == all same; -1 == all gexon ; -2 == all texon ; -3 some same; 0 == none my $ndouble=0; my $nskip=0; my $naltr=0; if($issame == 1) { $eq= $EQSAME; # keep but need exons test here } elsif($issame == 0) { # * genespan overlap but no same exons; test exon overlaps? # got lots of these: is it exon model disagreement? e.g. > 10 bases # * could be skipover of this gene; test tandems #was#? $eq= $EQNONE; # or overlap? better none if no shared exons $eq= $EQOVERLAP; foreach my $tdgene (@tdover) { my $tdexonv= $tdgenes{$tdgene}->[-1]; my $tdsame = _same_exons( $tdexonv, $otexonv, $tor); next if($tdsame==0); ## getting lots of NCBI alt-tr here for dmel-like; one alt misses gleanr; other matches $nskip++ if($tdsame==-3); # likely skipover $nskip++ if($tdsame==-2); $nskip++ if($tdsame==-1); # #$nskip++ if($tdsame==1); # is this possible? } if($nskip>0) { $eq= $EQSKIPOVER ; } # ?? ^^ is this suspect call; not likely if ?? $ngex < $ntex and/or $oerr == 0 elsif($oerr) { $eq= $oerr; } # is this valid? from blastnear err class } elsif ($issame == -1) { # all gexon match , but more texon > doubled ?? foreach my $tdgene (@tdover) { my $tdexonv= $tdgenes{$tdgene}->[-1]; my $tdsame = _same_exons( $tdexonv, $otexonv, $tor); next if($tdsame==0); $ndouble++ if($tdsame==-1); # likely doubled; need what more ? $ndouble++ if($tdsame==-3); # likely doubled; need what more ? # $ndouble++ if($tdsame==1); # shouldnt be possible # $ndouble++ if($tdsame==-2); # shouldnt be possible } # ** this test isnt working; really need blastp refgene x othergene result: errD for this if($ndouble>0) { $eq= $EQDOUBLED; } elsif($oerr) { $eq= $oerr; } # is this valid? from blastnear err class elsif($ngex>2 && $nsamex / $ngex >= 0.75) { $eq= $EQALTTR; } # another overlap category? } elsif ($issame == -2) { # all texon match; partial gene model? or g is doubled? $eq= $EQOVERLAP; # should be $EQINSIDE; use this for now $EQOVERLAP; # or call it $EQALTTR ? Overlap and Alttr groups not distinct w/o other tests if($oerr) { $eq= $oerr; } # is this valid? from blastnear err class elsif($ngex>2 && $nsamex / $ngex >= 0.75) { $eq= $EQALTTR; } # another overlap category? } elsif ($issame == -3) { # some same : skipover?, maybe ~doubled # other class: alt-tr like diff, 1 exon of several differs foreach my $tdgene (@tdover) { my $tdexonv= $tdgenes{$tdgene}->[-1]; my $tdsame = _same_exons( $tdexonv, $otexonv, $tor); next if($tdsame==0); $ndouble++ if($tdsame==-1); # likely doubled; need what more ? $nskip++ if($tdsame==-3); # likely skipover; need what more ? # $ndouble++ if($tdsame==1); # shouldnt be possible # $ndouble++ if($tdsame==-2); # shouldnt be possible } # should count total gene + tandem exon matches to other gene to distinguish these two if($ndouble>0) { $eq= $EQDOUBLED ; } elsif($nskip>0) { $eq= $EQSKIPOVER ; } elsif($oerr) { $eq= $oerr; } # is this valid? from blastnear err class elsif($ngex>2 && $nsamex / $ngex >= 0.75) { $eq= $EQALTTR; } # another overlap category? } if($eq != $EQSAME || $blerr != 0) { $testtan++; $debugtest= $debug>1 ? 1 : 0; #?? && $testtan < 200; # && $testtan > 200 && (($testtan % 10) == 3); } warn "# test ogene=$ogene, eq=$eq, issame=$issame, oblerr=$blerr, ngex=$ngex, noex=$ntex\n" if $debugtest; $ggroup{$ogroup}{source}{$sgene}{$ogene}= $eq; #? is this refgroup or $ogroup $ggroup{$ogroup}{other}{$ogene}{$sgene}= $eq; } } ## FIXME .. sub filterPairs { # extend here to all groups ... ## next unless($ggroup =~ m/$refgroup/); ## $nskip++ and next unless( exists $pairs{$refgroup}{$gid} ); # ignore other source genes ? # here or gffin: find alt-tr if possible and remove from pairs (which? shorter?) ## clustersize is not good with current input id pairs: limited to 15k near tandems ## need all dupl. prots to count cluster sizes, then restrict by distance. # if($clustersize) { # my $ndups= scalar( keys %{ $pairs{$refgroup}{$gid} } ); # my $skipit=0; # $skipit=1 if($clustersize>0 and $ndups < $clustersize); # $skipit=1 if($clustersize<0 and $ndups > -$clustersize); # $skipit=1 if($clustertop>0 and $ndups > $clustertop); # if($skipit) { # delete $pairs{$refgroup}{$gid}; # $nskip++ and next; # } # } } sub inputTandemidpairs { my($tandemidpairs)=@_; ## $pairs={} unless(ref $pairs); my $ok= ($tandemidpairs=~/\.gz$/) ? open(IN,"gunzip -c $tandemidpairs |") : open(IN,$tandemidpairs); die "open $tandemidpairs" unless($ok); my $ncrosspairs=0; my $altskip=0; # fixme for refgroup defined and idlist has many others while(){ if(m/^(\w+)[\,\.\s](\w+)/) { my($ga,$gb)=($1,$2); my($distance,$modelerr)= (0,""); $distance = $1 if (m/\t(\d+)/); $modelerr = $1 if (m/\,(err\w+)/); #???? # check; add err classes from blastnear: errDn == doubled; errOn == skipover # new pair scores from blastnear results: num is kb distance for near (0>num<15kb) # dist <=0 : num == 199 << skip these, alt tr # dist far : num == 299 # dist unlink: num == 399 my $modval= $EQSAME; if($modelerr =~ m/errD/) { $modval= $EQDOUBLED; } # dang ; this is ref x ref score; want ref x other elsif($modelerr =~ m/errO/) { $modval= $EQSKIPOVER; } # dang ; this is ref x ref score; want ref x other # if( m/,errD/) { $modval= $EQDOUBLED; } # elsif( m/,errO/) { $modval= $EQSKIPOVER; } my($agroup)= $ga =~ m/^(\D+)/; my($bgroup)= $gb =~ m/^(\D+)/; unless($refgroup) { ($refgroup)= ($agroup); # $ga =~ m/^(\D+)/; } # ##dont need here - see above after this call ; need refgroup to match # elsif ($agroup =~ /$refgroup/ and $agroup ne $refgroup) { # $refgroup= $agroup; # } $altskip++ and next if($distance == 199 && $agroup eq $bgroup); if($agroup eq $bgroup) { $pairs{$agroup}{$ga}{$gb}= $modval; $nidpairs++; } else { $xpairs{$agroup}{$ga}{$gb}= $modval; $ncrosspairs++;} # error is for ga, based on gb match } } close(IN); warn "# read pairs=$nidpairs; crossmatch=$ncrosspairs; altskip=$altskip\n" if $debug; return $nidpairs; ## ($pairs,$nidpairs); #$opairs, } sub print_options { my $ttype = "gene"; ## ($BY_GENE) ? "gene" : "exon"; my $cc = $dumpgenes ? "#t " : ""; my @flags=(); # push(@flags,reference_group => $refgroup) if($refgroup); # push(@flags,keep_only => $keeponly) if($keeponly); # push(@flags,overlap_filter => $overmarkfilter) if($overmarkfilter); # push(@flags,count_by => (($COUNTKEYS)? "ids" : "matches")); # push(@flags,min_align => $MINALIGN); # push(@flags,min_e_value => $MINEVAL); push(@flags,pct_overlap => $pctover) if($pctover); push(@flags,pct_inside => $pctinside) if($pctinside); push(@flags,tandem_pairs => $tandemidpairs); push(@flags,source_genes => $sourcegff); ## push(@flags,target_genes => $othergff); push(@flags,clusters => "$clustersize..$clustertop") if($clustersize); push(@flags,clusters => "all") unless($clustersize); push(@flags,distance_limit => "$neargroup[0]..$neargroup[1]") if(@neargroup); my %flags=@flags; my $flags= join ",", map{"$_=".$flags{$_}} sort keys %flags; print $cc,"Tandy-prot count of $ttype per predictor group\n"; print $cc,"Options: $flags\n"; print $cc,"Source : $refgroup\n"; print $cc,join("\t","Group ", "Stat ",@eqnames),"\n"; } sub print_group { # my($ggroup)= @_; my( $ggroup, $ggroupres, $refgroup, $sgenea)= @_; # my %ggroup= %$ggroupres; # global .. my $cc = $dumpgenes ? "#t " : ""; my ($pc,$pv,$tsame); # foreach my $ggroup (@ggroup) { my $sgroup= sprintf("%-15s", substr($ggroup,0,14)); my %ofound=(); # other found count classes? my %odone=(); my %tandone=(); # my @sgenes= sort keys %{$pairs{$refgroup}}; # grep pairs{a}{b}>0 or >1 ? my @sgenes= @$sgenea; my $nsourcegenes= 0; # scalar(@sgenes); my $ntandemgenes= 0; my $nalttr_ref = $altpairs{$refgroup} ? scalar( keys %{ $altpairs{$refgroup} } ) : 0; my $nalttr_other= $altpairs{$ggroup} ? scalar( keys %{ $altpairs{$ggroup} } ) : 0; foreach my $sgene (@sgenes) { next if($refgroup && $sgene !~ m/$refgroup/); $nsourcegenes++; my @standems= sort keys %{ $pairs{$refgroup}{$sgene} }; # recall we recip mapped these: dup count my @others = sort keys %{ $ggroup{$ggroup}{source}{$sgene} }; #.. includes {none} no srcmatch; add {duplmissed}? EQTDMISSED ? ## can we count sgenes,standems w/ no other overlaps ? == EQNONE ## always have others keys here, but below standems may be missing others foreach my $oid (@others) { next if($odone{same}{$sgene.$oid}++); my $eq= $ggroup{$ggroup}{source}{$sgene}{$oid} || $EQNONE; $ofound{same}{$eq}++; } foreach my $tid (@standems) { next if($refgroup && $tid !~ m/$refgroup/); $ntandemgenes++ unless($tandone{$tid}++); @others= sort keys %{ $ggroup{$ggroup}{source}{$tid} }; push(@others, "Missing") unless @others; foreach my $oid (@others) { next if($odone{tandem}{$tid.$oid}++); my $eq= $ggroup{$ggroup}{source}{$tid}{$oid} || $EQNONE; $ofound{tandem}{$eq}++; } } } my @ogenes= sort keys %{ $ggroup{$ggroup}{other} }; my $nothergenes= scalar(@ogenes); foreach my $oid (@ogenes) { my @sogenes = sort keys %{ $ggroup{$ggroup}{other}{$oid} }; ## here want to know how many standems an other covers ## total count wont do ? foreach my $sgene2 (@sogenes) { next if($odone{other}{$oid.$sgene2}++); my $eq= $ggroup{$ggroup}{other}{$oid}{$sgene2} || $EQNONE; $ofound{other}{$eq}++; } } my @same = map { $ofound{same}{$_} || 0; } @eqclass; my @tandem= map { $ofound{tandem}{$_} || 0; } @eqclass; my @other = map { $ofound{other}{$_} || 0; } @eqclass; ## ** with %odone correction above for already counted pairs, ## source>other == tandem>other but for extra Missing/None case of tandem>other print $cc,join("\t",$sgroup,"source.genes",$nsourcegenes, "source.alttr",$nalttr_ref),"\n"; print $cc,join("\t",$sgroup,"tandem.genes",$ntandemgenes),"\n"; print $cc,join("\t",$sgroup,"other.genes", $nothergenes,"other.alttr",$nalttr_other),"\n"; print $cc,join("\t",$sgroup,"source>other",@same),"\n"; print $cc,join("\t",$sgroup,"tandem>other",@tandem),"\n"; #fixme.now == same # print $cc,join("\t",$sgroup,"other>source",@other),"\n"; print "\n"; #dumpgenes($ggroup,\@sgenes,) if($dumpgenes); dumpgenes($ggroup, $sgenea, $refgroup) if($dumpgenes); } } sub dumpgenes { my($ggroup,$sgenes,$refgroup)= @_; my %ofound; my $cc=""; my $ig=0; print "#i other genes overlapping source+tandem genes \n"; print "#i", join("\t","Class", "Gene", "Location", "Other", "OverlapType", ),"\n"; my $noduplab=",".$eqnames{$EQTDMISSED}; # ggroup{$ogroup}{source}{$sgene}{duplmissed}= $EQTDMISSED; foreach my $sgene (@$sgenes) { next if($refgroup && $sgene !~ m/$refgroup/); $ig++; my @standems= sort keys %{ $pairs{$refgroup}{$sgene} }; # recall we recip mapped these: dup count my @others = grep !/duplmissed/, sort keys %{ $ggroup{$ggroup}{source}{$sgene} }; my $nodupv = $ggroup{$ggroup}{source}{$sgene}{duplmissed} ? $noduplab : ""; my($gref,$gsource,$gtype,$gb,$ge,$gor,$xtrasgid)=("") x 10; ($gref,$gsource,$gtype,$gb,$ge,$gor,$xtrasgid)= split "\t",$genespan{$refgroup}{$sgene} if $genespan{$refgroup}{$sgene}; push(@others,"Missing") unless(@others); # show no others... foreach my $oid (@others) { # next if($oid =~ /^duplmissed/); my $eq= $ggroup{$ggroup}{source}{$sgene}{$oid} || $EQNONE; print $cc, join("\t","source.$ig", $sgene, "$gref:$gb-$ge:$gor", $oid, $eqnames{$eq}.$nodupv, ),"\n"; } my $it=0; foreach my $tid (@standems) { next if($refgroup && $tid !~ m/$refgroup/); $it++; my @tothers= grep !/duplmissed/, sort keys %{ $ggroup{$ggroup}{source}{$tid} }; my $nodupv = $ggroup{$ggroup}{source}{$tid}{duplmissed} ? $noduplab : ""; push(@tothers,"Missing") unless(@tothers); # show no others... my ($tref,$tsource,$ttype,$tb,$te,$tor,)= ("") x 10; ($tref,$tsource,$ttype,$tb,$te,$tor,)= split "\t",$genespan{$refgroup}{$tid} if $genespan{$refgroup}{$tid}; foreach my $oid (@tothers) { # next if($oid =~ /^duplmissed/); my $eq= $ggroup{$ggroup}{source}{$tid}{$oid} || $EQNONE; print $cc, join("\t","tandem.$ig.$it", $tid, "$tref:$tb-$te:$tor", $oid, $eqnames{$eq}.$nodupv, ),"\n"; } } } print "\n"; } #....... sub _min { return ($_[1] < $_[0]) ? $_[1] : $_[0]; } sub _max { return ($_[1] > $_[0]) ? $_[1] : $_[0]; } sub _isoverlap { my($tb,$te, $qb,$qe)= @_; # return ($tb <= $qe && $te >= $qb) ? 1 : 0; my $over= ($tb <= $qe && $te >= $qb) ? 1 : 0; 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 (1,$rev) unless($over); } return $over; } sub _isoverlapfull { my($tb,$te,$to, $qb,$qe,$qo)= @_; if($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; return ($maxo/$leno < $pctover) ? 0 : 1; } return 1; } sub _sameloc { my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_; return ( $aor eq $bor && # strand/orient (abs($ab-$bb) < SAMEBASE) && (abs($ae-$be) < SAMEBASE) ) ? 1 : 0; } # sub _same_exons_OLD { # my( $aex, $bex)= @_; # # aex, bex are exon list [$b,$e] # # we are assuming exons are sorted by loc # ## this is bad for case of longer lead exon set for one gene # ## need to step thru all to look for alignment of exons # my $nax= @$aex; my $nbx= @$bex; # my $n= _min($nax,$nbx); # my @xs=(); # for(my $i=0; $i<$n; $i++) { # note these may be reverse sorted, if gene is rev # # my($ab,$ae)= @{$aex->[$i]};# ?? is this right # # my($bb,$be)= @{$bex->[$i]}; # my($ab,$ae)= ($aex->[$i]->[0], $aex->[$i]->[1]); # my($bb,$be)= ($bex->[$i]->[0], $bex->[$i]->[1]); # ##push(@xs,$i) if(abs($ab-$bb) <= SAMEBASE and abs($be-$be) <= SAMEBASE); # BUG be - be ** # push(@xs,$i) if(abs($ab-$bb) <= SAMEBASE and abs($ae-$be) <= SAMEBASE); # BUG be - be ** # } # my $issame= (@xs == 0) ? 0 # : ($nax == $nbx && @xs == $nax) ? 1 # : (@xs == $nax) ? -1 # : (@xs == $nbx) ? -2 # : -3; # some same # return wantarray ? ($issame, \@xs, $nax, $nbx, ) : $issame; # } sub _same_exons { my( $aex, $bex, $orient)= @_; # aex, bex are exon list [$b,$e] # exons are sorted by loc ** REVERSE sorted when -strand my $nax= @$aex; my $nbx= @$bex; my $n= _min($nax,$nbx); return 0 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=($orient eq "-")?1:0; # test data instead 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) ); } } $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]); push(@xs,$i) if(abs($ab-$bb) <= SAMEBASE and abs($ae-$be) <= SAMEBASE); } my $issame= (@xs == 0) ? 0 : ($nax == $nbx && @xs == $nax) ? 1 : (@xs == $nax) ? -1 # all of a : (@xs == $nbx) ? -2 # all of b : -3; # some same return wantarray ? ($issame, \@xs, $nax, $nbx, ) : $issame; } sub _isinside { #assume overlap? check for g-span >> q-span my($tb,$te,$to, $qb,$qe,$qo)= @_; if($qb > $tb and $qe < $te) { my $tlen= abs($te - $tb) || 1; my $qlen= abs($qe - $qb); return 1 if( $qlen/$tlen < $pctinside); # 0.50 # 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 #? what of _isinside ? my $over= ($tb <= $qe && $te >= $qb) ? 1 : 0; 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 (1,$rev) unless($over); } 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); } } __END__