#!/usr/bin/perl # package main; =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 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 to ggb gff (test) from: scaffold_4 tandy.common match 6194 6340 . + . ID=td_c1.Dappu1_FM5_96137; cid=1;ctype=far;nexons=1;ids=Dappu1_FM5_96137[1,1,0] scaffold_4 tandy.far match 6194 6340 . + . ID=td_Dappu1_FM5_96137;cid =1;gid=104;altid=Dappu1_FM5_96137 scaffold_4 tandy.equal match_part 6833 6914 1.2e-39 + . Parent=td_DP_DGIL_ SNO_00002352;cid=2;gid=1489;xid=DP_DGIL_SNO_00002342.1;xalt=DP_DGIL_SNO_00002342.4,DP_DGIL_SNO_00002352.3, DP_DGIL_SNO_00002352.6 to: scaffold_1 tandy region 131343 134630 . + . ID=td_c1.D00000023;tclass=common;c id=1;ctype=near;nexons=19;altexons=2;ids=D00000022[27,2,2],D00000023[36,5,4] scaffold_1 tandy HSP 131343 131457 1.2e-50 + . Parent=td_D00000023;tclass=near;ci d=1;gid=1020;xid=D00000023.2,D00000023.6 scaffold_1 tandy HSP 131873 131991 5.6e-51 + . Parent=td_D00000023;tclass=near;ci d=1;gid=1020;xid=D00000023.2,D00000023.6 scaffold_1 tandy HSP 132222 132336 5.3e-48 + . Parent=td_D00000023;tclass=near;ci d=1;gid=1020;xid=D00000023.2,D00000023.6 scaffold_1 tandy HSP 132691 133230 3.2e-312 + . Parent=td_D00000023;tclass =equal;cid=1;gid=1020;xid=D00000022.2,D00000022.4,D00000023.3,D00000023.7 cat scaffold_?/dpulex1_exons.nr4.gff | perl -ne \ 's/DP_DGIL\D+/D/g; s/NCBI_G\D+/G/g; s/Dappu1_FM5\D+/J/g; \ s/;xalt=/,/; s/\[(\d+),(\d+),(\d+)\]/[$1.$2.$3]/g; \ s/match_part/HSP/; \ s/tandy.common\tmatch/tandy.common\tregion/; \ if(s/tandy.(\w+)/tandy/) { $t=$1; s/;cid/;tclass=$t;cid/; } \ print unless(/tandy\tmatch/); ' \ > ! dpulex1_tandy4.gff =cut use strict; BEGIN{ warn "# FIXME: correct this perl for ref/scaffold with multi-ref input "; $| = 1; # print unbuffer } use constant SAMEBASE1 => 10; # for _sameloc, slop allowed in loca == locb our ($VERSION)= "1.4h"; use constant COMMON_MERGE => 1; our %EQlegend; our $commonidnum= 0; our %commonids=(); our %commonlist=(); our ($debug); # dang our, do we get main values here? our $BINSIZE ; our $GENEBINSIZE ; # use only where also checking gene ID our $NEARDIST; #? 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 _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; } # 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={}; 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; } } # format is NCBI_GNO_1182044.1:3009787-3009874 # sub split_exonid1 { # local $_= shift; # my $dropnum= shift or 0; # my ($exnum, $exloc)=(0,""); # s/[\+\*\@]//g; # drop marks at ^front # s/\:([\d\-\.]+)$// and $exloc=$1; # $dropnum && do { s/\.([\d]+)$// && $exnum= $1; }; # return (wantarray) ? ($_,$exnum,$exloc) : $_; # } sub clean_exonlist { my($alist, $dropnum, $keepmarks, $keeploc, $joiner, $wanthash)= @_; ##my @galt= grep /^\w/, split ",",$alist; # may be only 1 my @galt= grep( /^\w/, split(/\s*[,]\s*/, $alist)); # may be only 1 my %gids= map{ s/[\+\*\@]//g unless $keepmarks; s/\:[\d\-]+$// unless $keeploc; s/\.[\d]+$// if $dropnum; $_,1; } @galt; @galt= sort keys %gids; $joiner ||= ","; return ($wanthash) ? (\@galt, \%gids) : (wantarray) ? @galt : join($joiner, @galt); } sub print_genegroup1_gff { my( $geneid, $geneloc, $exons, $commonid, $geneidnum )=@_; return unless exists $geneloc->{$geneid}; my($gref,$gb,$ge,$gstrand, $geq, $galt, $geneidnum2, $gnexons, $gnaltexon)= @{ $geneloc->{$geneid} }; $geneidnum ||= $geneidnum2; $gref =~ s/[\_]+$//; # drop done flag my $matchID= "td_$geneid"; # FIXME: option my $gat= "ID=$matchID"; my $cid= $commonids{$commonid}; # commonid num # $commonids{$commonid}= $idmap->{$commonid}{num}= $commonidnum; my $ctype=""; (my $geqtype= $EQlegend{$geq}) =~ s/\W/_/g; if($commonid eq $geneid) { $ctype= $geqtype; $geqtype= "common"; print "\n"; } #?? $gat .= ";tdclass=$geqtype"; $gat .= ";cid=$cid"; $gat .= ";gid=$geneidnum" if($geneidnum); $gat .= ";ctype=$ctype" if($ctype); $gat .= ";nexons=$gnexons" if($gnexons); $gat .= ";altexons=$gnaltexon" if($gnaltexon); if($galt && $commonid eq $geneid) { $gat .= ";ids=$galt"; } elsif($galt) { $galt= clean_exonlist( $galt, 1, 0, 0, "."); $gat .= ";altid=$galt" if($galt); } print join("\t",$gref,"tandy.".$geqtype,"match", $gb, $ge,".",$gstrand,".",$gat),"\n"; return 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?? print_exon1( $exon, $xat); ## , $gref, $gstrand); # if(0) { # instead print as part of altid list for main exon -- same location # $xat.= ";altexon=1" unless($isalt); # foreach my $aexon (@altexons) { # print_exon1( $aexon, $xat); ## , $gref, $gstrand); # } # } } } 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($xid,$exnum)= ($qid,1); #if($qid =~ m/^(.+)\.(\d+)/) { $xid=$1; $exnum=$2; } (my $xeqtype= $EQlegend{$eq}) =~ s/\W/_/g; my $xat= $xattrib; #?? $xat .= ";tdclass=$xeqtype"; $xat .= ";xid=$qid"; # $xat.= ";ix=$exnum"; ## dont need this with xid.exnum my $alt2= $xalt; foreach my $ax (@altexons) { $alt2 .= ",".$ax->[1]; # qid $alt2 .= ",".$ax->[9]; # alt } $alt2= clean_exonlist( $alt2, 0, 0, 0); $alt2 =~ s/$qid[,]?//; $xat.= ";xalt=$alt2" if($alt2); # want this? debug only? print join("\t",$ref,"tandy.".$xeqtype,"match_part",$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 # } sub process_genegroup { my($ggroup, $commonid, $exon_table, $exon_table_row, $idmap)= @_; ## all of ggroup is presumptive common-gene-set ## need to 4. separate out gene models, add a common-id/tag, adjust exons? my $altexonids = $idmap->{altids}; my %geneloc=(); my %commoneq=(); my ($nexons, $naltexons, @altids, %didexons); foreach my $gg (@$ggroup) { my($ix, $geneid, $exons, $geneidnum)= @ $gg; 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 $qid= split_exonid($qid); my($trid,$exnum)= split_exonid($qid,1); $commoneq{$commonid}{$eq}++; $commoneq{$geneid}{$eq}++; $nexons++; if(@altexons) { $naltexons += scalar(@altexons); } if($ref =~ m/[^\_][\_]{2}$/) { $naltexons++; $nexons--; } # funky flag # my $alt3= $altexonids->{$qid} || ""; #?? shouldnt need this # $alt= $alt3 if(length($alt3) > length($alt)); my $alt2= "$qid,$alt"; # altexon_recursion( \$alt2, \@altexons, \%didexons) if(@altexons); foreach my $ax (@altexons) { $alt2 .= ",".$ax->[1]; # qid $alt2 .= ",".$ax->[9]; # alt } my @alt2= clean_exonlist( $alt2, 0, 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]; } else { $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; } } # #?? handle @altexons # foreach my $id ($trid) { # ?? , $trid # unless(exists $genes{$id}) { $genes{$id}=[]; } # push( @{ $genes{$id}}, [$exnum, $qid, $ref, $tb, $te, $p, $eq, $alt]); # ## FIXME: ^ with common ids, can have many exons at same location ; test # } } } # these are gene ids, see other altids == exon ids 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); #? if($commoneq{$id}{$eq} >= $eqmed) { $commoneq= $eq; $eqmed= $commoneq{$id}{$eq}; } } $geneloc{$id}->[4]= $commoneq; } my $genelist= genelist_fromalts( \%altids ); $geneloc{$commonid}->[5]= $genelist; $geneloc{$commonid}->[6]= 0; #$geneidnum $geneloc{$commonid}->[7]= $nexons; $geneloc{$commonid}->[8]= $naltexons; =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); } else { my $cid = $commonids{$commonid}; ## = $idmap->{$commonid}{num}; my @cloc = @ {$geneloc{$commonid}} [0,1,2,3] ; $commonlist{$cid}{commonid}= $commonid; $commonlist{$cid}{cid}= $cid; $commonlist{$cid}{cloc}= \@cloc; $commonlist{$cid}{altids}= \%altids; $commonlist{$cid}{geneloc}= \%geneloc; $commonlist{$cid}{genegroup}= $ggroup; # %{ $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 print_genegroup1_gff( $commonid, $genelocations, undef, $commonid); ## _geneloc1sort isnt working .. why? my @gg= sort _geneloc1sort @$ggroup; # ? sort group items by location; need global %geneloc # my @gg= @$ggroup; foreach my $gg (@gg) { my($ix, $geneid, $exons, $geneidnum)= @ $gg; print_genegroup1_gff( $geneid, $genelocations, $exons, $commonid, $geneidnum); } } 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 _isoverlap2 { 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 common2_overlap { my($cid1,$cid2)= @_; my $cloc1 = $commonlist{$cid1}{cloc}; my $cloc2 = $commonlist{$cid2}{cloc}; ##return _isoverlap2( $cloc1->[1], $cloc1->[2], $cloc2->[1],$cloc2->[2] ); ## -1 for 1 < 2 return _isoverlap2( @$cloc1[1,2], @$cloc2[1,2] ); ## -1 for 1 < 2 # also return -1 for cloc1 > cloc2, +1 for cloc1 < cloc2 ?? } # 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)= @_; # $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 $comid1= $commonlist{$cid1}{commonid}; my $comid2= $commonlist{$cid2}{commonid}; 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? $commonlist{$cid1}{genegroup}=$genegroup1; # foreach my $gg (@$g1group, @$g2group) { # my($ix, $geneid, $exons, $geneidnum)= @ $gg; # } my @cloc1 = @{$geneloc1->{$comid1}} [0,1,2,3]; $commonlist{$cid1}{cloc}= \@cloc1; delete $commonlist{$cid2}; } 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); warn "# ",('.' x 50),"\n" if $debug; warn "# commons_print n=$nc\n" if $debug; print "# commons_print n=$nc\n" unless $debug; #?? debug foreach my $cid (@cids) { my $genegroup = $commonlist{$cid}{genegroup}; my $commonid = $commonlist{$cid}{commonid}; my $geneloc = $commonlist{$cid}{geneloc}; print_genegroup( $genegroup, $commonid, $geneloc); } } sub commons_merge { my @cids= sort{$a <=> $b} keys %commonlist; my $nc= scalar(@cids); warn "\n\n# ",('.' x 50),"\n" if $debug; warn "# 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) { warn "# commons_merge $id > $jd\n" if $debug; common2_merge($jd,$id); } else { warn "# commons_merge $jd > $id\n" if $debug; common2_merge($id,$jd); } } } @cids= sort{$a <=> $b} keys %commonlist; my $nm= scalar(@cids); warn "# 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; # $nm - $nc ? } sub add_exon { my($exon, $exons, $mark, $didexon, $galt, $gids) = @_; my($jref,$jqid,$jqb,$jqe,$jb,$je,$jstrand,$jeq,$jp,$jalt)= @$exon; # return if ($jref =~/[\_]$/ || $jeq == -1); #mark did this one # return if (exists $didexon->{$jqid.$jb.$je}); ## 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 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; #? do we include alt exon's alt id list in galt ? return; # last; } } $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->{$_}++); } } sub tandem_grouper { my($exon_table, $genomefa, $queryfa, $location, $idmap)=@_; tandy_gffhead($exon_table, $genomefa, $queryfa, $location, $idmap); my $sizemap = $idmap->{sizemap}; #?? delete ; want this exon{exid} size list? my $altexonids = $idmap->{altids}; #? this is idlist/exonid #? 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 %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 # $xid =~ s/\:[\d\-]+$//; # drop :location; keep .exnum exists $exon_hash{$xid} or $exon_hash{$xid}=[]; push( @{$exon_hash{$xid}}, $_); } @$exon_table; 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 unless($eq == -2 || $eq == 0 || $eq == 4); # ok? next if($ref =~/[\_]$/); #mark did this one next unless($alt =~ m/[^+]?(\+\+)([^,+]+)/); # only stop at first mark my $xgid= $1.$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= @$galt_ref; my %gids= %$gids_ref; # my @galt= split ",",$alt ; # these are alternate exon ids, not gene ids # my %gids= map{ # s/[\+\*\@]//g; s/\:[\d\-]+$//; # drop loc, but keep .exnum # $_,1; # } @galt; while (my $xid= shift(@galt)) { # $xid =~ s/[\+\*\@]//g; $xid =~ s/\:[\d\-]+$//; # drop :location; keep .exnum $xid = split_exonid($xid,0); # drop loc, but keep .exnum next unless($xid =~ /^\w/); my($geneid,$exnum) = split_exonid($xid,1); #my($geneid,$exnum)= ($xid,1); #if($xid =~ m/^(.+)\.(\d+)/) { ($geneid,$exnum)= ($1,$2); } my $geneidnum= $idmap->{$geneid}{num} || 0; unless($commonid) { $commonidnum++; $commonid= $idmap->{$geneid}{common} || $geneid; $commonid= "c".$commonidnum.".".$commonid; $idmap->{$commonid}{num}= $commonids{$commonid}= $commonidnum; } my @exons=(); my @altexons= (); # 2. forestep collect all geneid exons ## my($ix0, $ix1)= ($ix-150, $ix+150); my($ix0, $ix1)= ($ix-500, $ix+500); $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 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; ? need flag ## will we miss some _issame due to above not yet processed ? else { if(1) { # do this check also in add_exon foreach my $ex (@exons) { my $issame= _sameloc1($jref,$jb,$je,$jstrand, @$ex[0,4,5,6]); # what of big exons that overlap all of smaller (and more ..) ? if($issame) { $didexon{$jqid.$jb.$je}++; $exon->[0].='__'; # mark alt done push @{$ex}, $exon; last; } } } #................ } } # 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"; 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 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 == 4 || $eq == 1 || _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 "___" # } # } # } # } # } #..................... # 3. exons at same/almost same locs # ** need these to stitch gene groups together ... better than this ## ? put these into @exons instead of each exon tail? ## do this check also in add_exon if(0) { 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 ($jref =~ m/[\_]$/ || $jeq == -1); #mark did this one next if (exists $didexon{$jqid.$jb.$je}); foreach my $ex (@exons,@altexons) { my $issame= _sameloc1($jref,$jb,$je,$jstrand, @$ex[0,4,5,6]); # what of big exons that overlap all of smaller (and more ..) ? if($issame) { $didexon{$jqid.$jb.$je}++; $exon->[0].='__'; # mark alt done push @{$ex}, $exon; push @altexons, $exon; last; } } } } #............................... push @genegroup, [$ix, $geneid, \@exons, $geneidnum ] if(@exons); } # ensure main geneid exons all in group ?? other alt?? # do we need this check ? # if (0 && @genegroup) { # ## @galt= map { $_ if($gids{$_}>2) } sort keys %gids; # @galt= grep { ($gids{$_}>2) } sort keys %gids; # while (my $xid= shift(@galt)) # { # my($geneid,$exnum)= ($xid,1); # if($xid =~ m/^(.+)\.(\d+)/) { ($geneid,$exnum)= ($1,$2); } # my @exons=(); my @altexons= (); # my $geneidnum= $idmap->{$geneid}{num} || 0; # # # 2. forestep collect all geneid exons # my($ix0, $ix1)= ($ix-500, $ix+500); # $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($jref =~/[\_]$/ || $jeq == -1); #mark did this one # # if ( # ($jqid =~ /$geneid\.(\d+)/) # main ids only here # ){ # # $exon->[0].="_"; # DONT mark done # next if($didexon{$jqid.$jb.$je}++ > 0); # push( @exons, $exon); # } # } # push @genegroup, [$ix, $geneid, \@exons, $geneidnum ] if(@exons); # } # } process_genegroup( \@genegroup, $commonid, $exon_table, $ix, $idmap) if @genegroup; } if( COMMON_MERGE ) { commons_merge(); commons_print(); } } 1; __END__ =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;