#!/usr/bin/perl # tandyx.perl =item notes - see similar genepairs.perl + add filter alt-tr: same predictor, same exon-loc (any case of shared CDS-exon == alt-tr?) * need also to separate full gene-exon set vs single exon dupls from tandy scores - not in HSP, but match, as ? fullgene:1, other flag =item usage setenv dp dmoj # pick tandy exon dupl w/ a nearby (-2) score # ?? add in match features; gzcat ${dp}1/${dp}_exons_tandy6jmd.gff.gz \ | egrep -v 'tandyweak|terepeat' | grep HSP | grep '=-2' \ > ${dp}1/${dp}_tandy6jmd_hspnear.gff # hspnear2: collect match !eqnear:0, collect all hsp/match gzcat ${dp}1/${dp}_exons_tandy6jmd.gff.gz | egrep -v 'tandyweak|terepeat|region|eqnear:0' | perl -ne \ 'if(/\tmatch/){ ($d)=m/ID=(\w+)/; print;} elsif(/\tHSP/){ ($p)=m/Parent=(\w+)/; print if($p eq $d); }'\ > ${dp}1/${dp}_tandy6jmd_hspnear2.gff set dspp=( dsec dsim dyak dere dana dpse dwil dvir dmoj dgri ) foreach sp ($dspp) echo =========[ $sp ]=============== setenv dp $sp ; cat ${dp}1/${dp}_tandy6jmd_hspnear.gff | perl $td/tandyx.perl echo end =item in data : xid= skip 1st when follow same id =eqval : want id cases w/o =eqval, they are identical exons to others collapsed in tandy .nr exons : but need to know which exon they match .. only in scaffold_nnn/*_exons.nr.gz scaffold_14624 tandy.near HSP 9538 9584 8.9e-18 - . Parent=td_s4c1g0.GLEAN_23726; tclass=near;cid=1;gid=5333;gxclass=g0j2i2; xid=GH_RGUI_GID_mRNA_40015172.7,GH_RGUI_GID_mRNA_40015172.7:15767-15826-=-2 scaffold_14624 tandy.near HSP 820631 820954 1.8e-183 - . Parent=td_s4c34g11.GH _RGUI_GID_mRNA_40015218;tclass=near;cid=34;gid=4025;gxclass=g11j0bei24; xid=GLEAN_23766.1,GH_BREN_NSC_50037609.1:823502-824017=-2, GH_BATZ_CNA_44034824.1:823502-824017-=-2,GH_DGIL_SNO_28256425.2:834714-834872, GH_BATZ_CNA_44034826.1:834259-834876,GH_BREN_NSC_50037610.1:827246-827827, GH_BREN_NSC_50037611.1:834259-834876,GH_DGIL_SNO_28256422.1:823502-823862,GH_RGUI_GID_mRNA_40015218.1:823502-823967, GH_RGUI_GID_mRNA_40015218.2:827196-827830,GH_RGUI_GID_mRNA_40015219.1:834271-834879, GLEAN_23766.1:823502-824020,dgri_GLEANR_14085.1:823502-824020, GH_BATZ_CNA_44034825.1:827246-827576,GH_DGIL_SNO_28256423.1:827246-827576,GLEAN_23767.1:827246-827576, GLEAN_23768.1:834259-834879,dgri_GLEANR_14086.1:827246-827576, dgri_GLEANR_14087.1:834259-834879 =cut use strict; use warnings; use Getopt::Long; my $align= 20; my( %gsum, %genes, %alttr, %ncalc, %in, %ggroup, %glocs, %gflags); my $criteria=""; my $pctgene=0; my $BYGENE=0; my $DUMPIDS=0; my $MINHSP=1; my $optok= GetOptions( "pctgene!", \$pctgene, "BYGENE!", \$BYGENE, "MINHSP=i", \$MINHSP, "dumpids!", \$DUMPIDS, ); die "usage: tandyx [-pctgene -bygene -minhsp=$MINHSP -dumpids] < exons_tandy6jmd.gff" unless($optok); while(<>) { next unless(/^\w/); my @v=split; my($tref,$src,$typ,$tb,$te)= @v[0..4]; my($terepeat)= ( m/terepeat=(\w+)/ ) ? $1 : 0; if( m/info=([^;]+)/) { # == $typ eq "match" my $in=$1; my @in=split",",$in; %in= map { split ":" }@in; } next unless($typ =~ /^(HSP|match)/); next if ($in{lowquality} or $terepeat or $src =~ /tandyweak/); $criteria .= " not(lowquality or TErepeat) and" unless($criteria =~ /lowqual/); next unless( $in{pctalign}>$align ); # $in{fullgenes} > 0 or $criteria .= " (pctalign>$align)" unless($criteria =~ /pctalign/); ## note lowquality includes fullgenes:0 my @at=split";",$v[-1]; my %at=map{my($k,$v)=split"=",$_,2;$k,$v;}@at; if($typ eq "HSP" and $at{xid}) { my @xid=split",",$at{xid}; my %gx=(); my %didex; my %sz=(); my $sz= $v[4] - $v[3]; ## use this odd key from tandy: gxclass=g1j5bei5; ## 'b' = begin-exon; 'e' = end-exon; if gene-hsp set has both, flag as fullgene ## other fullgenes: flag is tandem-cds-length / source cds-length my $gxclass= $at{gxclass} || ""; ##foreach (grep /\d=/, @xid) foreach (@xid) { my($gr,$gn,$xn,$eq,$xb,$xe); #unless( ($gr,$gn,$xn,$eq)= m/(\D+)(\d+)\.(\d+):[^=]+=([\d-]+)/ ) unless( ($gr,$gn,$xn,$xb,$xe,$eq)= m/(\D+)(\w+)\.(\d+)\:(\d+)\-(\d+)[^=]*=([\d-]+)/ ) { # calculate $eq =-2 if near; .. =1 if overlap; .. if( ($gr,$gn,$xn,$xb,$xe)= m/(\D+)(\w+)\.(\d+):(\d+)\-(\d+)/ ) { my($tb,$te)= @v[3,4]; if($tb <= $xe && $te >= $xb){ $eq=1; } # _isoverlap() else { my $bd= abs($tb - $xe); my $ed= abs($xb - $te); if($bd < 15000 || $ed < 15000) { $eq=-2; } else { $eq=0; } # far } $ncalc{$gr}++; } } # if( $gr=~/nada|RGUI|GLEAN_|BATZ|CEX|TR/) if($gn and $gr !~ m/nada|GLEAN_/){ next if($didex{"$gr.$gn.$xn.$eq"}++); $gx{$gr}{$eq}{$xn}++; $sz{$gr}{$eq} += $sz; $genes{$gr}{$gn}++; ## convert to genepairs hash:merge w/ that, using this geneid.exon count my $scid= $gr.$gn; ## this == gene id; ignoring tandy cluster id now. my $cxnum= $xn; my $exeq= $eq; $ggroup{$gr}{$scid}{$cxnum}{$exeq}++; $ggroup{"all"}{$scid}{$cxnum}{$exeq}++; #?? $ggroup{$gr}{"genes"}{$scid}++; # instead of $scgid use ggroup geneid $glocs{$gr}{$scid}{$cxnum}{$exeq} .= "$tref:$tb-$te,"; $gflags{$gr}{"start"}{$scid}{$exeq}++ if($gxclass =~ /b/); $gflags{$gr}{"stop"}{$scid}{$exeq}++ if($gxclass =~ /e/); # [$tref,$tb,$te]; # this will overwrite other near matches } }; ## are these gsums ok? this is still in 1-HSP loop; ns,nt are #exons/{group,eqclass,HSP} foreach my $gr (sort keys %gx) { my $ns= scalar(keys %{ $gx{$gr}{1} })||0; $gsum{$gr}{same} += $ns; my $nt= scalar(keys %{ $gx{$gr}{-2} })||0; $gsum{$gr}{near} += $nt; my($szsa,$szne); $gsum{$gr}{szsame} += $szsa= $sz{$gr}{1} || 0; $gsum{$gr}{sznear} += $szne= $sz{$gr}{-2} || 0; $gsum{$gr}{tsame}++ if ($ns>1); # was >0 .. cmp to tandem > 1; BIG Drop (>0: 5k,10k; >1: 50,400) ## if($nt>1 and $nt>$ns) { $gsum{$gr}{tandem}++; } #? is this useful if($nt>1 and $ns==0) { $gsum{$gr}{tandem}++; } #? tandem-only?? is this useful ## do per-gene ratios here, not below ?? $gsum{$gr}{ne_sa} += ( $nt / ($ns || 1) ); ###$gsum{$gr}{Sne_sa} += ($szne / ($szsa || 1)); #? ratio or diff near - same ? if($szne and $szsa) { # only when both exist # $gsum{$gr}{Sne_sa} += ($szne - $szsa); # $gsum{$gr}{NSne_sa} += $nt/$ns; # count for equal n exons $gsum{$gr}{Sne_sa} += ($szne / $szsa); $gsum{$gr}{NSne_sa} += 1; } } } } print_gsum(); print_groups(); # try both counts; which best? print_geneids() if($DUMPIDS); sub print_groups { my $species= $ENV{dp} || $ENV{species}; my $cc="# "; my $flag=""; $flag.= ($BYGENE) ? "Gene counts; " : "Exon counts; "; $flag.= "min_HSP=$MINHSP; "; $flag.= "species=$species; " if($species); print "\n"; print $cc,"Tandem gene clusters, options: $flag\n"; print $cc,"Criteria: $criteria\n"; # my @colv=qw(exons samenear same nearother near far); # my @col1=qw(N same+ same- near- near- far-); # my @col2=qw(exons near only other only only); my @colv=qw(exons samenear same near far sameall nearall nearfull); # my @col1=qw(N same+ same- near- far- same- near- near-); # my @col2=qw(exons near only only only all all full); my @col1=qw(N same+ only- only- only- same- near- near-); my @col2=qw(exons sanear same near far saall neall nefull); if($pctgene) { $col1[0]= "Exons"; $col2[0]="/gene"; } my $w; $w=6; print $cc; map{ printf "%".$w."s\t",$_; $w=6; } qw( _ _ N),@col1; print"\n"; ##$w=8; print $cc; map{ printf "%".$w."s\t",$_; $w=6; } qw(Species Group genes),@col2; print"\n"; $w=6; print $cc; map{ printf "%".$w."s\t",$_; $w=6; } qw(species group genes),@col2; print"\n"; $cc=""; # "#x "; my @ggroup= grep{$_ ne "all"} sort keys %ggroup; foreach my $ggroup (@ggroup) { my %cscore=(); my $ngene= scalar( keys %{ $ggroup{$ggroup}{"genes"} } ); my @clust= grep{ $_ ne "genes" } sort keys %{ $ggroup{$ggroup} }; my $nclus= @clust; # this == gene ids now for tandyx foreach my $scid (@clust) { my @cxnum= sort keys %{ $ggroup{$ggroup}{$scid} }; my $nx= @cxnum; my %score=(); if($BYGENE) { my %xscore=(); foreach my $cxnum (@cxnum) { my $eqsame= $ggroup{$ggroup}{$scid}{$cxnum}{1} || 0; my $eqnear= $ggroup{$ggroup}{$scid}{$cxnum}{-2} || 0; my $eqfar = $ggroup{$ggroup}{$scid}{$cxnum}{0} || 0; my $val=""; $val.= "same" if($eqsame); $val.= "near" if($eqnear); $val.= "far" if($eqfar && !$val); $xscore{$val}++; # << "samenear" score misleading?, 1/2 of count for same,near separately $xscore{"nearall"}++ if($eqnear); $xscore{"sameall"}++ if($eqsame); } # genescore++ for any exonscore >0 or >1 ? foreach my $sk (keys %xscore) { $score{$sk}++ if( $xscore{$sk} >= $MINHSP ); } $score{nearfull}++ if($gflags{$ggroup}{"start"}{$scid}{-2} and $gflags{$ggroup}{"stop"}{$scid}{-2}); } else { foreach my $cxnum (@cxnum) { my $eqsame= $ggroup{$ggroup}{$scid}{$cxnum}{1} || 0; my $eqnear= $ggroup{$ggroup}{$scid}{$cxnum}{-2} || 0; my $eqfar = $ggroup{$ggroup}{$scid}{$cxnum}{0} || 0; my $val=""; $val.= "same" if($eqsame); $val.= "near" if($eqnear); $val.= "far" if($eqfar && !$val); $score{$val}++; # << "samenear" score misleading?, 1/2 of count for same,near separately $score{"nearall"}++ if($eqnear); $score{"sameall"}++ if($eqsame); if( 0 && $val eq "near") { $eqsame= $ggroup{"all"}{$scid}{$cxnum}{1} || 0; $score{"nearother"}++ and $score{$val}-- if($eqsame); # && $score{$val}--; } } if($MINHSP > 1) { foreach my $sk (keys %score) { delete $score{$sk} if( $score{$sk} < $MINHSP ); } } $score{nearfull}++ if($gflags{$ggroup}{"start"}{$scid}{-2} and $gflags{$ggroup}{"stop"}{$scid}{-2}); } #? $cscore{"gdup"}++ if($score{"nearall"}); $cscore{"exons"} += $nx; foreach my $sk (keys %score) { $cscore{$sk} += $score{$sk} ; } } # $ggroup =~ s/_$//; $ggroup=substr($ggroup,0,14); my ($sp,$grp)= clean_group($species,$ggroup); print $cc; printf "%6s\t%10s\t",$sp,$grp; # printf "%15s\t%6d\t%6d\t",$ggroup,$nclus,$ngene; # nclus == ngene printf "%6d\t",$ngene; # nclus == ngene foreach my $sk (@colv) { if($pctgene) { printf "%6.3f\t",($cscore{$sk} || 0) / $nclus; } else { printf "%6d\t",($cscore{$sk} || 0); } } print "\n"; } } sub clean_group { my($species,$ggroup)=@_; my ($sp,$grp); unless( ($sp,$grp)= $ggroup =~ m/^(G.|CG|d...)_(\D+)/ ) { if($ggroup =~ m/^TR(d...)_/) { $sp=$1; $grp="OXFD_GPX"; } else { $grp=$ggroup; $sp=$species; } } $sp= $species if($sp =~ /^(CG|G.)$/ and $species); $grp =~ s/RGUI_GID_mRNA/RGUI_GID/; $grp =~ s/_$//; $grp= substr($grp,0,12); return ($sp,$grp); } sub print_geneids { print "# GeneID near exon matches\n"; print "# ",join("\t",qw(Group GeneID nsame nnear nearfull nearloc)),"\n"; #exon_same exon_near near_only #?? add tandy near locations to check my @ggroup= grep{$_ ne "all"} sort keys %ggroup; @ggroup = sort{ ($a =~ /NCBI/ and $b !~ /NCBI/) ? -1 : ($a !~ /NCBI/ and $b =~ /NCBI/) ? +1 : $a cmp $b } @ggroup; # put NCBI_GNO first foreach my $ggroup (@ggroup) { my @geneids= grep{ $_ ne "genes" } sort keys %{ $ggroup{$ggroup} }; my (%dumpid, %dumploc); foreach my $geneid (@geneids) { my @cxnum= sort keys %{ $ggroup{$ggroup}{$geneid} }; my ($nref,$nb,$ne)=("",0,0); foreach my $cx (@cxnum) { my $eqsame= $ggroup{$ggroup}{$geneid}{$cx}{1} || 0; my $eqnear= $ggroup{$ggroup}{$geneid}{$cx}{-2} || 0; # my $eqfar = $ggroup{$ggroup}{$geneid}{$cx}{0} || 0; $dumpid{$geneid}{-2} += $eqnear; $dumpid{$geneid}{1} += $eqsame; if($eqnear and $glocs{$ggroup}{$geneid}{$cx}{-2}) { my $gloc= $glocs{$ggroup}{$geneid}{$cx}{-2}; my @gloc= split ",",$gloc; foreach my $gex (@gloc) { my($tref,$tb,$te)= split(/[:-]/,$gex); # split ",", ... for many? $nref=$tref; $nb= $tb if(!$nb or $tb<$nb); $ne= $te if($te>$ne); } } } $dumpid{$geneid}{nearfull}++ if($gflags{$ggroup}{"start"}{$geneid}{-2} and $gflags{$ggroup}{"stop"}{$geneid}{-2}); my $nloc= "$nref:$nb-$ne"; $dumploc{$geneid}= $nloc if($nref and $ne); } print "# group $ggroup\n"; foreach my $geneid (sort keys %dumpid) { my $cnear= $dumpid{$geneid}{-2} || 0; my $csame= $dumpid{$geneid}{1} || 0; my $nearfull= $dumpid{$geneid}{nearfull} || 0; #my $xsame= $dumpid{$geneid}{xsame} || 0; #my $xnonly= $dumpid{$geneid}{xnearonly} || 0; if($cnear) { ## ||$csame my $loc= $dumploc{$geneid} || ""; print "#g ",join("\t",$geneid,$csame,$cnear, $nearfull, $loc),"\n"; #,$xsame,$xnear,$xnonly } } print "# ------------------------------------\n\n"; } } sub print_gsum { my $species= $ENV{dp} || $ENV{species}; my $cc="#s "; printf "# %6s\t%8s\t","species","source"; print "# ",join("\t", qw(Ngene near same sznear szsame tandem tsame ne/sa ta/ts Sne/sa)),"\n"; foreach my $gr (sort keys %gsum) { my $v; # my ($sp,$grp); # unless( ($sp,$grp)= $gr =~ m/^(G.|d...)_(\D+)/ ) { # if($gr =~ m/^TR(d...)_/) { $sp=$1; $grp="OXFD_GPX"; } # else { $grp=$gr; $sp=$species; } # } # $sp= $species if($sp =~ /^G.$/ and $species); # $grp =~ s/RGUI_GID_mRNA/RGUI_GID/; # $grp =~ s/_$//; $grp= substr($grp,0,12); my ($sp,$grp)= clean_group($species,$ggroup); print $cc; printf "%6s\t%10s\t",$sp,$grp; my $ngene= $gsum{$gr}{Ngene}= scalar( keys %{$genes{$gr}} ); $gsum{$gr}{szsame}= int($gsum{$gr}{szsame} / $ngene); $gsum{$gr}{sznear}= int($gsum{$gr}{sznear} / $ngene); foreach my $t (grep !/ne_sa|Sne_sa/, sort keys %{$gsum{$gr}}) { $v=$gsum{$gr}{$t}; print "$v\t"; } if(1) { $v= $gsum{$gr}{ne_sa} / $ngene; printf "%.2f\t",$v; $v= $gsum{$gr}{tandem} / $gsum{$gr}{tsame}; printf "%.3f\t",$v; $v= $gsum{$gr}{Sne_sa} / $gsum{$gr}{NSne_sa} ; printf "%.2f\t",$v; } else { $v= $gsum{$gr}{near} / $gsum{$gr}{same}; printf "%.2f\t",$v; $v= $gsum{$gr}{tandem} / $gsum{$gr}{tsame}; printf "%.2f\t",$v; $v= $gsum{$gr}{sznear} / $gsum{$gr}{szsame}; printf "%.2f\t",$v; } print "\n"; } } __END__ # ## not alttr found this way; what are all the extra cex,oxfd near counts? # if($gn and $gr !~ m/nada|GLEAN_/){ # if($alttr{$gr} and $alttr{$gr}{$xb.$xe} and not $alttr{$gr}{$xb.$xe}{$gn}) { # $alttr{$gr}{alttr}{$xb.$xe}{$gn}++; # $gr="nada"; next; # } # else { # $alttr{$gr}{$xb.$xe}{$gn}++; # 1st tr for this gene; keep it # } # } # foreach my $gr (sort keys %alttr) { # my @xk= keys %{ $alttr{$gr}{alttr} }; # my %atr= (); # foreach my $xk (@xk) { map { $atr{$_}++ } keys %{ $alttr{$gr}{alttr}{$xk} }; } # my @atr= sort keys %atr; # # @xk= grep !/alttr/, keys %{ $alttr{$gr} }; # my %agn= (); # foreach my $xk (@xk) { map { $agn{$_}++ } keys %{ $alttr{$gr}{$xk} }; } # my $nagn= scalar(keys %agn); # my $natr= scalar(@atr); # my $nerr=0; map{ $nerr++ if($agn{$_}); } @atr; # # printf "%6s\t%8s\t",$species,$gr; # print "ngene=$nagn; nalt-tr=$natr; nalt-err:$nerr; neq-calc:$ncalc{$gr}\n"; # } # debug =item R plots #...... R plots : all stats separate; one main plot? #...... line plots 'stead of bars? X=species; w/ all pred in one plot setwd("~/Desktop/dspp-work/dspp2work/gene-duplicates/") tandyx <- "~/Desktop/dspp-work/dspp2work/gene-duplicates/tandyx2.tab" df <- txt2 <- read.table(tandyx, header=T,strip.white=T) lsp <- levels(df$species) lpd <- levels(df$source) olpd <- c(1,2,5,4,3,6) # abinit; glnr; mappers olpd <- c(1,2,3,10,7, 6, 5,4,8,9) # abinit; glnr; mappers lpd <- lpd[olpd] j<-3 # 3..12 jlab<- colnames(df)[j] # ylab sublab <- paste(", Stat:",jlab) nd <- c() for (ipd in lpd) { dfseti <- df$source == ipd ndi <- matrix( df[dfseti,j],ncol=1) rownames(ndi) <- df[dfseti,1] colnames(ndi) <- ipd nd <- cbind( nd, ndi ) } cols <- heat.colors(nrow(nd)) # diff color range set / predictor barplot(nd,beside=T,col=cols, main = paste("Predictor x Tandem gene stats",sublab), ylab = jlab, xlab = "Predictor" ) legend( "topright", rownames(nd), horiz=T, col=cols, lwd = c(2,2), bty="n" ) #.... lines statplot( df, lpd, 3) # j 3..12 statplot <- function (df, lpd, j, ylims=NULL) { jlab<- colnames(df)[j] # ylab sublab <- paste(", Stat:",jlab) nd <- c() for (ipd in lpd) { dfseti <- df$source == ipd ndi <- matrix( df[dfseti,j],ncol=1) rownames(ndi) <- df[dfseti,1] colnames(ndi) <- ipd nd <- cbind( nd, ndi ) } nx <- factor(rownames(nd),levels=rownames(nd),ordered=T) cols <- rainbow(length(lpd)) pchs <- 1:length(lpd) # ncol(nd) if(is.null(ylims)) { mi<-min(nd); mx<-max(nd); ylims <- c( mi- sign(mi)*mi/10, mx+sign(mx)*mx/10) } plot( nd[nx,1], type="n", ylim = ylims, main = paste("Predictor x Tandem genes",sublab), ylab = jlab, xlab = "Species", xaxt="n" ) axis(1,at=nx,labels=levels(nx)) i<-0; for(ipd in lpd) { i<-i+1; lines( nd[nx,ipd], type="b", col=cols[i], pch=pchs[i])} legend( "topleft", lpd, horiz=F, col=cols, pch=pchs, cex=0.9, lwd = c(1,1), bty="n" ) } #... layout(mat, widths = rep(1, ncol(mat)), heights = rep(1, nrow(mat)), respect = FALSE) =cut