# tandem5.pm # package main; =item trial package version for tandem gene finding ; see main perl tandemgenes.pl =cut =item usage ./tandemgenes.pl -debug -act tandy \ -genome scaffold_4/dpulex1.fa -que scaffold_4/dpulex1_exons.nr \ -mineval 1e-30 -minalign 0.50 -loc=scaffold_4:830000-1030000 -igpartial > cyp934f2.tandy grep tandy.common cyp934f2.tandy | perl -ne\ 's/DP_D\D+/D/g; s/NCBI\D+/G/g; s/Dap\D+/J/g; @v=split; print join("\t",@v[1,3,8]),"\n";' =item NOTE on results The main result is identification of adjacent duplicated and conserved (multi-)coding-exon regions. Tandem gene calls can, likely do, include pseudogenes and other mixed models among this. =item single exon genes are a problem Many of the tests and grouping filters here are looking at multi-exon runs, so single exon genes are not handled well. This is a problem particularly for fruitfly data, where 1-exon genes are about 50% of genome, versus 10% or less for mouse, worm, daphnia. =item tandem algo4 1. find interesting ++ marked exon/gene 2. collect all exons of geneid from ++marked in region $begin - 10,000 , $end +10,000 save as geneid->@exons; this should include exons of tandem dupls. 3. collect all geneid exons from geneids/altids of 2 in same region 4. segregate exons to gene models in tandem region (how??) - using exon.nums from ids, break at change in exon.num sequence? 1,2,3,4,5 | 2,3,4,5 - using max intron size? max/ave of step b/n base exons? 5. mark as done all gene/exons in region and move on (step 1) =item genegroup ing this is better, but not quite good enough. -- have gene groups that overlap regions (probably exons also?) .. but ignore genes inside group region that dont share exons -- missing exons of normal/eq=equal gene calls in e.g. tandem gene groups add: 1. hash up exons by geneid and add to this mix any geneid exons that fall inside group region (but not outside) 2. scan all of group region for exons overlapping found set (but not non-overlap-tandem-exons) =item FIXME: merge_genegroup needs to create new geneids and need to add attributes better linking pairs/multigenes (good name would be ok, but also highlight altids that are pairs, distinguish from altids at same location). =item FIXME: flag where current predictions may be wrong added exons, changed gene models, etc. mark best/all gene model (if any) matching new-found tandem/dupl. models =item single exon genes are a problem many of the checks, grouping filters are looking at multi-exon runs, so single exon genes are not handled well. =item gff fixup; correct tclass in this package, using alt eq values; cat scaffold_?/dpulex1_exons.nr5o.gff | perl -pe'\ if(/^scaff/){ @v=split"\t"; $v[1]=~s/tandy.\w+/tandy/; \ if(/tclass=(equal|far)/ && $v[-1] =~ m/\=\-2/) { $nn++; $v[-1] =~ s/tclass=\w+/tclass=near/;} \ $_=join "\t",@v; } END{ warn "# Changed eq2near=$nn\n"; } ' > ! dpulex1_tandy5o.gff =item gene qual score parsing (match lines bestids=): grep bestids= scaffold_6541/dmoj_caf060210_exons_tandy.gff | perl -n genebest.perl grep bestids= scaf5q.tandy | perl -n genebest.perl # genebest.perl # usage: grep match scaffold.tandy.gff | perl -ne genebest.perl chomp; m/bestids=([^;]+)/ or next; $b=$1; @b=split",",$b; foreach (@b){ s,\[,/,; s/\]//; s/\ssltq//; ($d,@s)= split"/"; ($m=$d)=~s/\d//g; $ms{$m}{5} ++; foreach $k (0..2) { $ms{$m}{$k} += $s[$k]; } $ms{$m}{3} ++ if($s[3]<-1); $ms{$m}{4} ++ if($s[3]>0); } END{ print "# gene bestids\nscores/method\n"; print join"\t",qw(predictors num score lost tandem near equal),"\n"; foreach $m (sort keys %ms) { $n=$ms{$m}{5}; print"$m:\t$n"; foreach $k (0..4) { $v= sprintf("%.2f",$ms{$m}{$k}/$n); print"\t",$v; } print"\n"; } } # ............... grep bestids= scaffold_6541/dmoj_caf060210_exons_tandy.gff | perl -n $td/geneinfo.perl grep bestids= scaf5q.tandy | perl -n geneinfo.perl # geneinfo.perl # usage: grep match scaffold.tandy.gff | perl -ne geneinfo.perl chomp; m/info=([^;]+)/ or next; $b=$1; @b=split",",$b; %b= map { split ":" }@b; $tandemnew= $b{eqnear} > 1 && $b{eqsame} < 2; $tandemknow= !$tandemnew && $b{eqnear} > 1 && $b{fullgenes}>0; $bs{all}{n}++; map { $bs{all}{$_}+= $b{$_}} keys %b; $bs{all}{tandemnew} += $tandemnew; $bs{all}{tandemknow} += $tandemknow; ($bid)= m/bestids=([^;]+)/; %bid= map { s/\[.*//; s/\d//g; $_,1; } split",",$bid; @bid= sort keys %bid; foreach $bid (@bid) { $bs{$bid}{n}++; map{ $bs{$bid}{$_}+= $b{$_};} keys %b; $bs{$bid}{tandemnew} += $tandemnew; $bs{$bid}{tandemknow} += $tandemknow; } END{ @meth=sort keys %bs; @f= grep { $_ ne "n" } sort keys %{$bs{all}}; print "# gene info, all-methods\nscores/method\n"; print join"\t","predictor", "count",@f,"\n"; foreach $m (@meth) { $n=$bs{$m}{n}||1; printf "%-11s\t%3d",$m,$n; foreach $k (@f) { $v=$bs{$m}{$k}; $v= ($k =~ m/tandem(new|kno)/) ? $v : sprintf("%.2f",$v/$n); print"\t",$v; } print"\n"; } } # ............... ... interesting matches: $tandemnew = $b{eqnear} > 1 && $b{eqsame} < 2; $tandemknow= $b{eqnear} > 1 && $b{fullgenes}>0; ... only best id method $bs{all}{n}++; ($bid)= m/bestids=(\D+)/; $bs{$bid}{n}++; \ map{ $bs{$bid}{$_}+= $b{$_}; $bs{all}{$_}+= $b{$_}} keys %b; \ ... all id methods $bs{all}{n}++; map { $bs{all}{$_}+= $b{$_}} keys %b; \ ($bid)= m/bestids=([^;]+)/; %bid= map { s/\[.*//; s/\d//g; $_,1; } split",",$bid; \ @bid= sort keys %bid; \ foreach $bid (@bid) { $bs{$bid}{n}++; map{ $bs{$bid}{$_}+= $b{$_};} keys %b; } \ # ............... * gene bestids scores/method scaffold_4 predictors num score lost tandem near equal DP_DGIL_SNO_: 745 65.05 8.17 1.13 0.31 0.32 Dappu_FM_: 262 58.97 3.60 1.03 0.40 0.35 NCBI_GNO_: 561 62.45 4.09 1.67 0.42 0.32 # ............... # gene info, all-methods scores/method dpulex scaffold_4 predictor count eqfar eqnear eqsame fullgns gnids maxgns methods nexons tandkno tandnew tandems DP_DGIL_SNO_ 375 0.87 1.70 3.60 0.97 3.99 2.19 1.91 4.55 76 16 0.73 Dappu_FM_ 158 0.99 3.20 6.13 1.81 7.21 3.30 2.88 7.22 64 7 1.37 NCBI_GNO_ 245 0.85 2.43 4.98 1.36 5.40 2.60 2.45 5.91 70 13 1.04 all 421 0.81 1.58 3.32 0.91 3.72 2.09 1.85 4.24 78 19 0.70 =item data bug dang, looks like DP_DGIL_SNO_ CDS exons are in source dpulex1_predict.gff twice, and have diff exon ids for same location ... this throws off calcs here; shouldnt expect this problem from "good" input data. However such could be screened at =cut use strict; our ($VERSION)= "1.5q"; use constant SOURCE_HAS_EQTYPE => 1; # debug mode gff output use constant SAMEBASE1 => 10; # for _sameloc, slop allowed in loca == locb use constant COMMON_MERGE => 1; our $FULL_EXONLIST = 1; ## only for ($debug>1) output now our %EQlegend; our $GENEID_NUM; our $COMMONID_NUM; our %commonids=(); our %commonlist=(); our $IDmap={}; # useful package global our ($debug); our $NOTEOUT; # ($doprint) ? *STDOUT : *STDERR; our $IDPREFIX; our $REPEATSKIP; # user boolean # dang our, do we get main values here? our $BINSIZE ; our $GENEBINSIZE ; # use only where also checking gene ID our $NEARDIST; #? our $MIN_EXON_PRINT= 0; # gff SO types our $genegrouptype="region"; our $genetype = "match"; # region for common our $exontype = "HSP"; # "match_part"; # match_part or HSP for ggb ?? our $sourcetag = "tandy"; # modify with dataset tag ** our $sourcenear = "tandy.near"; our $sourcefar = "tandy.far"; #? change sourcetag or type for near/far difference? BEGIN{ $FULL_EXONLIST= 0 unless($debug>1); $IDPREFIX="td_" unless defined $IDPREFIX; # drop it? $NOTEOUT= *STDERR unless($NOTEOUT); ## warn "# FIXME: this perl wont yet handle multiple chromosome/scaffold input "; } sub _sameloc1 { # ref1,b1,e1 vs ref2,b2,e2; equal location with some +/- 10b slop my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_; $ar =~ s/[\_]+$//; $br =~ s/[\_]+$//; # $ar =~ s/_+$//; $br =~ s/_+$//; $aor ||= 0; $bor ||= 0; return # $ar eq $br && #?? is this bad? $aor eq $bor && (abs($ab-$bb) < SAMEBASE1) && (abs($ae-$be) < SAMEBASE1) ; } sub _isoverlap1 { my ($gb,$ge, $qb,$qe)= @_; if ($gb <= $qe && $ge >= $qb) { if(abs($ge-$gb) < abs($qe-$qb)) { return -1; } else { return 1; } } else { return 0; } } sub _isoverlapNotHuge1 { my ($gb,$ge, $qb,$qe)= @_; if (abs($ge-$gb) > 2*$NEARDIST || abs($qe-$qb) > 2*$NEARDIST) { return 0; } if ($gb <= $qe && $ge >= $qb) { if(abs($ge-$gb) < abs($qe-$qb)) { return -1; } else { return 1; } } else { return 0; } } sub _sameorlap1 { my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_; if ($ab <= $be && $ae >= $bb) { # $ar =~ s/[\_]+$//; $br =~ s/[\_]+$//; $aor ||= 0; $bor ||= 0; my $same= #? $aor eq $bor && (abs($ab-$bb) < SAMEBASE1) && (abs($ae-$be) < SAMEBASE1) ; return ($same) ? 1 : -1; # return 1 if($same); # return -1 if ($ab <= $be && $ae >= $bb); # overlap; return +2 instead? } return 0; } sub _isnear1 { # ($gb,$ge, $qb,$qe, $qid) my $gm= int(($_[0]+$_[1])/2); my $qm= int(($_[2]+$_[3])/2); return (abs($gm - $qm) < $NEARDIST) ? 1 : 0; } sub _isnear2 { my $gm= int(($_[0]+$_[1])/2); my $qm= int(($_[2]+$_[3])/2); my $dist= $_[4] || 2 * $NEARDIST; return (abs($gm - $qm) < $dist) ? 1 : 0; } # sort gene exon set by ref>start>stop # for exon array = ($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$xalt, @altexons) sub _exloc1sort { my $ar= $a->[0]; my $br= $b->[0]; $ar =~ s/[\_]+$//; $br =~ s/[\_]+$//; return # $ar cmp $br or # is this bad still? yes $a->[4] <=> $b->[4] or $a->[5] <=> $b->[5]; } # sort gene exon set by ref>start>stop # for genegroup array = ($ix, $geneid, $exons) # and genelocarray = [$ref,$tb,$te,$strand, $eq, $alt] my $genelocations={}; # temp global sub _geneloc1sort { my $ag= $genelocations->{$a->[1]}; # $a->geneid my $bg= $genelocations->{$b->[1]}; if ($ag && $bg) { return # $ag->[0] cmp $bg->[0] || ## this is bad ?? $ag->[1] <=> $bg->[1] || $ag->[2] <=> $bg->[2]; } elsif($ag) { return -1; } elsif($bg) { return 1; } else { return 0; } } sub clean_exonlist { my($alist, $dropnum, $keepmarks, $keeploc, $joiner, $wanthash, $wantfull)= @_; my %gfull=(); my @galt= grep( /\w/, split(/\s*[,]\s*/, $alist)); # may be only 1 my $inorder=1; my %gids= map{ my ($val,$full,$exeq,$exloc)=( $inorder++, $_, -1, ""); if($keepmarks == 2) { # sort by; is this working?? my $a= tr/@//d; my $p= tr/+//d; my $s= tr/*//d; my $o= tr/%#&//d; $val= 1 + $a * 2 + $p * 4 + $s * 1; #?? want alt matches 1st? #s/^[\+\*\@\%\#\&]*//; } elsif (!$keepmarks) { s/^[\+\*\@\%\#\&]*//; } # s/[\+\*\@]//g unless $keepmarks; s/\=([\d\-]+)$//; # and $exeq=$1 exeq added; testing s/\:([\d\-]+)$// unless $keeploc; # and $exloc=$1 s/\.[\d]+$// if $dropnum; $gfull{$_}= $full if($wantfull); $_,$val; } @galt; if($keepmarks == 2) { # sort by marks (++@ > ++* > other) ? @galt= sort {$gids{$b} <=> $gids{$a}} keys %gids; } else { # @galt= sort keys %gids; #? or try to keep input list order, esp. primary query id? @galt= sort {$gids{$a} <=> $gids{$b}} keys %gids; # input order **?? IS THIS BEST ??** } if($wantfull) { foreach (@galt) { $_= $gfull{$_} || $_; } } $joiner ||= ","; return ($wanthash) ? (\@galt, \%gids) : (wantarray) ? @galt : join($joiner, @galt); } sub print_genegroup1_gff { my( $geneid, $geneloc, $exons, $commonid, $geneidnum )=@_; my $nprint= 0; return $nprint unless ref $geneloc->{$geneid}; my($gref,$gb,$ge,$gstrand, $geq, $galt, $geneidnum2, $gnexons, $gnaltexon, $genescores, $geneinfo) = @{ $geneloc->{$geneid} }; my $gscore= "."; $geneidnum ||= $geneidnum2; $gref =~ s/[\_]+$//; # drop done flag return $nprint if($MIN_EXON_PRINT >0 && $gnexons && $MIN_EXON_PRINT >$gnexons); my $maintype= $genetype; # my $subtype= "match_part"; # or HSP for ggb my $matchID= $IDPREFIX.$geneid; # FIXME: IDPREFIX option my $gat= "ID=$matchID"; my $cid= $commonids{$commonid}; # commonid num my $ctype=""; (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g; if($commonid eq $geneid) { #? leave as is; use maintype diff # $ctype= $geqtype; $geqtype= "common"; $maintype= $genegrouptype; print "\n"; } my $source = $sourcetag; if(SOURCE_HAS_EQTYPE) { $source .= ".".$geqtype; # drop this? } else { if($EQlegend{$geq} =~ /far/) { $source= $sourcefar; } # elsif($EQlegend{$geq} =~ /near/) { $source= $sourcenear; } else { $source= $sourcenear; } # keep equal, altequal, .. w/ near group } # $gat .= ";ctype=$ctype" if($ctype); $gat .= ";tclass=$geqtype"; # replaces ctype= $gat .= ";cid=$cid"; $gat .= ";gid=$geneidnum" if($geneidnum); $gat .= ";nexons=$gnexons" if($gnexons); $gat .= ";altexons=$gnaltexon" if($gnaltexon); #? $gat .= ";info=$geneinfo" if ($geneinfo); if($galt && $commonid eq $geneid) { $gat .= ";ids=$galt"; } elsif($genescores) { $gat .= ";bestids=$genescores"; if($genescores =~ m/\[\s*([e\d.-]+)/) { $gscore=$1; } } elsif($galt) { $galt= clean_exonlist( $galt, 1, 0, 0, "."); # want to add gene_quality score here $gat .= ";altid=$galt" if($galt); } # move geqtype to $gat print join("\t",$gref, $source, $maintype, $gb, $ge, $gscore, $gstrand,".",$gat),"\n"; $nprint++; return $nprint unless (ref $exons); my($l_ref, $l_b, $l_e)=(0)x3; foreach my $exon (sort _exloc1sort @$exons) { # my($exnum, $xid, $ref, $tb, $te, $p, $eq, $xalt)= @$exon; my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$xalt, @altexons)= @{ $exon }; my $ismain= ($ref =~ m/[^\_][\_]{3}$/); # funky flag .. my $isalt = ($ref =~ m/[^\_][\_]{2}$/); # funky flag .. $ref =~ s/[\_]+$//; # drop done flag # next if($ref eq $l_ref && $tb eq $l_b && $te eq $l_e); #? need/want this? ; add $qid? # ($l_ref, $l_b, $l_e)= ($ref,$tb,$te); ## replace altpar=$galt with gid=$geneidnum like cid=$commonidnum # my $flag= $isalt ? "altexon=1" : ""; my $xat= "Parent=$matchID;cid=$cid"; $xat.= ";gid=$geneidnum" if($geneidnum); # $xat.= ";altpar=$galt" if($galt); $xat.= ";altexon=1" if($isalt); $xat.= ";mainx=1" if($ismain); # debug only?? ## HACK; FIXME: #see below# $exon->[3] .= ";gx=$gat"; if($qe =~ s/;gx=(\S+)//) { $xat.= ";gxclass=".$1; } print_exon1( $exon, $xat); ## , $gref, $gstrand); $nprint++; } return $nprint; } sub print_exon1 { my($exon,$xattrib)= @_; # ,$gref,$gstrand my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$xalt, @altexons)= @{ $exon }; $ref =~ s/[\_]+$//; # drop done flag (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g; my $xat= $xattrib; #revert to this# unless($xat =~ s/;cid=/;tclass=$xeqtype;cid=/) { $xat .= ";tclass=$xeqtype"; } if($qe =~ s/(;gx=\S+)//) {} # HACK; see above my $source = $sourcetag; if(SOURCE_HAS_EQTYPE) { $source .= ".".$xeqtype; # drop this? } else { if($xeqtype =~ /far/) { $source= $sourcefar; } # elsif($xeqtype =~ /near/) { $source= $sourcenear; } else { $source= $sourcenear; } # keep equal, altequal, .. w/ near group } #.. this is same as exon_allids() .. replace my $alt2= $xalt; foreach my $ax (@altexons) { $alt2 .= ",".$ax->[1]; # qid $alt2 .= ",".$ax->[9]; # alt } if ($FULL_EXONLIST) { # clean_exonlist($alist, $dropnum, $keepmarks, $keeploc, $joiner, $wanthash, $wantfull) $alt2= clean_exonlist( $alt2, 0, 0, 0, ",", 0, 1); $alt2 =~ s/[\+\*\@\%\#\&]*//g; # but dropmarks } else { $alt2= clean_exonlist( $alt2, 0, 0, 0); $alt2 =~ s/$qid[,]?//; } $xat .= ";xid=$qid"; $xat .= ",$alt2" if($alt2); # want this? # $xat .= ";xalt=$alt2" if($alt2); # want this? debug only? # $xat.= ";ix=$exnum"; ## dont need this with xid.exnum # move xeqtype to $xat print join("\t",$ref, $source, $exontype, $tb,$te,$p,$strand,".",$xat),"\n"; } ## this doesnt help find missing alt exon ids .. NCBI_GNO_462044 being one # sub altexon_recursion { # my($alt2, $altexons,$didexons)= @_; # # foreach my $ax (@$altexons) { # my($qid, $alt)= ($ax->[1],$ax->[9]); # next if($didexons->{$qid}++); # $$alt2 .= ",".$qid; # qid # $$alt2 .= ",".$alt; # alt # my $nx= scalar(@$ax)-1; # if($nx>9) { # my @more= $ax->[10..$nx]; # altexon_recursion($alt2, \@more, $didexons); # } # } # return $alt2; # ref, dont need # } ## ** gene-merge, like commons_merge, join gene/exon sets together where ## they overlap in region & gene ids ? ## require same-exon overlap? or just same-gene-region/id sub exon_allids { my( $exon, $exonorgeneid, $wanthash)= @_; # my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; my $asgeneid = ($exonorgeneid > 0) ? 1 : 0; my $fullexonid = ($exonorgeneid < 0) ? 1 : 0; my $realfullid = ($exonorgeneid < -1) ? 1 : 0; my $nx= @{$exon}; my $idlist = join ",", @{$exon}[1,9]; # could recursively make idlist here. see above altexon_recursion foreach my $ax ( @{$exon}[10 .. $nx-1] ) { $idlist .= ",".join ",", @{$ax}[1,9]; } return clean_exonlist( $idlist, $asgeneid, 0, $fullexonid,",",$wanthash, $realfullid); } sub merge_exons { my($iexons, $jexons, $domerge, $useidmatch, $jequalsi)= @_; my $debugx= $debug ? 20 : 0; my $ni= @$iexons; my $nj= @$jexons; $jequalsi= ($iexons == $jexons) unless($jequalsi); # my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; my ($nsame, $nolap, $nidsame)=(0,0,0); NEXTI: for(my $i=0; $i<$ni; $i++) { my ($igenea, $igeneh)= (); # wait for it# exon_allids($iexons->[$i], 1, 1); # warn "#D igenes[$i]=",join(",", sort keys %$igeneh),"\n" if($debugx-- > 0); #my $sameorlap = 0; my $j0= $jequalsi ? $i+1 : 0; for (my $j=$j0; $j<$nj; $j++) { my $sameorlap = _sameorlap1( @{$iexons->[$i]}[0,4,5,6], @{$jexons->[$j]}[0,4,5,6]); my $sameid= 0; if($sameorlap>0) { $nsame++; next NEXTI; } elsif($sameorlap<0) { $nolap++; next NEXTI; } elsif($useidmatch) { $igeneh or ($igenea, $igeneh)= exon_allids($iexons->[$i], 1, 1); my ($jgenes)= exon_allids( $jexons->[$j], 1, 1); # warn "#D jgenes[$j]=",join(",",@$jgenes),"\n" if($debugx-- > 0); foreach my $jgene (@$jgenes) { if($igeneh->{$jgene}) { $sameid=1; $nidsame++; next NEXTI; } } } } } #x ## requiring idsame + locsame doesnt seem to improve grouping; #x ## with either idsame or locsame, merge essentially all exons inside genegroup into one match model now #x ## need other way to separate out distinct gene models: by counting non-overlap dupls. of most common gene ids? my $canmerge= (($nsame + $nolap) >1 || $nidsame >0); #? # my $canmerge= ( ($nsame + $nolap) >1 || ($nidsame >0 && ($nsame + $nolap) >0) ); #? #warn "#D mx?$canmerge s=$nsame, o=$nolap, d=$nidsame\n" if $debugx-- > 0; return $canmerge unless($domerge && $canmerge); my @mexon=(); my @cjexons=(); push(@cjexons, @{$jexons}); # copy array ($nsame, $nolap, $nidsame)=(0,0,0); $debugx= $debug ? 1 : 0; my %ivalid= map{ $_,1} (0..$ni-1); for(my $i=0; $i<$ni; $i++) { next if($jequalsi && !$ivalid{$i}); my $sameorlap = 0; my $sameid= 0; my ($igenea, $igeneh)= (); # exon_allids($iexons->[$i], 1, 1); $nj= scalar(@cjexons); my $j0= $jequalsi ? $i+1 : 0; my $j= $j0; for (; $j<$nj; $j++) { next if($jequalsi && !$ivalid{$j}); $sameorlap = _sameorlap1( @{$iexons->[$i]}[0,4,5,6], @{$cjexons[$j]}[0,4,5,6]); if($sameorlap != 0) { last; } #x if($useidmatch) { $igeneh or ($igenea, $igeneh)= exon_allids($iexons->[$i], 1, 1); my ($jgenes)= exon_allids( $jexons->[$j], 1, 1); foreach my $jgene (@$jgenes) { if($igeneh->{$jgene}) { $sameid++; last; } } last if $sameid; } } if($sameorlap>0) { # merge jexon[j] into iexon[i] my $jexon= $jequalsi ? $jexons->[$j] : splice( @cjexons,$j,1); push( @{$iexons->[$i]}, $jexon); # add to altexon list push( @mexon, $iexons->[$i]); $ivalid{$j}= 0; # push(@mexon, $jexon ); # need to update each exon geneid, ... $nsame++; } elsif ($sameorlap<0) { # keep both? or pick longest/best eval/ my $jexon= $jequalsi ? $jexons->[$j] : splice( @cjexons,$j,1); my $ieval= $iexons->[$i]->[8]; my $jeval= $jexon->[8]; if($jeval < $ieval) { push( @$jexon, $iexons->[$i]); # add to altexon list push( @mexon, $jexon ); } else { push( @{$iexons->[$i]}, $jexon); # add to altexon list push( @mexon, $iexons->[$i]) } $ivalid{$j}= 0; $nolap++; } elsif ($sameid>0) { # but keep both exons push(@mexon, $iexons->[$i]); unless($jequalsi) { my $jexon= $jequalsi ? $jexons->[$j] : splice( @cjexons,$j,1); push(@mexon, $jexon ); $ivalid{$j}= 0; } $nidsame++; } else { # save iexon[i] , go to next i; what jexon to move? push(@mexon, $iexons->[$i]); } } #warn "#D mx: s=$nsame, o=$nolap, d=$nidsame\n" if $debugx-- > 0; # copy all remaining cjexons push(@mexon, @cjexons) unless($jequalsi); return \@mexon; # didmerge } =item merge_genegroup algo1 -- merges exons in gene-group based on (1) shared exon overlap (nex-overlap > 1 or >2) and (2) shared gene-ids -- essentially merges all exons in a genegroup now with this criteria (test more?) -- probably should replace this with algo2 =item merge_genegroup algo2 ** Before this, use algo1 for same-loc exon merging? -- start with all exons in gene group -- find common gene ids, and count non-overlapping dupl. exons in region -- max # genes in region == max nov-dupl. exons? -- need to allow for exons from alternate gene models as part of common gene set -- i.e. any predictor can err, rely on dupl. exon groups to build models =cut sub merge_genegroup { # algo2 my($ggroup, $commonid)= @_; my $debugx = $debug ? 20 : 0; my ($nmerged, $nsplit)= (0,0); my ($irow1, $geneid1, $geneidnum1, %ingeneids); my $mgroup = $ggroup; # #** Before this, use algo1 for same-loc exon merging? see instead merge_exons # ($mgroup, $nmerged)= merge_genegroup1($mgroup, $commonid, 1); # if($nmerged>0) { $nmerged=0; } # reset; report tho ? my $ng= @$mgroup; my @mexon= (); #*** FIXME: use existing gene-groups/gene-models and gene_quality scores to # guide gene_classifier, gene_split. either don't merge all exons, or keep # $ggroup around for comparing (run thru gene_quality here?) #1. merge all exons in genegroup foreach my $ig (0..$ng-1) { my($irow, $geneidi, $exonsi, $geneidnumi)= @ {$mgroup->[$ig]}; push(@mexon, @$exonsi); $ingeneids{$geneidi}++; ($irow1, $geneid1, $geneidnum1)= ($irow, $geneidi, $geneidnumi) unless($geneid1); $nmerged++; # == all } # ** try merge exons only not merge_genegroup1 .. my $mergedexons= merge_exons( \@mexon, \@mexon, 1, 0, 1); if ($mergedexons && ref($mergedexons) =~ /ARRAY/) { my $xinx= @mexon; my $xout= @$mergedexons; print $NOTEOUT "# merge_exons in=$xinx out=$xout \n" if $debug; @mexon = @$mergedexons; } @mexon= sort _exloc1sort @mexon; # sort by location my $nx= @mexon; print $NOTEOUT "\n# merge_genegroup[v2] $commonid in=$ng nx=$nx \n" if $debug; ($nx > 1) or return (wantarray) ? ($mgroup, $nmerged) : $mgroup; my( $n_dupl_nov, $n_diff_nov, $n_same_ov, $n_diff_ov)=(0) x 10; my (%dupgenes, %dupsamegene, %allgenes, %xlist, %gene2ix); for(my $ix=0; $ix<$nx; $ix++) { my ($iexona, $iexonh)= exon_allids($mexon[$ix], 0, 1); # want notFULL/loc, but exon ids here not gene? foreach my $ixid (@$iexona) { my($geneid,$exnum)= split_exonid($ixid,1); $exnum ||= 1; $allgenes{$geneid}{$exnum} ++; $allgenes{$geneid}{0} ++; $xlist{indx}{$geneid}{$ix}= $ix unless exists $xlist{indx}{$geneid}{$ix}; # same ix num $xlist{exnum}{$geneid}{$ix}= $exnum; } $xlist{indx}{all}{$ix}= $ix unless exists $xlist{indx}{all}{$ix}; # same ix num for (my $jx=$ix+1; $jx<$nx; $jx++) { my $sameorlap = _sameorlap1( @{$mexon[$ix]}[0,4,5,6], @{$mexon[$jx]}[0,4,5,6]); my $sameid = 0; #? do we also want gene-id matching? ##? want $altids{$geneid}{$exnum} : see process_genegroup if ($sameorlap == 0) { my ($jexids)= exon_allids( $mexon[$jx], 0, 1); # notFull exonid my $allj= 0; foreach my $jexid (@$jexids) { if( $iexonh->{$jexid} ) { $sameid++; if ($sameorlap == 0) { # ?? need to use j->qid instead of j-altids here? # need to use this instead of $geneid from any alternate exon id? # or force set mexon.1 to be matching $jexid/iexid ? my $jqid= $mexon[$jx]->[1]; # qid; my $iqid= $mexon[$ix]->[1]; # qid my($geneid,$exnum)= split_exonid($jexid,1); $exnum ||= 1; $dupgenes{$geneid}{$exnum} ++; $dupgenes{$geneid}{0} ++; for my $gn ('all',$geneid) { if($gn ne 'all' || $allj++ == 0) { $xlist{indx}{$gn}{$jx}= $ix; # same ix num ## $xlist{dup}{$gn}{$ix} ++; exists $xlist{alt}{$gn}{$ix} or $xlist{alt}{$gn}{$ix}=[]; exists $xlist{alt}{$gn}{$jx} or $xlist{alt}{$gn}{$jx}=[]; push( @{$xlist{alt}{$gn}{$ix}}, $jx); push( @{$xlist{alt}{$gn}{$jx}}, $ix); } } ## ** save this code ** may need again # if($iqid ne $jexid && $iqid =~ m/^$geneid/) { # # other exon of same gene at other location is duplicate: gene model combining 2+ genes ? # my($igeneid,$ixnum)= split_exonid($iqid,1); # unless( !$ixnum || $exnum == $ixnum # || ($dupsamegene{$geneid}{$ixnum} && $dupsamegene{$geneid}{$ixnum} =~ m/\b$exnum,/) # ) { # $dupsamegene{$geneid}{$ixnum} .= $exnum."," ; # $dupsamegene{$geneid}{$exnum} .= $ixnum."," ; # } # } } } } } if ($sameorlap == 0 && $sameid > 0) { # duplicate non-overlapping $n_dupl_nov ++; } elsif ($sameorlap == 0 && $sameid == 0) { # diff exons non-overlapping $n_diff_nov ++; } elsif ($sameorlap != 0 && $sameid > 0) { # same (altid)exon overlapping ** IGNORE, none $n_same_ov ++; # these are all 0; skip sameid test } elsif ($sameorlap != 0 && $sameid == 0) { # diff exons overlapping $n_diff_ov ++; } } } # move this out to own sub? input: dupgenes, allgenes, dupsamegene # $mgroup, @mexon, (%dupgenes, %dupsamegene, %allgenes); my @merged= (); my ($GC_in, $GC_gx, $GC_ix, $GC_jx, $GC_jbe) = gene_classifier( \@mexon, \%xlist, "all"); # 2. now split @mexon into $maxdup @merged gene models ? # max n_dupl_nov may not be best genemodel choice # mark gene models that have internal duplicated exons (and many) as likely combined tandems # that need splitting .. delete $ingeneids{$geneid1}; ## temp till get geneid renames fixed my $didsplit = gene_split (\@mexon, \@merged, $GC_gx, $GC_ix, $GC_jx, $GC_jbe, \%ingeneids); $nsplit += $didsplit; ## should == 1+max value @gx push( @merged, [$irow1, $geneid1, \@mexon, $geneidnum1 ] ) if(@mexon); # any remaining @mexon go here $mgroup= \@merged; print $NOTEOUT "# merge_genegroup $commonid nmerge=$nmerged nsplit=$nsplit\n" if $debug; return (wantarray) ? ($mgroup, $nmerged) : $mgroup; # replaces $ggroup } =item cut from merge_genegroup keep this code for a bit, was useful for classifying duplicate genes may want to reuse ##### see gene_classifier() ; drop # ............ drop this part as not useful ?? use some of it for gene_classifier ?? # if(0) { #........................................................... # my ($maxg,$maxe,$maxdup,$maxexlist)=(0) x 10; # my @dupgenes= sort keys %dupgenes; # my $ndg= @dupgenes; # my %dupscore=(); # foreach my $gn (@dupgenes) { # # check that $gn == common id? # my @dupex = sort{$a<=>$b} keys %{$dupgenes{$gn}}; shift(@dupex); # 0 # my @allex = sort{$a<=>$b} keys %{$allgenes{$gn}}; shift(@allex); # 0 # # my $exshow= join ", ", map{ # # "$_:". ($dupgenes{$gn}{$_}||0) ; # "$_:". ($dupgenes{$gn}{$_}||0) ."/". $allgenes{$gn}{$_} ; # } @allex; # my $comgn = $IDmap->{$gn}{common} || $gn; # # my $cgid= $commonids{$comgn} || $geneidnum1; # my $exonids= $IDmap->{$gn}{exons} || $IDmap->{$comgn}{exons} || []; #? # my $totexons= scalar( @$exonids ); # my $note= ($comgn ne $gn) ? " (comgn=$comgn)" : ""; # # my $score_dupsame= 0; # if($dupsamegene{$gn}) { # my $smax=1; # my @exonsame= map { # my $m= $dupsamegene{$gn}{$_}; # $m =~ s/,*$//; # my $nm= $m =~ tr/,/,/; $nm++; # $smax= $nm if($nm>$smax); # "$_=$m"; # } sort keys %{$dupsamegene{$gn}}; # $score_dupsame = scalar(@exonsame) * $smax; # $dupscore{$gn}{nsame}= $smax; # best n-genes here # # $note .= " (exdup: "; # $note .= join ", ", @exonsame; # $note .= ")"; # } # # my ($ndupe, $nexons)= (scalar(@dupex), $totexons); # my $score_allhere = scalar(@allex) / $totexons ; # far models may have only a few located here # my $score_dupcover = ($ndupe / $nexons); #? # # ## another good score: median number of dupl. exons, or max with 2,3+ exons? # ## * should look for exon clusters (ordered exnums, separated by intergene space) # # my $dmax= 0; # foreach my $ex (@dupex) { # my $nd= $dupgenes{$gn}{$ex}; # $dmax= $nd if ($nd > $dmax); # ($maxg,$maxe,$maxdup,$maxexlist)= ($gn,$ex,$nd,join(",", @dupex)) # if( $nd > $maxdup ); # } # my $score_dupmax= $dmax; # # my $score= $dupscore{$gn}{score}= 1 # + $score_dupmax / $nexons # 1..n # + $score_dupcover * 5 # 0..1 * nexons? // Weight high # + $score_allhere # 0..1 * $nexons? # - $score_dupsame / $nexons # n dups in gene! (0,1,2) * nexons // reduce by this? # ; # my $sv= join ",", map{ sprintf "%.2g",$_; } # ($score_dupmax/ $nexons, $score_dupcover, $score_allhere, -$score_dupsame/ $nexons); # $score= sprintf "%.2g",$score; # print $NOTEOUT "# "," $gn sc=$score ($sv); ex=[ $exshow] $ndupe of $nexons $note\n" if $debug; # # if($debug) { # look at each exon in this group??? how to cluster into gene models # # print $NOTEOUT "# Exons for $gn .............. \n"; # for(my $ix=0; $ix<$nx; $ix++) { # # next unless exists $xlist{indx}{$gn}{$ix}; # my $jx= $xlist{indx}{$gn}{$ix}; # my @ixalt= ( exists $xlist{alt}{$gn}{$ix} ) ? @{$xlist{alt}{$gn}{$ix}} : (); # if( exists $xlist{alt}{$gn}{$jx} ) { push(@ixalt, @{$xlist{alt}{$gn}{$jx}}); } # my %ixalt= map{$_=>1} @ixalt; my $ixalt= "ija:". join ",", sort keys %ixalt; # # my $exon= $mexon[$ix]; # my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; # my $alte= ""; # # $alte= clean_exonlist( $alt, 0, 0, 0); # foreach my $dgn (@dupgenes) { # next if($qid =~ /^$dgn/); # if($alt =~ /($dgn[\d\.]+)/) { $alte .= $1.","; } # } # # my $MET= mark_startstop($strand,$alt); # print $NOTEOUT "# ",join("\t", "i".$ix, "j".$jx.$MET, $ixalt, $qid,$tb,$te,$strand,$eq, $alte), "\n"; ## ,$alt # } # } # # } # # print $NOTEOUT " # # ndupl-genes: $ndg ; max dupl. exon: $maxg.$maxe n=$maxdup; exons=$maxexlist; # # overlaps: n_dupl_nov=$n_dupl_nov; n_diff_nov=$n_diff_nov; n_same_ov=$n_same_ov; n_diff_ov=$n_diff_ov \n" # if $debug; # } #............................................................. =cut sub mark_startstop { my ($strand,$jalt)=@_; my @jalt= clean_exonlist( $jalt, 0, 0, 0); my ($nstart,$nstop)=(0,0); my $rstrand= ($strand eq '+') ? -1 : ($strand eq '-') ? +1 : -$strand; foreach my $xid (@jalt) { my $xo= $IDmap->{exonstrand}->{$xid} || 0; # or next; my $rev= ($xo && $xo == $rstrand); if($IDmap->{exonstart}{$xid}) { $rev ? $nstop++ : $nstart++; } if($IDmap->{exonstop}{$xid} ) { $rev ? $nstart++ : $nstop++; } } return wantarray ? ($nstart,$nstop) : (('^') x $nstart) . (('*') x $nstop); } =item GENE MODEL CLASSIFIER ($gx,$ix,$jx,$jbe) = gene_classify(\@mexon, \%xlist, ); mark gene model transition ($gx num) by my ($in, @gx, @jx, @or, @jbe, @ix) -- ix = multi-exon index, sorted by location; -- jx = index to base exon (where duplicates) -- transition at jx <= jx-1 with run of jx+1..n < jx-1..-n -- must distinguish for, reverse transitions: -- use start/stop exons (jbe); problematic data just now -- add in qid,altid exon numbers : changes in runs suggest new gene model input: @mexon [ix=0..nx-1] : all exons in gene group (overlapping redundants removed) $xlist{indx}{$GROUP}{$ix} : (where $GROUP = "all"; also have per dup-gene model id) $xlist{alt}{$GROUP}{$ix} : alternate ix indices (e.g. dupl. exons) =cut sub gene_classifier { my( $mexon, $xlist, $GROUP )= @_; my( $in, @gx, @ix, @jx, @jbe, @or, @gexn, @dx); # return all but @or ? $GROUP ||= "all"; # can we have this as param: want per geneid ? my $gx= 0; my $nx= scalar(@$mexon); $in= 0; ## should we guess for #genes/predictor-type? my @geneids = sort keys %{$xlist->{exnum}}; my $ningenes= scalar(@geneids); # not to be confused with #genes we find here my ($ltb,$lte)=(0,0); my %exnums; for(my $ix=0; $ix<$nx; $ix++) { next unless exists $xlist->{indx}{$GROUP}{$ix}; my $jx= $xlist->{indx}{$GROUP}{$ix}; #?? use alt/ija: also, some jx dont have lowest base exon value my @ixalt= ( exists $xlist->{alt}{$GROUP}{$ix} ) ? @{$xlist->{alt}{$GROUP}{$ix}} : (); if( exists $xlist->{alt}{$GROUP}{$jx} ) { push(@ixalt, @{$xlist->{alt}{$GROUP}{$jx}}); } # ix,mexon should be location sorted here my $exon= $mexon->[$ix]; my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; my $minj= $jx; foreach my $kx (@ixalt) { $minj= $kx if ($kx < $minj); } $jx[$in]= $minj; $ix[$in]= $ix; $gx[$in]= 0; ## set all exons to same gene id 0 to start; $or[$in]= $strand; #? convert to number $exnums{$minj}++ if ($minj<5); $dx[$in]= ($in==0) ? 0 : ($tb - $lte); my ($jb,$je)= mark_startstop($strand,$alt); $jbe[$in][0]= $jb; $jbe[$in][1]= $je; # add table of gexn[ix][igene]=exnum to find runs / transitions from avail gene models # this should be comparable to gx[ix] as alternate gene models to test against # want to use existing prediction models where they are good (maybe adding exons?) for (my $ig=0; $ig < $ningenes; $ig++) { my $gn= $geneids[$ig]; my $exnum= $xlist->{exnum}{$gn}{$ix} || 0; $gexn[$in][$ig]= $exnum; } $in++; ($ltb,$lte)= ($tb,$te); } return ( $in, \@gx, \@ix, \@jx, \@jbe) if ($in<2); # or more? my ($mingene1,$mingene2)= sort{$b<=>$a} values %exnums; my $mindupgenes= $mingene2 || 0; my $ntransit= 0; my $mintr= ($in < 6) ? 2 : 3; #?? prefer many groups on weak info? my $lastjmax= 0; # num exons in base gene my $basejmax=0; my $basec= 0; #? use base gene model exon count? my ($avintron, $nintron, $maxintron, $intergenei, $intergened)= (0) x 10; if($in > 1) { my @dxinter= sort{$b <=> $a} @dx; $maxintron= $dxinter[0]; my $mode = int($nx/3); if($mode>0) { for (my $m= $mode; $m < 2*$mode; $m++) { $avintron += $dxinter[$m]; $nintron++; } $avintron= int( $avintron/ $nintron); for (my $m= 1; $m < $mode; $m++) { # biggest interexon dists; look for big step my $dxm = $dxinter[$m] or last; next if ($dxm > $avintron * 9); if( ( $dxinter[$m-1] / $dxm ) > ( $dxm / $dxinter[$m+1] ) ) { $intergened= $dxinter[$m-1]; $intergenei= $m; last; } # my $ratio= $dxinter[$m-1] / $dxm ; # # my $diff = $dxinter[$m-1] - $dxm; # #? check maxintron .. dmx .. aveintron # if($ratio > 3) { $intergened= $dxinter[$m-1]; $intergenei= $m; last; } # good guess? } } if($intergened) { for (my $i=1; $i<$in; $i++) { my $genei= $gx[$i]; my $disti= $dx[$i]; my $rev = ($or[$i] ne $or[$i-1]); # + > - or - > + transition if($disti >= $intergened) { my $gx= 1 + abs($genei); if( $rev ) { $gx= -$gx; } # ?? will this work for( my $k=$i; $k<$in; $k++) { $gx[$k]= $gx; } $ntransit++; } } } } my $note=""; $note= "intergene($ntransit);" if $ntransit>0; my $done= 0; $done=1 if ($mindupgenes>0 && $ntransit >= $mindupgenes); # or look for more? ## ** SHOULD add intron/intergene distance as important model transition criteria ## ** better algo may be to est. # genes in group, and start split at biggest inter-exon spaces unless($done) { for (my $i=1; $i<$in; $i++) { my $tr = 0; # transition quality my $rev = ($or[$i] ne $or[$i-1]); # + > - or - > + transition my $genei= $gx[$i]; $note .= "."; next if($genei != $gx[$i-1]); # already did transit; any calcs to make ? next if($i<$in-1 && $genei != $gx[$i+1]); # coming to transit; any calcs to make ? my $priori= ($genei == 0) ? 0 : abs($genei)-1; my $exonc = 0; for (my $m=$i; $m>=0; $m--) { $exonc++; last if($gx[$m]!=$genei); } my $priorc= 0; for (my $m=$i; $m>=0; $m--) { next if($gx[$m]==$genei); $priorc++; last if(abs($gx[$m])!=$priori); } my $genedir=0; if (($rev && $genei < 0) || (!$rev && $genei > 0)) { $genedir= 1; } elsif (($rev && $genei > 0) || (!$rev && $genei < 0)) { $genedir= -1; } if ( $rev ) { $tr++; $note .= "rev;"; } elsif ( $genei == 0 ) { $tr++ if( $jx[$i] <= $jx[$i-1] ); } # 1st gene elsif ( $genei < 0 ) { $tr++ if( $jx[$i] <= $lastjmax && $jx[$i] > $jx[$i-1] ); } # rev gene 2+ elsif ( $genei > 0 ) { $tr++ if( $jx[$i] <= $lastjmax && $jx[$i] <= $jx[$i-1] ); } # for gene 2+ #?? how to weight and when? if($avintron > 0 && $dx[$i] > $avintron * 9) { $tr += 9; $note .= "intergene;"; } elsif($avintron > 0 && $dx[$i] > $avintron * 4) { $tr += 1; $note .= "inter4;"; } #? higher wt if ( $tr != 0 ) { # 1. find max jx for prior gene models; all? # ensure $genei>0 doesnt have spurious transits w/in a gene model run # count exons in base model; maybe want average/range basec for all so far? $basec=0; for (my $m=0; $m<$i; $m++) { $basec++; last if($gx[$m]!=0); } $basejmax= 0; for (my $m=0; $m<$i; $m++) { $basejmax= $jx[$m] if($gx[$m]==0 && $jx[$m]>$basejmax); } for (my $m=0; $m<$i-5; $m++) { $lastjmax= $jx[$m] if($jx[$m]>$lastjmax); } ## dont allow bouncing around from old to new exons w/in gene model ## cause transition to new $gx: e.g. for basec ~10, # i18:j18^; i19:j0 ; i20:j1 ; i21:j21< ; i22:j2; i23:j23; .. # # add in use of runs and transitions in geneid .exnum < # foreach my $k (0..5) { # my ($f,$r) = ($i + $k, $i - 1 - $k); # last if ($f >= $in || $r < 0); # # @runs[$f]= @{gexn[$f]}; #? # # for (my $ig=0; $ig < $ningenes; $ig++) { $runs[$f][$ig]= $gexn[$f][$ig]; } # } #? use start,stop mark to boost $tr #? distinguish b,e for i,i-1? if ($i < 2) { } # skip jbe #? elsif( $jbe[$i][0]>0 || $jbe[$i][1]>0 ) { $tr++; $note .= "be;"; } elsif( $jbe[$i-1][0]>0 || $jbe[$i-1][1]>0 ) { $tr++; $note .= "be;";} ## give tr start at j0,j1,j2 high tr weight? #? only do this for $i, not $i+k? if ($genei != 0 && $i > $basejmax && $jx[ $i] < 4) { $tr += 2; $note .= "j0;"; } KLOOP: foreach my $k (1..5) { #? include 0,$i for tests? my ($f,$r) = ($i + $k, $i - 1 - $k); last if ($f >= $in || $r < 0); if ($genei != 0 && $jx[ $f] == $ix[ $f]) { # new exon; no info? cant compare to jx[ r] if( $k>2 && $jbe[$k][0] + $jbe[$k][1]>0 ) { $tr++; $note .= "be;";} # ?? next; } #? check also that 1..5 are run of ascending/descending jx values if($genedir < 0) { # or $rev? reversed if($jx[ $r] <= $basejmax && $jx[ $f] > $jx[ $r]) { $tr=0; $note .= "r0;"; last KLOOP; } # spurious*?? if($jx[ $f] < $lastjmax && $jx[ $f] > $jx[ $r]) { $tr++; $note .= "r 0) { # want to see base exons repeated # .. need to test jx.r against all jx[i..k] for this if($jx[ $r] <= $basejmax) { for (my $f1= $i; $f1<=$f; $f1++) { if( $jx[ $f1] < $basejmax && $jx[ $f1] > $jx[ $r]) { $tr=0; $note .= "r0;"; last KLOOP; } # spurious*?? } } if($jx[ $f] < $jx[ $i]) { $tr=0; last KLOOP; } # better transition to follow? elsif($jx[ $f] < $lastjmax && $jx[ $f] <= $jx[ $r]) { $tr++; $note .= "f $jx[ $i] && $jx[ $f] <= 2+$k) { $tr++; $note .= "j1;"; } ## a base exon? } else { if($jx[ $f] <= $jx[ $r]) { if (!$rev && $jx[ $f] < $jx[ $i]) { $tr--; $note .= "fi-;";} #? elsif ($rev && $jx[ $f] > $jx[ $i]) { $tr--; $note .= "fi-;";} else { $tr++; $note .= "fr1;"; } } } if($rev) { if($or[ $f] ne $or[ $r]) { $tr++; $note .= "frev;"; } #? only if $rev } else { if($or[ $f] ne $or[ $r]) { $tr--; $note .= "rrev;"; } #? reduce because $i is not best tr? } } # transit if iexon > basec and jbe start/stop ? if($genei != 0 && $tr < $mintr && $basec>0 && $exonc >= $basec) { ## || $exonc >= $priorc ? if( $jbe[$i][0]>0 || $jbe[$i][1]>0 ) { $tr++; $note .= "be2;"; } # += $mintr; # count both if( $jbe[$i-1][0]>0 || $jbe[$i-1][1]>0 ) { $tr++; $note .= "be2;"; } # += $mintr; } # check we don't transit after very partial gene model .. # if ($priorc > 0 && $basec > 0 && $priorc < $basec/4) { $tr= 0; } # is this useful? if($tr >= $mintr) { my $gx= 1 + abs($genei); for (my $m=0; $m<$i-5; $m++) { $lastjmax= $jx[$m] if($jx[$m]>$lastjmax); } $note .= "/g$gx,i$i "; if( $rev ) { $gx= -$gx; } # ?? will this work # fixme; preserve old gx -- add to gx? for( my $k=$i; $k<$in; $k++) { if( abs($gx[$k]) > abs($genei) ) { my $grev= ($gx[$k]<0) ? -1 : 1; $gx[$k]= $grev * ( 1 + abs( $gx[$k]) ); } else { $gx[$k]= $gx; } } $ntransit++; } } } } # loop; done if($debug>1) { my @dupgenes= sort keys %{$xlist->{alt}}; print $NOTEOUT "# Exons for $GROUP .............. [method=$note] \n"; for(my $i=0; $i<$in; $i++) { my $ix= $ix[$i]; my $jx= $jx[$i]; ## $xlist{indx}{$GROUP}{$ix}; my $gx= $gx[$i]; my @ixalt= ( exists $xlist->{alt}{$GROUP}{$ix} ) ? @{$xlist->{alt}{$GROUP}{$ix}} : (); if( exists $xlist->{alt}{$GROUP}{$jx} ) { push(@ixalt, @{$xlist->{alt}{$GROUP}{$jx}}); } my %ixalt= map{$_=>1} @ixalt; my $ixalt= "ija:". join ",", sort keys %ixalt; my $exon= $mexon->[$ix]; my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; my $alte= clean_exonlist( $alt, 0, 0, 0); # $qid should be in $alte now; drop as col my $MET= mark_startstop($strand,$alt); print $NOTEOUT "# ", join("\t", "i".$ix, "j".$jx.$MET, "g".$gx, $ixalt, $tb,$te,$strand,$eq,$alte), "\n"; } } return ( $in, \@gx, \@ix, \@jx, \@jbe); } =item sub gene_split $nsplit = gene_split (\@mexon, \@merged, \@gx, \@ix, \@jx, \@jbe, \%ingeneids) input: @mexon (all group exons) from gene_classifier: @gx (gene number); @ix (@mexon index), @jx, @jbe \%ingeneids == temp fix: list of valid geneids to re-use output: @merged genegroups = list of [ $irow1, $geneid, \@gexon, $geneidnum ] @mexon reduced by all @gexon group exons =cut sub gene_split { my( $mexon, $merged, $gx, $ix, $jx, $jbe, $ingeneids)= @_; my $nsplit= 0; return 0 if ($gx->[-1] == 0); # only 1 gene model in group; my $in= scalar(@$ix); # same as scalar(@$mexon) ? my %geneids=(); my %gexon=(); for(my $i=0; $i<$in; $i++) { my $ix= $ix->[$i]; my $jx= $jx->[$i]; ## $xlist{indx}{$GROUP}{$ix}; my $gx= $gx->[$i]; # should be ordered by $i $gx= -$gx if($gx<0); # drop - from geneindex; should be ok my $exon= $mexon->[$ix] or next; $mexon->[$ix]= undef; ## add attribute to exon w/ gene class info: # my $MET= (('^') x $jbe->[$i][0]) . (('*') x $jbe->[$i][1]); # for gff output; change 'g4j0^**i44' to 'g4j0bei44' ^^^^=b; ***=e my $mat= ($jbe->[$i][0] ? 'b' :'') . ($jbe->[$i][1] ? 'e' : ''); my $gat= join("", "g".$gx, "j".$jx.$mat, "i".$ix, ); $exon->[3] =~ s/;gx=\S+//; # drop old $exon->[3] .= ";gx=$gat"; # FIXME: append to $p/8 ; split out for print; alt is $qe/3 ? my $geneid= split_exonid($exon->[1],1); $geneids{$gx}{$geneid}++; ## FIXME: here or elsewhere, measure each alt geneid(exons) for consistency with found gene model ## count exons-matched / total-exons exists $gexon{$gx} or $gexon{$gx}=[]; push( @{$gexon{$gx}}, $exon); } foreach my $gx (sort keys %gexon) { my @gexon= @{ $gexon{$gx} }; #?? or next; my ($geneid)= sort{ $geneids{$gx}{$b} <=> $geneids{$gx}{$a} } keys %{$geneids{$gx}}; $geneid = "g$gx.".$geneid; # new id here .. $gx num? ## FIXME: temp-reuse input ingeneids or ****FIX common_merge # my ($inid)= sort keys %$ingeneids; # if ($inid) { $geneid= $inid; delete $ingeneids->{$inid}; } exists $IDmap->{$geneid}{num} or $IDmap->{$geneid}{num}= ++$GENEID_NUM; my $geneidnum= $IDmap->{$geneid}{num}; push( @$merged, [ 1, $geneid, \@gexon, $geneidnum ] ); # irow == 0 ok? not used? $nsplit++; } # now remove undefs from @mexon my $mx= scalar(@$mexon); for (my $ix= $mx-1; $ix >= 0; $ix--) { splice(@$mexon, $ix, 1) unless(defined $mexon->[$ix]); } return $nsplit; } =item gene_quality measure each alt geneid(exons) for consistency with found gene model count exons-matched / total-exons ; exons-lost, exons-gained, start/stop exon? add as geneloc field? =cut use constant QUALSUMID => "QUALSUMID"; # or $groupgeneid sub gene_quality { my($ix, $groupgeneid, $exons, $geneidnum)= @_; my %genequal=(); my %geneids=(); my $ixnum= 0; my $niexons= scalar(@$exons); my ($neqsame, $neqnear, $neqfar, $nrepeat)= (0) x 10; # count only if exon has such flag, not all ids w/ it foreach my $exon (@$exons) { $ixnum++; # my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; my ($eqsame, $eqnear, $eqfar, $isrepeat)= (0) x 10; # count only if exon has such flag, not all ids w/ it my ($iexona, $iexonh)= exon_allids( $exon, -2, 1); # full exonids foreach my $xid (@$iexona) { my ($gid,$exnum,$exloc,$exeq)= split_exonid($xid,1); ## keep exloc also to match overlap? or use exeq $geneids{$gid}{exons}{$exnum}++; $geneids{$gid}{exeq}{$exeq}++; # only per gene? $eqsame= 1 if($exeq == 1); $eqnear= 1 if($exeq == -2); $eqfar= 1 if($exeq == 0); # problem w/ undef? -1 is missing val $isrepeat=1 if($REPEATSKIP && $IDmap->{repeatlist}->{"$gid.$exnum"}); # is this exon score or gene score # $nrepeat= $IDmap->{repeatlist}->{$gid} # if($REPEATSKIP && $IDmap->{repeatlist}->{$gid}); # is this exon score or gene score # ^ count be one exon match from rep group among otherwise non-rep exons; # drop value here for isrepeat, i.e. nexons marked as repetitive # $geneids{$gid}{matchto}{$ixnum}++; # {$exnum}{$ixnum}++ ? # $geneids{$gid}{full} .= "$xid,"; } ## this is a per exon score; want count/exongroup $neqsame += $eqsame; $neqnear += $eqnear; $neqfar += $eqfar; $nrepeat += $isrepeat; } $genequal{QUALSUMID}{nexons}= $niexons; # useless $genequal{QUALSUMID}{eqsame}= $neqsame; $genequal{QUALSUMID}{eqnear}= $neqnear; $genequal{QUALSUMID}{eqfar }= $neqfar; $genequal{QUALSUMID}{repeats}= $nrepeat if($REPEATSKIP); $genequal{QUALSUMID}{fullgenes}= 0; # my $exonsizemap = $IDmap->{sizemap}; # exon{$xid1} size my @geneids= sort keys %geneids; my $igene=0; foreach my $gid (@geneids) { my($hasstart,$hasstop,$exlost,$exhas)=(0) x 10; $igene++; my $exonids = $IDmap->{$gid}{exons} || []; my $ngexons = @$exonids || 1; my $exoncover = $niexons / $ngexons; foreach my $xnum (@$exonids) { my $gxid= "$gid.$xnum"; $hasstart=1 if($IDmap->{exonstart}{$gxid}); #? count multiples? $hasstop= 1 if($IDmap->{exonstop}{$gxid} ); if( $geneids{$gid}{exons}{$xnum}) { $exhas++; } else { $exlost ++; } ## my $xo= $IDmap->{exonstrand}->{$xid} || 0; # or next; #? count exon-sizes vs match-lengh? } my $exgained= ($niexons < $exhas) ? 0 : $niexons - $exhas; my $hasends= $hasstart + $hasstop; my @eqtype= sort{ $geneids{$gid}{exeq}{$b} <=> $geneids{$gid}{exeq}{$a} } keys %{$geneids{$gid}{exeq}}; my $eqtype= $eqtype[0]; $genequal{$gid}{eqtype}= $eqtype; $genequal{$gid}{cover} = $exoncover; $genequal{$gid}{gained}= $exgained; $genequal{$gid}{lost} = $exlost; $genequal{$gid}{ends} = $hasends; =item FIXME: qual score x 1. weight loss, gain separately: good score minimize both DONE x 2. use related eq -2 near tandems to boost score , DONE BUT use this only for PREDICT_TYPE groups (? which we don't know here?), e.g. NCBI_GNO* vs Dappu1* vs DGIL_SNO e.g. eq=1 + 4 x eq=-2 >> eq=1 + 1 x eq=-1 ?? include entire gene group in qual score calcs ? ** should use gene_quality to influence gene_split =cut my $eqscore= ($eqtype == 1) ? 1 : ($eqtype == -2) ? 0.5 : 0; my $nolostscore= ($ngexons - $exlost) / $ngexons; # 0..1 my $nogainscore= ($niexons - $exgained) / $niexons; # 0..1 ; weight this less than lostscore for quality? if($exoncover > 1) { $exoncover= 1/$exoncover; } # for score ?? $genequal{QUALSUMID}{fullgenes}++ if ($nolostscore > 0.75); #? what score, 3/4 of gene exons sounds good # my $score0= $genequal{$gid}{score}= ($exoncover + $hasends/2 + $eqscore) / 3; # 0..1 + 0..1 > 0..1 my $score= $genequal{$gid}{score}= ( ($nolostscore * 0.8 + $nogainscore * 0.2) + $hasends/2 + $eqscore) / 3; } $genequal{QUALSUMID}{geneids}= $igene; use constant PREFER_TANDEMS => 1; # regularity in nature is a signal; these are useful if(PREFER_TANDEMS) { # look at all genes from prediction groups, score higher if several eq-2 from same group; my %gmethod=(); my $ntand=0; my $nmodel=0; foreach my $gid (@geneids) { (my $gmeth= $gid) =~ s/\d//g; # drop all digits; will this work? my $eq= $genequal{$gid}{eqtype}; #? only care about -2? # $gmethod{$gmeth}{eqtype}{$eq}++; $gmethod{$gmeth}{tandems}++ if($eq == -2); ## or +4 ? or other? $gmethod{$gmeth}{models}++; $ntand++ if($eq == -2); $nmodel++; } my ($tanmax, $genemax)= (0,0); foreach my $gid (@geneids) { (my $gmeth= $gid) =~ s/\d//g; my $models= $gmethod{$gmeth}{models} || 0; my $tandems= $gmethod{$gmeth}{tandems} || 0; $genequal{$gid}{tandems} = $tandems; $tanmax = $tandems if($tanmax<$tandems); $genemax= $models if($genemax<$models); next unless($ntand>0); my $tscore= $genequal{$gid}{tanscore}= $tandems / $ntand; # / $gmethod{$gmeth}{models}; #? or mtand / ntand ? #?? add in $mscore= $gmethod{$gmeth}{models} / $nmodel; ## standardize to 1? my $oscore= $genequal{$gid}{score}; my $nscore= $genequal{$gid}{score} = ($oscore * 0.7 + $tscore * 0.3 ); } $genequal{QUALSUMID}{methods} = scalar(keys %gmethod); $genequal{QUALSUMID}{maxgenes} = $genemax; $genequal{QUALSUMID}{tandems} = $tanmax; # est. max tandems to main gene in region for any method; +1 for main gene # so ($tanmax + 1) / $genemax is some est. of num tandem genes in these regions (picked for their tandemness) } my @bestgenes= sort{ $genequal{$b}{score} <=> $genequal{$a}{score} } @geneids; my $debugout= $debug>1; # add as param if($debugout) { print $NOTEOUT "# ",(".") x 50,"\n"; print $NOTEOUT "# gene_quality for $groupgeneid nexons=$niexons; ","\n" ; my $igene=0; foreach my $gid (@geneids) { $igene++; my $eqtype= $genequal{$gid}{eqtype}; my $exoncover= $genequal{$gid}{cover} ; my $exgained= $genequal{$gid}{gained}; my $exlost= $genequal{$gid}{lost} ; my $hasends= $genequal{$gid}{ends} ; my $tanscore= $genequal{$gid}{tanscore} || 0; $tanscore= sprintf("%.2f",$tanscore); my $score= $genequal{$gid}{score}; $score= sprintf("%.3f",$score); $exoncover= sprintf("%.3f",$exoncover); print $NOTEOUT "# gene_quality $gid($igene): " ."s=$score; ends=$hasends; loss/gain=$exlost/$exgained; eq=$eqtype; tand=$tanscore; cov=$exoncover\n"; } print $NOTEOUT "# ",(".") x 50,"\n" ; } return (\@bestgenes, \%genequal); } =item merge_genegroup1 : old # old algorithm; not as good as algo2 =cut # sub merge_genegroup1 { # algo1 # my($ggroup, $commonid, $noidmatch)= @_; # # my $debugx = $debug ? 20 : 0; # my $nmerged= 0; # my $mgroup = $ggroup; # my $imerged= 0; # my $iloop=0; # my $ng= @$mgroup; # # print $NOTEOUT "# merge_genegroup[v1] $commonid in=$ng skipid=$noidmatch\n" if $debug; # # LOOPIJ: do { $iloop++; # $ng= @$mgroup; # $imerged= 0; # my @merged= (); # my %unmerged= map{$_,1} (0..$ng-1); # # NEXTI: foreach my $ig (0..$ng-1) { # my($ix, $geneidi, $exonsi, $geneidnumi)= @ {$mgroup->[$ig]}; # #warn "#D mrg[$iloop]? $ig/$geneidi $unmerged{$ig} \n" if $debugx > 0; # next unless($unmerged{$ig}); # # NEXTJ: for my $jg ($ig+1 .. $ng-1) { # my($jx, $geneidj, $exonsj, $geneidnumj)= @ {$mgroup->[$jg]}; # next unless($unmerged{$jg}); # # # 1. check for common gene/alt ids, next if none ? or ignore ids and only check exon locations? # # 2. pairwise compare all exonsi x exonsj for overlap: _sameloc1 and _overlap1 ? # # if $issame for >= 2? exons, merge genes; merge issame exons; dont merge isoverlap (!same) exons # # my $mergedexons= merge_exons( $exonsi, $exonsj, 1, !$noidmatch); # ## this is unsafe: with exons from diff. predictors that can string them # ## together in widely diff. gene models, exon overlap alone wont say if these belong together # ## need exonoverlap + geneid overlap # # #warn "#D mrg[$iloop]? $ig/$geneidi x $jg/$geneidj = ",($mergedexons)?1:0,"\n" if $debugx > 0; # # if ($mergedexons && ref($mergedexons) =~ /ARRAY/) { # $imerged++; # push( @merged, [$ix, $geneidi, $mergedexons, $geneidnumi ] ); # ## but want to try more merges into merged set ... # ## may want to run thru I again using @merged # $unmerged{$ig}= $unmerged{$jg}= 0; # next NEXTI; # NEXTI; # } # } # } # # $nmerged += $imerged; # #warn "#D mrg[$iloop] in=$imerged tn=$nmerged\n" if $debugx-- > 0; # # unless($nmerged>0) { # return (wantarray) ? ($ggroup,0) : $ggroup # } # foreach my $ig (sort keys %unmerged) { # push( @merged, $mgroup->[$ig] ) if $unmerged{$ig}; # } # $mgroup= undef; $mgroup= \@merged; # } while($imerged>0); # # print $NOTEOUT "# merge_genegroup $commonid nmerge=$nmerged\n" if $debug; # return (wantarray) ? ($mgroup, $nmerged) : $mgroup; # replaces $ggroup # } =item process_genegroup ## all of ggroup is presumptive common-gene-set ## need to 4. separate out gene models, add a common-id/tag, adjust exons? ## FIXME: dont turn far, wide match regions into (gff) region? or mark differently ## so they don't command wide views .. maybe segregate gene groups by near, far types? ## but some are mixed types =cut sub process_genegroup { my($ggroup, $commonid, $exon_table, $exon_table_row, $idmap)= @_; my($gmerged, $nmerged)= merge_genegroup($ggroup, $commonid); $ggroup= $gmerged if($gmerged && $nmerged>0); # my $altexonids = $idmap->{altids}; my %geneloc=(); my %commoneq=(); my ($nexons, $naltexons, @altids, %didexons); my $igroup=0; foreach my $genegroup (@$ggroup) { my($ix, $geneid, $exons, $geneidnum)= @ $genegroup; my $g_qualscore = 0; my $g_nexons = @$exons; my $g_naltexons = 0; # count foreach my $exon (@$exons) { my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt, @altexons)= @{ $exon }; ## qid may have location: ID=Dappu1_FM5_96539.2:169647-169896 ; drop it # my($trid,$exnum)= split_exonid($qid,1); $qid= split_exonid($qid); # recode into near, far, only? ;; but see gene_quality for eqsame,eqnear,eqfar counts my $ceq= $eq; if($EQlegend{$ceq} =~ /near/) { $ceq= -2; } elsif($EQlegend{$ceq} =~ /altequal/) { $ceq= -2; } # include equal, altequal by default w/ .near elsif($EQlegend{$ceq} =~ /far/) { $ceq= 0; } $commoneq{$commonid}{$ceq}++; $commoneq{$geneid}{$ceq}++; $nexons++; if(@altexons) { my $na= scalar(@altexons); $naltexons += $na; $g_naltexons+=$na; } if($ref =~ m/[^\_][\_]{2}$/) { $naltexons++; $nexons--; } # funky flag # my $alt3= $altexonids->{$qid} || ""; #?? shouldnt need this # $alt= $alt3 if(length($alt3) > length($alt)); my @alt2= exon_allids( $exon, 0, 0); push @altids, @alt2; $alt= join(",", @alt2); foreach my $id ($commonid, $geneid) { # was $trid unless(exists $geneloc{$id}) { (my $gref = $ref) =~ s/[\_]+$//; # drop done flag $geneloc{$id} = [$gref,$tb,$te,$strand, $eq, $alt, $geneidnum, $g_nexons, $g_naltexons, # for common region $g_qualscore, #? replace this or $galt[5] with geneid[qualscore], list ]; } else { # FIXME: for far exons; not sure we want to have huge common ranges for # 'genes' and regions; use _isnear1() to decide? # _isnear2($tb, $te, $geneloc{$id}->[1], $geneloc{$id}->[2] , 2 * $NEARDIST); $geneloc{$id}->[1]= $tb if( $tb<$geneloc{$id}->[1]); $geneloc{$id}->[2]= $te if( $te>$geneloc{$id}->[2]); my $newalt= $geneloc{$id}->[5] . "," . $alt; $newalt= clean_exonlist( $newalt, 0, 0, 0); $geneloc{$id}->[5]= $newalt; } } } ## FIXME: here or elsewhere, measure each alt geneid(exons) for consistency with found gene model ## count exons-matched / total-exons ; exons-lost, exons-gained, start/stop exon? ## add as geneloc field? my($bestlist, $genequal)= gene_quality($ix, $geneid, $exons, $geneidnum); # see also $genequal->{input} scores for $exons ; or $geneid ? # change $geneid to $bestlist->[0] ?? # insert qual scores into $geneloc{$geneid}->[n] my ($baseid,$preid)= ($geneid,""); if( $baseid =~ s/^((c|g)\d+\.)//) { $preid=$1; } if($bestlist->[0] ne $baseid) { my ($oldid, $newid)= ($baseid, $bestlist->[0]); if($debug>1) { my $sco= sprintf "%.2f",($genequal->{$oldid}{score} || -1); my $scn= sprintf "%.2f",($genequal->{$newid}{score} || -1); print $NOTEOUT "# process_genegroup: quality changes id: $oldid ($sco) => $newid ($scn) \n"; } $baseid= $newid; $geneid= $preid.$baseid; exists $IDmap->{$geneid}{num} or $IDmap->{$geneid}{num}= ++$GENEID_NUM; $geneidnum= $IDmap->{$geneid}{num}; ## this was bad because we lost original geneloc{$geneid} my $oldid2= $preid.$oldid; $geneloc{$geneid}= delete $geneloc{$oldid2}; $genegroup->[1]= $geneid; $genegroup->[3]= $geneidnum; } $g_qualscore= $genequal->{$baseid}{score} || 0; $geneloc{$geneid}->[8]= $g_naltexons; $geneloc{$geneid}->[9]= $g_qualscore; # .. replace [5]/alt with new altgene list + [qualscore], or stuff as [9] my $genescores= join ",", map { my $sco= sprintf "%.0f",( 100 * $genequal->{$_}{score} || -1); my $geq= $genequal->{$_}{eqtype} || 0; my $glo= $genequal->{$_}{lost} || 0; my $gtd= $genequal->{$_}{tandems} || 0; $_."[$sco/$glo/$gtd/$geq]"; #? add more score parts? exons-lost is good } @$bestlist; $geneloc{$geneid}->[9]= $genescores; my $geneinfo= join ",", map { $_.":".$genequal->{QUALSUMID}{$_}; } sort keys %{ $genequal->{QUALSUMID} } ; $geneloc{$geneid}->[10]= $geneinfo; print $NOTEOUT "# process_genegroup info= $geneinfo \n" if $debug>1; $igroup++; } # my $eqtype= $genequal{$gid}{eqtype}; # my $exoncover= $genequal{$gid}{cover} ; # my $exgained= $genequal{$gid}{gained}; # my $exlost= $genequal{$gid}{lost} ; # my $hasends= $genequal{$gid}{ends} ; # my $tandems= $genequal{$gid}{tandems}; # my $score= $genequal{$gid}{score}; # these are gene ids, see $altexonids = $idmap->{altids}; my %altids=(); map{ # my($exnum)= (s/\.([\d]+)$//) ? $1 : 1; my($geneid,$exnum)= split_exonid($_,1); $exnum ||= 1; $altids{$geneid}{$exnum}++; $altids{$geneid}{0}++; } @altids; foreach my $id (sort keys %geneloc) { my $eqmed=0; # do this for all main gene 'match' rows? my $commoneq=1; # want median exon $eq value . foreach my $eq (sort {$b<=>$a} keys %{$commoneq{$id}}) { next if( $eq == -1 || $eq == 1 || $eq == 2); #? # recode into near, far, only? do above if($commoneq{$id}{$eq} >= $eqmed) { $commoneq= $eq; $eqmed= $commoneq{$id}{$eq}; } } $geneloc{$id}->[4]= $commoneq; } my $genelist= genelist_fromalts( \%altids ); $geneloc{$commonid}->[5]= $genelist; # == $alt list revised for common region output $geneloc{$commonid}->[6]= 0; #$geneidnum << use for $cid ?? $geneloc{$commonid}->[7]= $nexons; $geneloc{$commonid}->[8]= $naltexons; # $geneloc{$commonid}->[9]= $g_qualscore; =item MERGE genegroups ... REVISE: before print here, collect *all* genegroups and stitch together those with ... major overlaps (how? by exons or is that exhausted?) *** use %altids/$genelist common genes accross gene groups; need $nx,$kx values to decide (need $kx>1 ?) merge cid & gid (within cid?) =cut unless( COMMON_MERGE ) { print_genegroup( $ggroup, $commonid, \%geneloc); return ($ggroup, $commonid, 0, \%geneloc); #? } else { my $cid = $commonids{$commonid}; ## = $idmap->{$commonid}{num}; my @cloc = @ {$geneloc{$commonid}} [0,1,2,3] ; $geneloc{$commonid}->[6]= $cid; #$geneidnum $commonlist{$cid}{commonid}= $commonid; $commonlist{$cid}{cid}= $cid; $commonlist{$cid}{cloc}= \@cloc; $commonlist{$cid}{altids}= \%altids; $commonlist{$cid}{geneloc}= \%geneloc; $commonlist{$cid}{genegroup}= $ggroup; return ($ggroup, $commonid, $cid, \%geneloc, \%altids); #? # %{ $commonlist{$cid} } = ( # commonid => $commonid, # cid => $cid, # cloc => \@cloc, # altids => \%altids, # geneloc => \%geneloc, # genegroup => $ggroup, # ); } } sub print_genegroup { my( $ggroup, $commonid, $geneloc, )= @_; $genelocations= $geneloc; # global for _sort my $nprint= print_genegroup1_gff( $commonid, $genelocations, undef, $commonid); return $nprint if ($nprint==0); #?? flag for ignore small groups? unless(@$ggroup) { print $NOTEOUT "# missing group genes for $commonid\n"; return $nprint; } ## _geneloc1sort isnt working .. why? # my @gg= @$ggroup; my @gg= sort _geneloc1sort @$ggroup; # ? sort group items by location; need global %geneloc foreach my $gg (@gg) { my($ix, $geneid, $exons, $geneidnum)= @ $gg; $nprint += print_genegroup1_gff( $geneid, $genelocations, $exons, $commonid, $geneidnum); } return $nprint; } sub genelist_fromalts { my($altids)= @_; my @genelist= map{ my($nx)= $altids->{$_}{0}; my @kx = sort keys %{$altids->{$_}}; my($kx)= scalar(@kx) - 1; my $k2 =0; foreach my $k (@kx) { $k2++ if($k>0 && $altids->{$_}{$k}>1); } $_ ."[$nx.$kx.$k2]"; } sort keys %$altids; return (wantarray) ? @genelist : join ",", @genelist; } sub common2_overlap { my($cid1,$cid2)= @_; my $cloc1 = $commonlist{$cid1}{cloc}; my $cloc2 = $commonlist{$cid2}{cloc}; # FIXME: dont' score *far* ranges (or fix those) as overlaps to merge # when region is composed of quite distant genes or exons. ## really need a disjunct region location(s) return _isoverlapNotHuge1( @$cloc1[1,2], @$cloc2[1,2] ); ## -1 for 1 < 2 } # see genelist_fromalts sub common2_genes { my($cid1,$cid2)= @_; my $altids1= $commonlist{$cid1}{altids}; my $altids2= $commonlist{$cid2}{altids}; my $score= 0; my @id1= keys %$altids1; foreach my $id (@id1) { next unless(exists $altids2->{$id}); my($nx)= $altids2->{$id}{0}; my @kx= sort keys %{$altids2->{$id}}; my($kx)= scalar(@kx) - 1; next unless($nx>2 && $kx>1); $score++; } if($score) { my @id2= keys %$altids2; # $score= -$score if(scalar(@id1) < scalar(@id2)); # -n for 1 < 2 return (scalar(@id1) < scalar(@id2)) ? -1 : 1; } return $score; } sub common2_merge { my($cid1,$cid2)= @_; my $comid1= $commonlist{$cid1}{commonid}; my $comid2= $commonlist{$cid2}{commonid}; if(1) { my $genegroup1= $commonlist{$cid1}{genegroup}; my $genegroup2= $commonlist{$cid2}{genegroup}; push( @$genegroup1, @$genegroup2); my($exon_table, $exon_table_row)= (undef, 0); # these are not used .. drop from parlist my($new_ggroup, $new_commonid, $new_cid, $new_geneloc, $new_altids) = process_genegroup( $genegroup1, $comid1, $exon_table, $exon_table_row, $IDmap); ## ^^ updates all desired $commonlist{$cid} fields. . merged from common genegroup # # may have been changed: $commonlist{$cid1}{genegroup}= $genegroup1; # if( $commonlist{$cid1}{genegroup} ne $new_ggroup) { # warn "# common2_merge: odd change for commonlist {$cid1}: $new_commonid, $new_cid \n"; # $commonlist{$cid1}{genegroup}= $new_ggroup; # shouldnt need this # } delete $commonlist{$cid2}; } else { #..... all following steps should be handled by rerun of # # process_genegroup($ggroup, $commonid, $exon_table, $exon_table_row, $idmap) # $commonlist{$cid}{commonid}= $commonid; # $commonlist{$cid}{cid}= $cid; # $commonlist{$cid}{cloc}= \@cloc; # $commonlist{$cid}{altids}= \%altids; # $commonlist{$cid}{geneloc}= \%geneloc; # $commonlist{$cid}{genegroup}= $ggroup; my $altids1= $commonlist{$cid1}{altids}; my $altids2= $commonlist{$cid2}{altids}; my @genelist1= genelist_fromalts($altids1); my @genelist2= genelist_fromalts($altids2); foreach my $gni (@genelist2) { (my $gg=$gni) =~ s/\[.*\]//; unless( grep(/^$gg/,@genelist1) ) { push(@genelist1,$gni); foreach my $x (keys %{$altids2->{$gg}}) { $altids1->{$gg}{$x} = $altids2->{$gg}{$x}; } } } my $genelist= join(",",@genelist1); my $geneloc1= $commonlist{$cid1}{geneloc}; my $geneloc2= $commonlist{$cid2}{geneloc}; my($tb,$te)= @ { $geneloc2->{$comid2} }[1,2]; $geneloc1->{$comid1}->[1]= $tb if( $tb < $geneloc1->{$comid1}->[1]); $geneloc1->{$comid1}->[2]= $te if( $te > $geneloc1->{$comid1}->[2]); $geneloc1->{$comid1}->[5]= $genelist; $geneloc1->{$comid1}->[6]= 0; #$geneidnum my($nexons,$naltexons)= @ { $geneloc2->{$comid2} }[7,8]; $geneloc1->{$comid1}->[7] += ($nexons || 0); $geneloc1->{$comid1}->[8] += ($naltexons || 0); # also add to geneloc1 all gene locs from 2 foreach my $id (keys %{ $geneloc2 }) { next if($id eq $comid1 || $id eq $comid2); next if(exists $geneloc1->{$id}); ## FIXME merge 1,2 here; not a problem? - no gid overlap $geneloc1->{$id}= $geneloc2->{$id}; } $commonlist{$cid1}{geneloc}= $geneloc1; # done this merge my $genegroup1= $commonlist{$cid1}{genegroup}; my $genegroup2= $commonlist{$cid2}{genegroup}; push( @$genegroup1, @$genegroup2); # safe? ##? need to rerun merge_genegroup after this ; and redo process_genegroup ## to get it right .. my($gmerged, $nmerged)= merge_genegroup( $genegroup1, $comid1); $genegroup1= $gmerged if($gmerged && $nmerged>0); $commonlist{$cid1}{genegroup}= $genegroup1; my @cloc1 = @{$geneloc1->{$comid1}} [0,1,2,3]; $commonlist{$cid1}{cloc}= \@cloc1; delete $commonlist{$cid2}; } #.................. merge_genegroup replacement .......... } sub show_old_common { my($cid)= @_; print "#c.$cid# "; print "\n" and return unless exists($commonlist{$cid}); my $commonid = $commonlist{$cid}{commonid}; my $geneloc = $commonlist{$cid}{geneloc}; print_genegroup1_gff( $commonid, $geneloc, undef, $commonid); } sub commons_print { my @cids= sort{$a <=> $b} keys %commonlist; # sort by loc here instead my $nc= scalar(@cids); my $nout=0; print $NOTEOUT "# ",('.' x 50),"\n" if $debug; print $NOTEOUT "# commons_print in=$nc\n" if $debug; foreach my $cid (@cids) { my $genegroup = $commonlist{$cid}{genegroup}; my $commonid = $commonlist{$cid}{commonid}; my $geneloc = $commonlist{$cid}{geneloc}; my $np= print_genegroup( $genegroup, $commonid, $geneloc); $nout++ if($np>0); } print $NOTEOUT "# commons_print out=$nout\n" if $debug; } sub commons_merge { my @cids= sort{$a <=> $b} keys %commonlist; my $nc= scalar(@cids); print $NOTEOUT "\n\n# ",('.' x 50),"\n" if $debug; print $NOTEOUT "# commons_merge n=$nc\n" if $debug; # if($debug) { # print "# original commons\n"; # foreach my $id (@cids) { show_old_common($id); } # print "# ",('.' x 50),"\n" if $debug; # } for(my $ic=0; $ic<$nc; $ic++) { my $id= $cids[$ic]; for(my $jc=$ic+1; $jc<$nc; $jc++) { my $jd= $cids[$jc]; next unless(exists($commonlist{$jd}) && exists($commonlist{$id})); my $olap= common2_overlap($id,$jd); next unless($olap); my $ogene= common2_genes($id,$jd); next unless($ogene); my $oscore= $olap+$ogene; if($oscore<0) { print $NOTEOUT "# commons_merge $id > $jd\n" if $debug; common2_merge($jd,$id); } else { print $NOTEOUT "# commons_merge $jd > $id\n" if $debug; common2_merge($id,$jd); } } } @cids= sort{$a <=> $b} keys %commonlist; my $nm= scalar(@cids); print $NOTEOUT "# commons_merge done n=$nm\n" if $debug; # if($debug) { # print "# new commons\n"; # foreach my $id (@cids) { show_old_common($id); } # print "# ",('.' x 50),"\n\n"; # } return $nm - $nc; # ? } sub append_sameloc_exon { my( $exon, $exons, $didexon) = @_; my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; foreach my $ex (@$exons) { my $issame= _sameloc1($jref,$jb,$je,$jstrand, @$ex[0,4,5,6]); if($issame) { $didexon->{$jqid.$jb.$je}++; $exon->[0] .= '__'; # mark alt done push @{$ex}, $exon; return 1; } } return 0; } sub add_exon { my($exon, $exons, $mark, $didexon, $galt, $gids) = @_; my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; ## do before same check ? or after? my @jalt= clean_exonlist( $jalt, 0, 0, 0); foreach (@jalt) { push(@$galt, $_) unless($gids->{$_}++); } # 1st check for existing overlap to add as alt id return if( append_sameloc_exon($exon, $exons, $didexon) ); $didexon->{$jqid.$jb.$je}++; $exon->[0] .= $mark; # "__"; # mark done push( @$exons, $exon); # my @jalt= clean_exonlist( $jalt, 0, 0, 0); # foreach (@jalt) { push(@$galt, $_) unless($gids->{$_}++); } } =item exon-near remarking revise here or in calling exon_table maker, to remark or add mark for any alt-exon ids that fall in 'near' category, though main exon id may be on this spot, knowing it has also near alt-exon is important. See exon_finder at this point: my ($altnear_id,$alteq_id)=('') x 9; my $altnear=0; foreach my $aid (@alt) { my($ab,$ae); ($aid,$ab,$ae)= split(/[:-]/,$aid); next unless(defined $ae); ($ab,$ae)=($ae,$ab) if($ab>$ae); my $aeq= _isoverlap($tb,$te,$ab,$ae,$aid); if($aeq!=0) { $alteq_id=$aid; $alt =~ s/$aid/\*$aid/; # mark all matches $eq= 2 unless($eq>0 || $aeq<=0); # last; } else { if( _isnear($tb,$te,$ab,$ae,$aid)) { $altnear= 1; $altnear_id=$aid; } } } =cut sub tandem_grouper { my($exon_table, $genomefa, $queryfa, $location, $idmap)=@_; tandy_gffhead($exon_table, $genomefa, $queryfa, $location, $idmap); $IDmap= $idmap; # fixme? mostly for debug output my $sizemap = $idmap->{sizemap}; #?? delete ; want this exon{exid} size list? my $altexonids = $idmap->{altids}; #? this is idlist/exonid my $exonidmap = $idmap->{exonidmap}; # xonid <> xonid:loc my $alltable = $idmap->{alltable}; # xonid <> xonid:loc #? need this my $old_table= $exon_table; my @xsort= sort _exloc1sort @$exon_table; $exon_table= \@xsort; my $nx= scalar(@$exon_table); ## exon id / $qid *should* be unique per row ? wait, maybe not, many matches/exon ## using @list/xid doesnt seem to affect result my $etable= (ref $alltable) ? $alltable : $exon_table; my %exon_hash= (); map{ # my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt)= @{ $_ }; # my $xid= $_->[1]; # $xid = split_exonid($xid,0); # drop loc, but keep .exnum my $exon= $_; my $xidlist = join ",", @{$exon}[1,9]; my @xid = clean_exonlist( $xidlist, 0, 0, 0); foreach my $xid (@xid) { # do we get any dups this way? need some $qid.$tb didexon check? exists $exon_hash{$xid} or $exon_hash{$xid}=[]; push( @{$exon_hash{$xid}}, $exon); } } @$etable; foreach my $ix (0..$nx-1) { my $thisexon= $exon_table->[$ix]; my($ref,$qid,$qb,$qe,$tb,$te,$strand,$eq,$p,$alt)= @{ $thisexon }; # $thisexon->[7] == $eq; next if($ref =~/[\_]$/); # done this one ## keep also where eq == 1,2 but altid == -2/near ?? #? next unless($eq <= -2 || $eq == 0 || $eq >= 4); # ok? ## below step thru exon_table should find all same geneid with ++marks ## but maybe safer to also pass all here? #? next unless($alt =~ m/[^+]?(\+\+)([^,+]+)/); # only stop at first mark ?? next unless($alt =~ m/(\+{1,})([^,+]+)/); ## now have +1 entries for blat-matched genes my ($pluss,$xgid)= ($1,$2); my $whichexon = length($pluss); my ($ibefore,$iafter); #? 150, 500, 1500 if ($whichexon == 1) { ($ibefore,$iafter)= (10,500); } # dont really need -$ix elsif ($whichexon == 2) { ($ibefore,$iafter)= (100,400); } # most marks start here else { ($ibefore,$iafter)= (100,100); } # these should already be caught by mark 1 or 2 my $commonid=""; my @genegroup=(); # only collect over this marked gene? then print @genegroup my %didexon; # this tracks only for this gene group, _ mark tracks over all my %gidmark; my ($galt_ref, $gids_ref)= clean_exonlist( $alt, 0, 0, 0,",",1); ##my ($galt_ref, $gids_ref)= clean_exonlist( $alt, 0, 2, 0,",",1); # ^^ keep marks here on ids and pick ++@ (altnear) before ++* before others ? # ^^ 2 == sort @galt by marks, but drop them my @galt= @$galt_ref; my %gids= %$gids_ref; my $ixid= 0; while (my $xid= shift(@galt)) { $xid = split_exonid($xid,0); # drop loc, but keep .exnum next unless($xid =~ /\w/); $ixid++; # my($xid1,$xloc,$xeq) = split_exonid($xid,0); # drop loc, but keep .exnum # my $xid0= $exonidmap->{$xid}; # full exonid = geneid.exnum:loc=exeq ?? << no =exeq here ## need to parse $alt for =exeq value # my $xeq=-1; if($alt =~ m/$xid\:([\d\-]+)=(\d+)/) { ($xloc,$xeq)= ($1,$2); } # #xeq: expect only -2/near,0/far,1/eq my($geneid,$exnum) = split_exonid($xid,1); #? exists $IDmap->{$geneid}{num} or $IDmap->{$geneid}{num}= ++$GENEID_NUM; my $geneidnum= $idmap->{$geneid}{num}; ## should exist unless($commonid) { $COMMONID_NUM++; $commonid= $idmap->{$geneid}{common} || $geneid; $commonid= "c".$COMMONID_NUM.".".$commonid; $idmap->{$commonid}{num}= $commonids{$commonid}= $COMMONID_NUM; } my @exons=(); my @altexons= (); # 2. forestep collect all geneid exons # how many steps in array is enough? sorted by loc here (_exloc1sort) my($ix0, $ix1)= ($ix-$ibefore, $ix+$iafter); $ix0= 0 if($ix0<0); $ix1= $nx if($ix1>$nx); for (my $jx= $ix0; $jx<$ix1; $jx++) { my $exon = $exon_table->[$jx]; my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; next if (exists $didexon{$jqid.$jb.$je}); next if ($jref =~ m/[\_]$/ || $jeq == -1); #skip marked done ?? or check for alt match #?? keep distant exons from group? next unless _isnear2( $tb,$te, $jb,$je, 3 * $NEARDIST); # == 45kb ? next if ($jeq == 9); # REPEATSKIP ?? or add it anyway with flag if ( ($jx == $ix) || $jqid =~ /$geneid\.(\d+)/ || $jalt =~ /$geneid\.(\d+)/ ){ add_exon($exon, \@exons, '_', \%didexon, \@galt, \%gids); } ## add in _sameloc exons here? not below as altexons; ## will we miss some _issame due to above not yet processed ? else { my $didappend= append_sameloc_exon( $exon, \@exons, \%didexon); } } # end: step thru exon_table for marked genes # n. collect all exons (in region & near) for marked genes # ?? want this also for altids ? maybe only where $eq == 1 / equal == gene on this spot # not sure about this use here; maybe need to do AFTER all scan for -2,4,.. near exons my $exonids= $idmap->{$geneid}{exons} || []; #?? foreach my $xnum (@$exonids) { my $xid1= "$geneid.$xnum"; # if($debug) { # are we missing some? # my $nx=0; my $nh=0; # $nx= scalar(@$exonids); # exists $exon_hash{$xid1} and $nh= scalar(@{$exon_hash{$xid1}}); # print $NOTEOUT "#D $xid1 nxon=$nx; exon_hash=$nh\n" if ($nh==0); # } exists $exon_hash{$xid1} or next; #?? want even exons skipped by tandy exon_finder my @hhexons = @{$exon_hash{$xid1}}; foreach my $exon (@hhexons) { my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; next if ($jref =~ m/[\_]$/ || $jeq == -1); #mark did this one next if (exists $didexon{$jqid.$jb.$je}); # should we keep all if some exons/gene are near, but others not? unless( $eq <= -2 || $eq == 1 || $eq == 4 || $eq == 5 || _isnear1( $jb, $je, $tb, $te) ) { next; } # warn "# not near: $jqid:$jb-$je <> $xid/$qid:$tb-$te\n" ; add_exon($exon, \@exons, '_', \%didexon, \@galt, \%gids); # was "___" } } # if(0) { # # nb. altid exons that match equal should include all exons? # foreach my $kexon (@exons) { # my($kref,$kqid,$kqb,$kqe,$kb,$ke,$kstrand,$keq,$kp,$kalt, @kaltexons)= @{ $kexon }; # # next unless($keq == 1); # is this right? # # my @kalt= clean_exonlist( $kalt, 1, 0, 0); # foreach my $kgeneid (@kalt) { # my $exonids= $idmap->{$kgeneid}{exons} || []; #?? # foreach my $xnum (@$exonids) { # my $xid1= "$kgeneid.$xnum"; # exists $exon_hash{$xid1} or next; # my @exons = @{$exon_hash{$xid1}}; # foreach my $exon (@exons) { # my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; # next unless($jeq == 1); # is this right? # next if ($jref =~ m/[\_]$/ || $jeq == -1); #mark did this one # next if (exists $didexon{$jqid.$jb.$je}); # add_exon($exon, \@exons, '', \%didexon, \@galt, \%gids); # was "___" # } # } # } # } # } #..................... push @genegroup, [$ix, $geneid, \@exons, $geneidnum ] if(@exons); } # ensure main geneid exons all in group ?? other alt?? process_genegroup( \@genegroup, $commonid, $exon_table, $ix, $idmap) if @genegroup; } if( COMMON_MERGE ) { my $nreduce= commons_merge(); commons_print(); } } 1; __END__ =item Daphnia pulex tests 2007 june 14 melon:/bio/bio-grid/mb/EVidenceModeler/daphd scaffold_1 size=4193030 tandy HSP=739 ; match=161 ; region=90 scaffold_2 size=3740169 tandy HSP=1068 ; match=268 ; region=153 scaffold_3 size=3777634 tandy HSP=1281 ; match=304 ; region=153 scaffold_4 size=3075709 tandy HSP=1785 ; match=421 ; region=191 scaffold_5 size=2511979 tandy HSP=1028 ; match=172 ; region=76 scaffold_6 size=2406117 tandy HSP=1475 ; match=367 ; region=173 scaffold_7 size=2324446 tandy HSP=960 ; match=227 ; region=101 scaffold_8 size=2335496 tandy HSP=695 ; match=157 ; region=82 scaffold_9 size=2251199 tandy HSP=1259 ; match=353 ; region=201 26 MB total Gene counts cat scaffold_?/dpulex1_predict.gff | grep mRNA | perl -ne'($r,$s,$t)=split; print "$s\n" if($t eq "mRNA");' | s ort | uniq -c 11196 DGIL_SNO << this is 2x data stutter; 5598 is right 3315 JGI 4332 NCBI_GNO 2381 twinscan ?? drop these use 4332 as median gene count # genebest for Dpulex scaffold 1-9 # gene bestids scores/method (num is all ids x all gene matches, ~ 2/match) predictors num score lost tandem near equal DP_DGIL_SNO_: 3945 65.92 6.77 0.70 0.23 0.33 Dappu_FM_: 1819 63.80 2.47 0.63 0.24 0.31 NCBI_GNO_: 3110 64.18 2.93 0.88 0.26 0.30 # geneinfo for Dpulex scaffold 1-9 # gene info, all-methods scores/method (count is all gene matches) predictor count eqfar eqnear eqsame fullgns gnids maxgns methods nexons tandkno tandnew tandems DP_DGIL_SNO_ 2049 0.89 1.41 3.77 1.72 4.08 2.12 1.87 4.67 374 98 0.50 Dappu_FM_ 911 1.23 2.25 6.09 3.36 7.12 3.05 2.71 7.08 284 44 0.76 NCBI_GNO_ 1322 1.04 1.90 4.96 2.51 5.66 2.60 2.37 5.92 332 81 0.68 all 2430 0.87 1.26 3.32 1.50 3.65 1.97 1.76 4.23 387 114 0.47 # ............................................................. Using 501 for all known/new tandems, and 4332 gene predictions for these 9 scaffolds gives tandem rate of ~ 12% (is this calc right? known tandems are included in base count of 4000, but not 'new' tandem models) echo "# Tandy tandem genes stats" foreach sc (scaffold_*) echo echo "# genebest for Dpulex $sc" grep bestids= $sc/dpulex1_exons_tandy.gff | perl -n $td/genebest.perl echo "# geneinfo for Dpulex $sc" grep bestids= $sc/dpulex1_exons_tandy.gff | perl -n $td/geneinfo.perl echo "# ............................................................." end #set scs=scaffold_? #echo "# Tandy tandem genes stats for Dpulex $scs " foreach one (1) set scs="Dros.moj scaffold 6473 6541 6680" echo "# Tandy tandem genes stats for $scs " echo echo "# genebest for $scs" cat scaffold_*/*_exons_tandy.gff | grep bestids= | perl -n $td/genebest.perl echo "# geneinfo for $scs" cat scaffold_*/*_exons_tandy.gff | grep bestids= | perl -n $td/geneinfo.perl echo "# ............................................................." end =item Celegans tests ??? -- see gene-struct stats: Celegans is closer to Daphnia than Ffly, Mouse, .. -- need another organism with multi=exon genes similar to Daphnia -- 2-exon average in fflies has problems w/ some of the classifying that relies on >1 exon/gene Use this genome data, 7 chromosomes (1 is MtDNA) $sc/cele1/celegans-dna-WS167.fa.gz microbe% zgrep mRNA $sc/cele1/celegans-slim-WS167.gff.gz | perl -ne'($r,$s,$t,@x)=split; print "$s\n" if($t eq "mRNA");' | sort | uniq -c 27049 Coding_transcript << these mRNA 21182 Genefinder 221 Transposon_CDS 8641 history 21681 twinscan microbe% zgrep Coding_transcript $sc/cele1/celegans-slim-WS167.gff.gz | perl -ne'($r,$s,$t,@x)=split; print "$t\n";' | sort | uniq -c 171044 CDS << tandy data set, why more than exons? 136862 exon 17515 five_prime_UTR 20089 gene 27049 mRNA 14810 three_prime_UTR cele1: # gene info, all-methods scores/method predictor count eqfar eqnear eqsame fullgns gnids maxgns methods nexons tandkno tandnew tandems I: all 569 0.29 1.39 3.28 0.61 1.66 1.31 1.31 3.95 95 47 0.55 II: all 1253 0.27 1.43 2.71 0.63 1.63 1.33 1.26 3.57 218 148 0.60 III: all 462 0.29 1.39 3.23 0.64 1.60 1.34 1.23 3.90 85 31 0.60 IV: all 1018 0.35 1.31 2.99 0.78 1.69 1.33 1.28 3.69 163 103 0.57 V: all 2363 0.34 1.20 2.61 0.60 1.44 1.26 1.17 3.37 356 258 0.54 X: all 511 0.14 0.99 3.12 0.52 1.42 1.20 1.18 3.74 54 41 0.46 ......... total: all 6176 0.30 1.28 2.84 0.64 1.55 1.29 1.23 3.59 971 628 0.55 ......... Using 1599 for known/new tandems, and 20077 genes (6 chrs) gives tandem rate of ~ 8% cat {I,II,III,IV,V,X}/*_exons_tandy.gff | grep bestids= | perl -n $td/geneinfo.perl | grep -v '_' =item Drosophila mel tests set dpid=dmel4 cp -p $sc/dmel3/dmel_r430.fa.gz ${dpid}.fa.gz ; gunzip ${dpid}.fa.gz echo '##gff-version 3' > ${dpid}_predict.gff # Release 4.3 GFF gene model like this, use 'exon' not CDS; cant recode as CDS. gzcat $sc/dmel3/gff/dmel-*-r4.3.0.gff.gz | grep FlyBase | perl -ne \ '($r,$s,$t)=split"\t"; print if($t =~ /gene|mRNA|exon|CDS|protein|chromosome_arm/ && $s eq "FlyBase"); ' \ >> ${dpid}_predict.gff dmel4: # gene info, all-methods scores/method predictor count eqfar eqnear eqsame fullgns gnids maxgns methods nexons tandkno tandnew tandems chr2L: chr2R: chr3L: chr3R: chrX : chr4 : all 16 1.44 1.94 3.31 0.38 1.94 1.94 1.00 4.12 2 0 0.69 cat {2L,2R,3L,3R,4,X}/*_exons_tandy.gff | grep bestids= | perl -n $td/geneinfo.perl | grep -v '_' =item Drosophila moj tests melon:/bio/bio-grid/mb/EVidenceModeler/dmoj scaffold_6541/dmoj_caf060210.fa scaf-size=2543558; Sizes: scaffold_6308 size=3356042 scaffold_6328 size=4453435 scaffold_6359 size=4525533 scaffold_6473 size=16943266 tandy HSP=2371 ; match=629 ; region=301 scaffold_6482 size=2735782 scaffold_6496 size=26866924 scaffold_6498 size=3408170 scaffold_6500 size=32352404 scaffold_6540 size=34148556 scaffold_6541 size=2543558 tandy HSP=579 ; match=212 ; region=115 scaffold_6654 size=2564135 scaffold_6680 size=24764193 tandy HSP=4477 ; match=1112 ; region=498 ** scaffold_6680 takes ~12 hrs for full tandy run (blat in ~8 hr?) Gene counts (3 scaffs): gzgrep mRNA dmoj_caf060210_predict.gff.gz | egrep 'scaffold_6473|scaffold_6541|scaffold_6680' | \ perl -ne'($r,$s,$t)=split; print "$s\n" if($t eq "mRNA");' | sort | uniq -c 9737 DGIL_SNO 8146 EISE_CEX * spurious multi-transcripts/gene 10313 EISE_CGW * spurious multi-transcripts/gene 4805 GLEAN 4108 GLEANR 3826 NCBI_GNO 3862 OXFD_GPX 3441 RGUI_GID .... == ~ 4000 genes 44 MB total for 3 scaffs; @ st = 24764193 + 16943266 + 2543558 melon.% ll scaffold_6*/dmoj_caf060210_exons_tandy.gff 1547179 Jun 14 09:34 scaffold_6473/dmoj_caf060210_exons_tandy.gff 1098772 Jun 13 22:52 scaffold_6541/dmoj_caf060210_exons_tandy.gff 3400935 Jun 14 11:18 scaffold_6680/dmoj_caf060210_exons_tandy.gff melon.% grep -c bestids= scaffold_6*/dmoj_caf060210_exons_tandy.gff scaffold_6473/dmoj_caf060210_exons_tandy.gff:629 scaffold_6541/dmoj_caf060210_exons_tandy.gff:212 melon.% grep -c region scaffold_6*/dmoj_caf060210_exons_tandy.gff scaffold_6473/dmoj_caf060210_exons_tandy.gff:301 scaffold_6541/dmoj_caf060210_exons_tandy.gff:115 melon.% grep -c HSP scaffold_6*/dmoj_caf060210_exons_tandy.gff scaffold_6473/dmoj_caf060210_exons_tandy.gff:2371 scaffold_6541/dmoj_caf060210_exons_tandy.gff:579 # genebest for Dros.moj scaffold 6473 6541 6680 # gene bestids scores/method predictors num score lost tandem near equal GI_BATZ_CNA_: 1013 65.70 1.98 0.21 0.07 0.62 GI_BREN_NSC_: 1280 63.87 2.14 0.24 0.09 0.57 GI_DGIL_SNO_: 2666 68.56 1.34 0.43 0.09 0.51 GI_EISE_CEX_: 2256 75.30 1.76 1.07 0.18 0.71 GI_EISE_CGW_: 2676 76.71 1.32 2.05 0.20 0.72 GI_NCBI_GNO_: 1302 67.12 3.09 0.42 0.15 0.61 GI_PACH_GMP_: 1165 73.96 3.74 0.29 0.10 0.84 GI_RGUI_GID_: 1795 61.31 3.11 0.35 0.09 0.50 GLEAN_ : 1977 63.73 1.50 0.45 0.11 0.49 TRdmoj_: 1481 72.42 3.12 0.92 0.23 0.68 dmoj_GLEANR_: 1825 63.36 1.88 0.47 0.12 0.48 # geneinfo for Dros.moj scaffold 6473 6541 6680 # gene info, all-methods scores/method predictor count eqfar eqnear eqsame fullgns gnids maxgns methods nexons tandkno tandnew tandems GI_BATZ_CNA_ 635 0.60 2.66 7.14 10.66 22.21 4.75 9.79 7.39 234 9 0.78 GI_BREN_NSC_ 732 0.57 2.51 6.62 9.94 20.83 4.58 9.34 6.87 265 9 0.84 GI_DGIL_SNO_ 1400 0.42 1.41 4.40 5.65 12.11 3.06 5.59 4.69 274 19 0.49 GI_EISE_CEX_ 778 0.24 2.19 5.58 7.78 16.01 3.78 8.46 5.92 238 15 0.85 GI_EISE_CGW_ 685 0.16 2.56 6.36 9.32 18.71 4.28 9.77 6.63 245 17 1.09 GI_NCBI_GNO_ 831 0.56 2.25 5.87 9.29 19.49 4.47 8.74 6.16 265 15 0.83 GI_PACH_GMP_ 635 0.18 2.40 6.52 8.89 18.64 4.09 9.92 6.76 211 12 0.85 GI_RGUI_GID_ 993 0.51 1.99 5.34 8.42 17.80 4.26 7.93 5.62 278 14 0.80 GLEAN_ 928 0.51 2.16 5.68 9.23 19.37 4.54 8.67 5.94 285 18 0.88 TRdmoj_ 775 0.32 2.31 5.85 8.62 17.55 4.15 8.85 6.16 246 17 0.94 all 1953 0.41 1.12 3.47 4.58 9.95 2.78 4.76 3.80 302 27 0.47 dmoj_GLEANR_ 893 0.53 2.24 5.73 9.50 19.93 4.66 8.85 6.01 284 18 0.91 # ............................................................. Using 329 for all known/new tandems, and ~4000 gene predictions for these 3 scaffolds gives tandem rate of ~ 8% =item togff2 test ggb # #!/usr/bin/perl # $v=9; # while(){ # s/tandy\d/tandy$v/g; s/tandy\.\w+/tandy$v/g; # s/^/#/ if(m/match\t/); # s/(ID|Parent)=([^;]+);/gene "$2" ; /; # s/tclass=(\w+);/Note "$1" ; /; # s/;cid=.*$//; # s/\t[\d\.]+e-/\t/; # print; # } # # __DATA__ ##gff-version 2 #reference = scaffold_4:939037-957688 [tandy6] feature = alignment:tandy6 region:tandy6 #bgcolor = #A0A0C0 glyph = graded_segments height = 6 connector = dashed box_subparts = 0 =item rough count of dupl. genes/tandem regions cat scaf4.tandy | perl -ne \ 'if(/tandy.near/){($c)=m/cid=(\d+)/; ($g)=m/gid=(\d+)/; $c{$c}{$g} .= "$c\t$g\n" if($g);} \ elsif(/tandy.common/){ ($c)=m/cid=(\d+)/; chomp; @v=split"\t"; $g=join(";",@v[0,1,3,4,6,8]); \ $g=~s/tdclass=\w+//; $c{$c}{0}= "$c\t0\t$g\n";} \ END{foreach $c (sort keys %c){ @g=sort keys %{$c{$c}}; \ if(@g>1){ print @{$c{$c}}{@g}; }}}' \ | sort | uniq -c | sort -k2,2n -k3,3n 1 6 0 scaffold_4;tandy.common;72618;81800;+;ID=td_c.DP_DGIL_SNO_00002359;;cid=6;ids=Dappu1_FM5_18890[3,1,1];nexons=7;altexons=3 7 6 936 ^^ 4 dupl. genes, but 3 predictors mix up which are proper gene models 1 7 0 scaffold_4;tandy.common;72320;104289;+;ID=td_c.DP_DGIL_SNO_00002357;;cid=7;ids=DP_DGIL_SNO_00002358[13,2,2],DP_DGIL_SNO_00002360[20,2,2],Dappu1_FM5_95896[9,1,1],NCBI_GNO_210044[24,2,2],NCBI_GNO_214044[10,1,1],NCBI_GNO_216044[2,2,0];nexons=41;altexons=12 1 7 623 4 7 1019 .. mixup with above, maybe due to messed up preditor models (1 joins above + left other gene) 1 13 0 scaffold_4;tandy.common;285155;287317;+;ID=td_c.DP_DGIL_SNO_00002397;;cid=13;ids=Dappu1_FM5_221092[3,1,1],NCBI_GNO_256044[2,1,1];nexons=4;altexons=1 3 13 682 .. probably 2 tandem genes 1 21 0 scaffold_4;tandy.common;564347;566156;+;ID=td_c.DP_DGIL_SNO_00002467;;cid=21;ids=Dappu1_FM5_95983[3,1,1],NCBI_GNO_348044[2,1,1];nexons=5;altexons=1 2 21 1082 .. tandem pair, JGI, GNOMON get right, snap joins to 1 1 28 0 scaffold_4;tandy.common;726750;769626;-;ID=td_c.Dappu1_FM5_191599;;cid=28;ids=DP_DGIL_SNO_00002510[18,12,6],DP_DGIL_SNO_00002511[3,2,1],DP_DGIL_SNO_00002513[33,12,12],DP_DGIL_SNO_00002514[8,3,2],DP_DGIL_SNO_00002515[1,1,0],Dappu1_FM5_191599[35,10,7],Dappu1_FM5_42513[38,9,6],NCBI_GNO_396044[66,14,12],NCBI_GNO_398044[50,10,9];nexons=108;altexons=22 3 28 694 .. tandem pair probably; 1 30 0 scaffold_4;tandy.common;788933;790821;+;ID=td_c.NCBI_GNO_412044;;cid=30;ids=NCBI_GNO_412044[4,2,2];nexons=4 2 30 638 .. rev tandem pair 1 38 0 scaffold_4;tandy.common;969489;974809;+;ID=td_c.Dappu1_FM5_191628;;cid=38;ids=DP_DGIL_SNO_00002561[4,3,1],DP_DGIL_SNO_00002563[21,7,4],Dappu1_FM5_191628[15,6,3],NCBI_GNO_470044[11,5,3],NCBI_GNO_472044[9,4,1];nexons=42;altexons=6 7 38 225 1 38 443 .. Cyp tandem pair+ 1 81 0 scaffold_4;tandy.common;2366621;2390272;+;ID=td_c.DP_DGIL_SNO_00002892;;cid=81;ids=DP_DGIL_SNO_00002891[10,4,2],DP_DGIL_SNO_00002892[55,7,6],DP_DGIL_SNO_00002893[54,5,5],DP_DGIL_SNO_00002894[93,8,7],DP_DGIL_SNO_00002895[127,12,11],DP_DGIL_SNO_00002897[163,7,7],DP_DGIL_SNO_00002898[57,2,2],Dappu1_FM5_207029[84,6,6],Dappu1_FM5_221319[69,3,3],Dappu1_FM5_230334[54,1,1],Dappu1_FM5_42627[101,3,3],Dappu1_FM5_96311[27,3,3],NCBI_GNO_1418043[66,3,3],NCBI_GNO_1426043[78,3,3],NCBI_GNO_1448043[54,1,1],NCBI_GNO_914044[27,4,3],NCBI_GNO_920044[84,5,4],NCBI_GNO_924044[103,6,6],NCBI_GNO_928044[105,3,3];nexons=251;altexons=120 9 81 25 23 81 220 .. hemoglobin 8-pack 10 81 304 3 81 420 30 81 518 15 81 996 1 82 0 scaffold_4;tandy.common;2376449;2389371;+;ID=td_c.NCBI_GNO_926044;;cid=82;ids=NCBI_GNO_926044[6,3,1];nexons=16;altexons=3 7 82 608 =cut =item gene groups and altexons These 3 gene groups share exons via near-same (altexons); should they become 1 group? run a risk of making huge groups via shared exons # query: scaffold_4/dpulex1_exons.nr; genome: scaffold_4/dpulex1.fa # Location: scaffold_4:830000-1030000 scaffold_4 tandy.common match 939510 952854 . - . ID=td_c. Dappu1_FM5_191620;tdclass=common;cid=6;ids=DP_DGIL_SNO_00002554[42, 12,9],Dappu1_FM5_191620[2,2,0],NCBI_GNO_454044[3,2,1], NCBI_GNO_456044[9,4,2],NCBI_GNO_458044[13,3,2];nexons=85;altexons=48 scaffold_4 tandy.common match 944216 957688 . - . ID=td_c. DP_DGIL_SNO_00002556;tdclass=common;cid=7;ids=Dappu1_FM5_42415[1,1,0 ],NCBI_GNO_460044[3,3,0],NCBI_GNO_462044[2,2,0];nexons=30;altexons=1 scaffold_4 tandy.common match 939037 950039 . - . ID=td_c. NCBI_GNO_452044;tdclass=common;cid=8;ids=DP_DGIL_SNO_00002553[2,1,1] ,NCBI_GNO_452044[7,5,1];nexons=20 =cut 1;