#!/usr/bin/perl # genepairs.perl =item usage cat scaffold.tandy.gff | perl genepairs.perl count tandy regions with near dups, find association with predictor groups does predictor1,2,3 find/miss dupl. exons? options: -align=$align -[no]genecount -[no]skipoverlap -keepgroup=DGIL -skipgroup=DGIL -marked=proteinhsp|estmatch .. keep only tandy.gff with this attribute =item tandyx update * see related tandyx.perl: combine any parts from that? + size stats; =item batch #!/bin/tcsh set td=/bio/bio-grid/dpulex/prots/ set tdgff=`ls */*_exons_tandy6jm.gff.gz` set out=genepairs-proteinhsp.txt; touch $out foreach gf ($tdgff) echo "# $gf" >> $out gzcat $gf | perl $td/genepairs.perl -mark proteinhsp >> $out echo '#............................................' >> $out echo >> $out end =item R barplot, see below =cut use Getopt::Long; our $NEARDIST = 15000; # was 15k my $SKIPOVERLAP=0; # ?? want this to weed out tandy overlapping clusters? my $pctgene= 0; my $groupok= 1; my $align= 33; my $skipgroup=""; my $keepgroup=""; my $BYGENE= 0; my $DUMPIDS=0; my $marked=""; my $notmarked=""; my $nc= $ENV{nc} || 2; my $optok= GetOptions( "ngene=i", \$nc, # this is obsolete for genepairs ?? "align=i", \$align, "genecount!", \$BYGENE, #was# $pctgene, "dumpids!", \$DUMPIDS, "groupok!", \$groupok, "skipoverlap!", \$SKIPOVERLAP, "skipgroup=s", \$skipgroup, "keepgroup=s", \$keepgroup, "marked=s", \$marked, # proteinhsp=1 ; estmatch=1; ... filter by some mark "notmarked=s", \$notmarked, # invert marked : proteinhsp=1 ; estmatch=1; ... filter by some mark ); die " usage : cat scaffolds/*tandy.gff | perl genepairs.perl options: -align=$align -[no]genecount -[no]skipoverlap -keepgroup=XXX -skipgroup=DGIL -marked=proteinhsp|estmatch|... -notmarked=invert-marked \n" unless($optok); my (%gs,%in,%ggroup,%ggenes,%gexons,%cxnums,%skipov,%skiprgn); my $criteria=""; $marked= $notmarked if($notmarked); while(<>){ next unless(/^\w/); chomp; my($ref,$src,$typ,$tb,$te,@gff)= split "\t"; if($keepgroup && $typ =~ /region/) { my($ids) = m/ids=([^;]+)/; my($cid) = m/cid=(\d+)/; $skiprgn{"$ref.$cid"}++ unless($ids =~ m/$keepgroup/); } next unless($typ =~ /^(HSP|match)/); my($eqclass)= ( m/tclass=(\w+)/ ) ? $1 : ""; if( m/info=([^;]+)/) { # == $typ eq "match" my $in=$1; my @in=split",",$in; %in= map { split ":" }@in; } my($terepeat)= ( m/terepeat=(\w+)/ ) ? $1 : 0; next if ($in{lowquality}>0 or $terepeat or $src =~ /tandyweak/); $criteria .= " not(lowquality or TErepeat) and" unless($criteria =~ /lowqual/); next unless($in{fullgenes} > 0 or $in{pctalign}>$align ); # $criteria .= " (pctalign>$align or fullgene) and" unless($criteria =~ /pctalign/); if($marked) { my($ismarked)= ( m/$marked=(\w+)/ ) ? $1 : 0; my $not= ($notmarked) ? "not ": ""; $criteria .= " $not$marked and" unless($criteria =~ /$marked/); if($notmarked) { next if($ismarked); } else { next unless($ismarked); } } my $ishsp= ($typ eq "HSP"); my($sid)= $ref; #m/^(\w+)/; $sid =~ s/scaffold_/s/; my($cid)= m/cid=(\d+)/; my($gid)= m/gid=(\d+)/; my($xid)= m/xid=([^;]+)/; # ishsp only my $scid= "$sid.c$cid"; my $scgid= "$scid.g$gid"; next if ($skiprgn{"$ref.$cid"}); if($SKIPOVERLAP) { my $bin= 10000 * int(($tb+$te)/20000); if (!$ishsp && $overlapped{$ref.$bin}++) { $criteria .= " not overlapped and" unless($criteria =~ /overlapped/); my @oldloc= @ { $saveloc{$ref.$bin} || [] }; my $ov= 0; foreach my $rloc (@oldloc) { next if($scid eq $rloc->[3]); #?? same cluster ok $ov=1 and last if _isoverlap( $rloc->[1], $rloc->[2], $tb, $te); } if($ov) { $skipov{$scgid}++; $skipov++; } # mark this gene and all hsps as skipped } next if($skipov>0 and $skipov{$scgid}); push( @{$saveloc{$ref.$bin}}, [$ref,$tb,$te,$scid,$scgid]); } if($ishsp && $xid) { my $cxnum= ++$cxnums{$scid}; # count all hsp/cluster # want to know if ggroup has both eq = 1, -2 at each hsp match, per cid cluster my @xid= split",",$xid; foreach my $xid (@xid) { my ($geneid,$exnum,$xloc,$exeq)= split_exonid($xid,1); my ($ggroup)= $geneid =~ m/^(\D+)/; $ggroup =~ s/DP_DGIL_SNO/SNAP/; # hack fixup; is this only place we collect group from geneid? $ggroup="any" unless($groupok); # for Apis/amel4,1 model set next if ($skipgroup && $ggroup =~ m/$skipgroup/); # above# next if ($keepgroup && $ggroup !~ m/$keepgroup/); ## if missing $exeq test from $xloc _isoverlap/_isnear $tb,$te if( $exeq == -1 ) { if($xloc) { my($xb,$xe)= split(/[-\.]+/,$xloc); if($xe>0) { if(_isoverlap($xb,$xe, $tb,$te)) { $exeq=1; } elsif(_isnear($xb,$xe, $tb,$te)) { $exeq=-2; } else { $exeq= 0; } # far } } } next if( $exeq == -1 ); #? also count uniq geneids / ggroup ? $ggroup{$ggroup}{"genes"}{$geneid}++; # instead of $scgid use ggroup geneid if($pctgene) { $gxnums{$ggroup}{$scid}{$geneid}++; # is this wrong?? use $scgid? my $gxnum= scalar(keys %{$gxnums{$ggroup}{$scid}} ); $cxnum= $gxnum; } $ggroup{$ggroup}{$scid}{$cxnum}{$exeq}++; $ggroup{"all"}{$scid}{$cxnum}{$exeq}++; #?? # count gene ids; BYGENE; But dont want ++ per HSP, per match? $scgid #?? want also exnum count/geneid same/near/far -- for id dumps $ggenes{$ggroup}{$scid}{$geneid}{$exeq}++; $ggenes{"all"}{$scid}{$geneid}{$exeq}++; #?? $gexons{$ggroup}{$scid}{$geneid}{$exeq}{$exnum}++ if($exnum>0); # want link from geneid to all cxnum dupl hsps $gexons{$ggroup}{$scid}{$geneid}{hsps}{$cxnum}{$exeq}++; } } # if(0) { # # check, filter by nexons, pctalign, fullgenes ? # $gs{ngene}{$scid}++ unless($ishsp); # $gs{geneid}{$scid}{$scgid}++ unless($ishsp); # # $gs{nexon}{$scid}++ if($ishsp); # if($eqclass && $ishsp) { $gs{eqnear}{$scid}++ if($eqclass =~ /^near$/); } # # if($eqclass && $ishsp) { $gs{eqfar}{$scid}++ if($eqclass =~ /^far$/); } # # if($eqclass && $ishsp) { $gs{eqsame}{$scid}++ if($eqclass =~ /^equal$/); } # # # $gs{duponly }{$scid}++ #?? fixme / drop? # # if(!$ishsp && $in{eqsame} == 0 && ($in{eqnear} > 0 || $in{eqfar } > 0)); # } #......................... } $criteria=~s/ and$//; show_pairedgroups(); # show_geneclust(); # exit; #---------------------------------- sub show_pairedgroups { my $flag=""; $flag.= ($BYGENE) ? "gene counts" : "exon counts"; print "\nTandem gene clusters: paired same-near exons; $flag\n"; print "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); if($pctgene) { $col1[0]= "Exons"; $col2[0]="/gene"; } ##if($BYGENE) { $col1[0]= "N"; $col2[0]="genes"; } #print join("\t",qw(__________ _N___),@col1),"\n"; #print join("\t",qw(Group_____ clust),@col2),"\n"; $w=15; map{ printf "%".$w."s\t",$_; $w=6; } qw(Gene N N),@col1; print"\n"; $w=15; map{ printf "%".$w."s\t",$_; $w=6; } qw(Group clust genes),@col2; print"\n"; my @ggroup= grep{$_ ne "all"} sort keys %ggroup; foreach my $ggroup (@ggroup) { my %cscore=(); my $ngene= scalar( keys %{ $ggroup{$ggroup}{"genes"} } ); # $cscore{"genes"} = $ngene; my @clust= grep{ $_ ne "genes" } sort keys %{ $ggroup{$ggroup} }; my $nclus= @clust; if($BYGENE) { foreach my $scid (@clust) { my $nx= 0; my %score=(); my @geneids= sort keys %{ $ggenes{$ggroup}{$scid} }; foreach my $geneid (@geneids) { my $eqsame= $ggenes{$ggroup}{$scid}{$geneid}{1} || 0; my $eqnear= $ggenes{$ggroup}{$scid}{$geneid}{-2} || 0; my $eqfar = $ggenes{$ggroup}{$scid}{$geneid}{0} || 0; my $nhsp = $eqsame || 1; ## $ggenes{$ggroup}{$geneid}{"hsp"} || 1; # only same $nx += $nhsp; my $sc= int( _min($eqsame,$eqnear) / $nhsp); $score{"samenear"} += $sc if($sc>0); my $sc= int( ($eqsame - $eqnear) / $nhsp); $score{"same"} += $sc if($sc>0); my $sc= int( ($eqnear - $eqsame) / $nhsp); $score{"near"} += $sc if($sc>0); my $sc= int( ($eqfar) / $nhsp); $score{"far"} += $sc if($sc>0); # my $val=""; # $val.= "same" if($eqsame); # $val.= "near" if($eqnear); # $val.= "far" if($eqfar && !$val); # $score{$val}++; # # fixme for BYGENE # if($val eq "near") { # $eqsame= $ggenes{"all"}{$geneid}{1} || 0; # $score{"nearother"}++ and $score{$val}-- if($eqsame); # } } $cscore{"exons"} += $nx; #? add in hsp counts foreach my $sk (keys %score) { $cscore{$sk} += $score{$sk} ; } } } else { foreach my $scid (@clust) { my @cxnum= sort keys %{ $ggroup{$ggroup}{$scid} }; my $nx= @cxnum; my %score=(); foreach my $cxnum (@cxnum) { ## my @eq= sort keys %{ $ggroup{$ggroup}{$scid}{$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 if($val eq "near") { $eqsame= $ggroup{"all"}{$scid}{$cxnum}{1} || 0; $score{"nearother"}++ and $score{$val}-- if($eqsame); # && $score{$val}--; } } $cscore{"exons"} += $nx; foreach my $sk (keys %score) { ## $cscore{$sk} += $score{$sk} / $nx; #?? or not $cscore{$sk} += $score{$sk} ; } } } $ggroup =~ s/_$//; $ggroup=substr($ggroup,0,14); printf "%15s\t%6d\t%6d\t",$ggroup,$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"; } print "# Overlapping regions skipped: $skipov\n" if $skipov; if($DUMPIDS) { print "# GeneID near exon matches\n"; print "# ",join("\t",qw(Group GeneID nsame nnear exon_same exon_near near_only)),"\n"; my @ggroup= grep{$_ ne "all"} sort keys %ggenes; @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 @clust= grep{ $_ ne "genes" } sort keys %{ $ggroup{$ggroup} }; my @clust= sort keys %{ $ggenes{$ggroup} }; my %dumpid; foreach my $scid (@clust) { my @geneids= sort keys %{ $ggenes{$ggroup}{$scid} }; foreach my $geneid (@geneids) { my $eqsame= $ggenes{$ggroup}{$scid}{$geneid}{1} || 0; my $eqnear= $ggenes{$ggroup}{$scid}{$geneid}{-2} || 0; # my $eqfar = $ggenes{$ggroup}{$scid}{$geneid}{0} || 0; $dumpid{$geneid}{-2} += $eqnear; $dumpid{$geneid}{1} += $eqsame; ($eqsame,$eqnear)=(0,0); $gexons{$ggroup}{$scid}{$geneid}{1} and $eqsame= scalar(keys %{$gexons{$ggroup}{$scid}{$geneid}{1}}); $gexons{$ggroup}{$scid}{$geneid}{-2} and $eqnear= scalar(keys %{$gexons{$ggroup}{$scid}{$geneid}{-2}}); $dumpid{$geneid}{xnear} += $eqnear; $dumpid{$geneid}{xsame} += $eqsame; # add ids with near-only hsps: my @chsps= keys %{$gexons{$ggroup}{$scid}{$geneid}{hsps}}; foreach my $cxnum (@chsps) { my $cnear= $gexons{$ggroup}{$scid}{$geneid}{hsps}{$cxnum}{-2} || 0; my $csame= $gexons{$ggroup}{$scid}{$geneid}{hsps}{$cxnum}{1} || 0; $dumpid{$geneid}{xnearonly} += $cnear if($cnear>0 && $csame==0); } } } ## (my $glab=$ggroup) =~ s/_$//; $glab = substr($glab,0,14); print "# group $ggroup\n"; foreach my $geneid (sort keys %dumpid) { my $cnear= $dumpid{$geneid}{-2} || 0; my $csame= $dumpid{$geneid}{1} || 0; my $xnear= $dumpid{$geneid}{xnear} || 0; my $xsame= $dumpid{$geneid}{xsame} || 0; my $xnonly= $dumpid{$geneid}{xnearonly} || 0; if($cnear) { ## ||$csame print "#g ",join("\t",$geneid,$csame,$cnear,$xsame,$xnear,$xnonly),"\n"; } } print "# ------------------------------------\n\n"; } } } sub show_geneclust { my (%ns, %ts, %ngn); @gn= sort keys %{$gs{ngene}}; @g3= grep { $gs{ngene}{$_} >= $nc } @gn; my $ngnt=0; foreach $gn (@g3) { $nonly{eqnear}++ if ($gs{eqnear}{$gn} > 0 && ! $gs{eqfar }{$gn}); $nonly{eqfar }++ if ($gs{eqfar }{$gn} > 0 && ! $gs{eqnear}{$gn}); $ngnt += $gs{ngene}{$gn}; } print "# criteria: $criteria\n"; print "# gene cluster (n>=$nc) counts\n"; print "total regions: ",scalar(@g3),"\t"; foreach $eq (qw(eqnear eqfar eqsame )) { foreach $gn (@g3) { $ns{$eq}++ if ($gs{$eq}{$gn} > 0); $ts{$eq}+= $gs{$eq}{$gn}; $ngn= $gs{ngene}{$gn}; # = $gs{ngeneid}{$gn}; $ngn{$eq}+= $ngn if ($gs{$eq}{$gn} > 0); #??? } print "$eq: ",$ns{$eq},"\t"; } print"\n"; print "genes/region: $ngnt\t"; foreach $eq (qw(eqnear eqfar eqsame )) { $ns= $ngn{$eq} || 0; print "$eq: ",$ns,"\t"; } print"\n"; print "only class regions: \t"; foreach $eq (qw(eqnear eqfar)) { $ns= $nonly{$eq}; print "$eq: ",$ns,"\t"; } print"\n"; print "average HSPs/region: \t"; foreach $eq (qw(eqnear eqfar eqsame )) { $ns= $ns{$eq} || 1; $ts= sprintf("%.1f",$ts{$eq} / $ns); print "$eq: ",$ts,"\t"; } print"\n"; print "# Overlapping regions skipped: $skipov\n"; } sub _isoverlap { my($gb,$ge, $qb,$qe)= @_; return ($gb <= $qe && $ge >= $qb) ? 1 : 0; } sub _min { my($gb,$ge)= @_; return ($ge >= $gb) ? $gb : $ge; } sub _isnear { # my($gb,$ge, $qb,$qe)= @_; my $gm= int(($gb+$ge)/2); my $qm= int(($qb+$qe)/2); return (abs($gm - $qm) < $NEARDIST) ? 1 : 0; } sub split_exonid { local $_= shift; my $dropnum= shift; $dropnum ||= 0; my ($exnum, $exloc, $exeq)=(0,"",-1); s/\=([\d\-]+)$// and $exeq=$1; s/\:([\d\-]+)$// and $exloc=$1; if($dropnum && s/\.(\d+)$//) { $exnum= $1; }; return (wantarray) ? ($_,$exnum,$exloc,$exeq) : $_; } __END__ ("dsec" "dsim" "dyak" "dere" "dana" "dpse" "dper" "dwil" "dmoj" "dvir" "dgri") set spp=(dpse dyak dana dere dmoj dsec dgri dwil) set spa=(GA GE GF GG GI GM GH GK) foreach sp ( $spp ) echo -n $sp : grep -c "#g.${sp}_GLEANR" genepairs-dmelhspid.txt end foreach sp ( $spa ) echo -n $sp : grep -c "#g.${sp}_NCBI" genepairs-dmelhspid.txt end set tdgff=`ls */*_exons_tandy6jmd.gff.gz` set out=genepairs-dmelhspid.txt; touch $out foreach gf ($tdgff) echo "# $gf" >> $out gzcat $gf | perl $td/genepairs.perl -dumpids -mark dmelhsp >> $out end cat genepairs-dmelhspid.txt | perl -ne\ '$sp=$1 if (m,^# (\w+)/,); $sp=$1 if(/(\w+)_GLEANR/); next unless(/^#g/); \ next if(/RGUI_GID|EISE_CEX|DGIL_SNO|BATZ_CNA|GLEAN_|TRd\w+_/);\ ($c,$id,$n,$s)=split; ($g)= $id=~m/^(\D+)/; next unless($s>=0);\ $gn{$sp}{$g}++ if($n); $gs{$sp}{$g}++ if($s); \ $cn{$sp}{$g} += $n; $cs{$sp}{$g} += $s; \ END{ print join("\t",qw(species group gnear gsame xnear xsame)),"\n"; \ foreach $sp (sort keys %gn) { \ foreach $g (sort keys %{$gn{$sp}}) { \ $n=$gn{$sp}{$g}; $s=$gs{$sp}{$g}; \ $cn=$cn{$sp}{$g}; $cs=$cs{$sp}{$g}; \ print join("\t",$sp,$g,$n,$s,$cn,$cs),"\n"; }}}' species group gnear gsame xnear xsame dsec GM_BREN_NSC_ 1395 1256 4631 6182 dsec GM_EISE_CGW_ 3965 3673 9200 11722 dsec GM_NCBI_GNO_ 2158 2038 7953 11140 dsec GM_PACH_GMP_ 1505 1419 4998 8340 dsec dsec_GLEANR_ 2842 2710 9857 14438 dyak GE_BREN_NSC_ 1565 1427 5751 6820 dyak GE_EISE_CGW_ 2548 2362 7390 8712 dyak GE_NCBI_GNO_ 2477 2353 9561 12312 dyak GE_PACH_GMP_ 1827 1742 6625 9575 dyak dyak_GLEANR_ 2489 2386 9295 12804 dere GG_BREN_NSC_ 1375 1253 5469 6781 dere GG_EISE_CGW_ 3328 3099 9966 10645 dere GG_NCBI_GNO_ 2248 2143 10825 12479 dere GG_PACH_GMP_ 1739 1669 7094 9999 dere dere_GLEANR_ 2200 2098 8302 12249 dana GF_BREN_NSC_ 1278 1105 5130 5473 dana GF_EISE_CGW_ 2030 1874 6187 8398 dana GF_NCBI_GNO_ 2277 2135 9834 12007 dana GF_PACH_GMP_ 1842 1769 6892 11097 dana dana_GLEANR_ 2356 2222 9444 12596 dpse GA_BREN_NSC_ 1757 1586 7283 7547 dpse GA_EISE_CGW_ 3481 3254 11494 12373 dpse GA_NCBI_GNO_ 2943 2805 14570 16947 dpse GA_PACH_GMP_ 2299 2238 9711 14775 dpse dpse_GLEANR_ 2793 2686 11879 16120 dwil GK_BREN_NSC_ 998 852 3901 4027 dwil GK_EISE_CGW_ 1763 1637 6469 8241 dwil GK_NCBI_GNO_ 2553 2409 13768 15105 dwil GK_PACH_GMP_ 2039 1961 7097 13075 dwil dwil_GLEANR_ 2273 2141 9817 12478 dmoj GI_BREN_NSC_ 1194 1048 4714 5069 dmoj GI_EISE_CGW_ 3025 2787 10233 11037 dmoj GI_NCBI_GNO_ 2733 2618 15093 16945 dmoj GI_PACH_GMP_ 1861 1781 6814 12079 dmoj dmoj_GLEANR_ 2403 2277 9955 14029 dgri GH_BREN_NSC_ 1191 1039 5178 5079 dgri GH_EISE_CGW_ 1793 1700 7708 9047 dgri GH_NCBI_GNO_ 2735 2649 18817 17328 dgri GH_PACH_GMP_ 1682 1642 6724 11084 dgri dgri_GLEANR_ 2277 2185 10908 13830 dmel4d CG_BATZ_CON_ 746 610 3433 2003 dmel4d CG_NCBI_GNO_ 1541 1424 6335 6668 dmel4d FBtr 1399 1336 5900 6503 daphe Dappu 1254 1128 5025 6746 daphe NCBI_GNO_ 1554 1398 6962 8280 R stats: setwd("~/Desktop/dspp-work/genomesoft/tandy/") gp2 <- read.table("genepairs.dat", header=T, strip.white=T) gp4 <- read.table("genepairs-prots.dat", header=T, strip.white=T) title<-"tandem genes, homolog filter" gp <- gp4 # ^^ cant use total genes adjust here as -prots subset differs / species gp5 <- read.table("genepairs-notprots.dat", header=T, strip.white=T) title<-"tandem genes, NOT homolog filter" gp <- gp5 gp6 <- read.table("genepairs-dmelhsp.dat", header=T, strip.white=T) title<-"tandem genes, DroMel homologs" gp <- gp6 gp7 <- read.table("genepairs-dmelnot.dat", header=T, strip.white=T) title<-"tandem genes, NOT DroMel homologs" gp <- gp7 lev1 <- c("caeele_Coding_tr",grep("NCBI",as.character(gp$Group),value=T)) gfac <- factor(gp$Group, levels=lev1, ordered=T,exclude=NULL) gpa <- gp[!is.na(gfac),-1] rownames(gpa) <- gsub("Coding_tr","WB",gsub("NCBI_GNO","NCBI",gp[!is.na(gfac),1])) ex2gn <- 1 ; ylabel <- "duplicate exons" # ex2gn <- 1/gpa[,"tgenes"]; ylabel<- "duplicates / total genes" # ex2gn <- gpa[,"ngenes"]/gpa[,"nexons"] # no good; doesnt transform exon counts to gene c layout(as.matrix(1)) par(mar=c(4,4,3,2)+0.1) barplot(as.matrix(t( ex2gn * gpa[,c(4,6,7,8)])),ylab=ylabel, main=title, beside=T, legend=T, col = c("lightblue", "mistyrose", "lightcyan", "lavender") ) # xset <- c(5,4,6,7,8); xcol<-c("grey90","lightblue", "mistyrose", "lightcyan", "lavender") # xset <- c(4,6,7,8); xcol<-c("lightblue", "mistyrose", "lightcyan", "lavender") # .. predictor comparison # drop dromel, dappux, caele from pred.compare; add drosim, drovir, droper sppset <- c("drosec","drosim","droyak","droere","droana","dropse","droper", "drowil","dromoj","drovir","drogri") spat <- paste(sppset,collapse="|") predset <- c("NCBI_GNO","GLEANR","EISE_CGW","PACH_GMP") # BREN_NSC ? PACH_GMP? others? np <- length(predset) #> par("mar") == 5.1 4.1 4.1 2.1 == bot,lef,top,ri layout(nf <- matrix(1:np, np, 1, byrow=TRUE) ) layout.show(n=np) par(cex=1) par(mar=c(2,4,3,2)+0.1) ymax <- 3000 # 8000 # 2000 ; max(gp[,5]) for (pd in predset) { levs <- grep(spat,grep(pd,as.character(gp$Group),value=T),value=T) gfac <- factor(gp$Group, levels=levs, ordered=T,exclude=NULL) gpa <- gp[!is.na(gfac),-1] rown <- gp[!is.na(gfac),1] # for only 1 predictor/x-axis, put predname on xlab, species in rownames rown <- gsub("_.*","",rown) # drop predictor rown <- gsub("dro","D",rown) rownames(gpa) <- rown xlabel <- pd xlabel <- gsub("NCBI_GNO","Gnomon",xlabel) xlabel <- gsub("EISE_CGW","GeneWise",xlabel) xlabel <- gsub("PACH_GMP","GeneMapper",xlabel) xlabel <- gsub("GLEANR","GleanR",xlabel) predtitle <- paste(title,", Predictor: ",xlabel) ; # was Precictor :( ex2gn <- 1 ; ylabel <- "duplicate exons" #tparm <- c(4,6,7,8) tparm <- c(4,6,7) tcols <- c("lightblue", "mistyrose", "lightcyan", "lavender") [1:length(tparm)] bpx <- barplot( as.matrix(t( ex2gn * gpa[,tparm])),ylab=ylabel, main=predtitle, ylim=c(0,ymax), beside=T, legend=F, col = tcols ) # add line from dsec/same_near top left to right abline(h=ex2gn * max(gpa["Dsec",tparm]), lty="dotted") } # one legend at end: legend( x="topleft", colnames(gpa[,tparm]), lwd=3, col=tcols) #......... levgr <- c(grep("GLEANR",as.character(gp$Group),value=T)) # caeele_Coding_tr, dappul_JGI, dromel_FBtr grfac <- factor(gp$Group, levels=levgr, ordered=T,exclude=NULL) gra <- gp[!is.na(grfac),-1] rownames(gra) <- gsub("Coding_tr","WB",gsub("NCBI_GNO","NCBI",gp[!is.na(grfac),1])) barplot(as.matrix(t(gra[,c(5,4,6,7,8)]/gra[,"tgenes"])),ylab="tandems / total genes", beside=T,legend=T,col = c("grey90","lightblue", "mistyrose", "lightcyan", "lavender") ) #............ # droyak_EISE_CGW levgw <- c(grep("CGW",as.character(gp$Group),value=T)) gwfac <- factor(gp$Group, levels=levgw, ordered=T,exclude=NULL) gwa <- gp[!is.na(gwfac),-1] rownames(gwa) <- gsub("Coding_tr","WB",gsub("EISE_CGW","CGW",gp[!is.na(gwfac),1])) bwa <- barplot(as.matrix(t(gwa[,c(5,4,6,7,8)]/gwa[,"tgenes"])),ylab="tandems / total genes", beside=T,legend=T,col = c("grey90","lightblue", "mistyrose", "lightcyan", "lavender") ) #...... gpa <- gp[!is.na(gfac),-1] rownames(gpa) <- gsub("Coding_tr","WB",gsub("NCBI_GNO","NCBI",gp[!is.na(gfac),1])) gpa nclust ngenes nexons same_near same near_other near_only far_only tgenes texons caeele_WB 4451 5992 36694 6687 20900 1652 3972 4177 21616 123981 dappul_NCBI 6861 8652 60597 12437 29998 2616 7039 9644 27084 137475 dromel_NCBI 1903 2850 13217 3026 8633 694 943 213 15723 56422 barplot(as.matrix(t(gpa[,c(5,4,6,7,8)]/gpa[,"tgenes"])),ylab="tandems / total genes", beside=T,legend=T,col = c("grey90","lightblue", "mistyrose", "lightcyan", "lavender") ) .... gpn <- gp[!is.na(gfac),c(6,5,7,8)] rownames(gpn) <- gsub("Coding_tr","WB",gsub("NCBI_GNO","NCBI", gp[!is.na(gfac),"Group"])) > gpn same_near same near_other near_only caeele_WB 6687 20900 1652 3972 dappul_NCBI 12437 29998 2616 7039 dromel_NCBI 3026 8633 694 943 droyak_NCBI 5598 15735 1368 2223 droere_NCBI 5807 17469 1354 2484 dropse_NCBI 8721 23236 1539 2485 dromoj_NCBI 9410 25166 1851 3443 drogri_NCBI 6725 15083 1239 2262 drowil_NCBI 9327 27060 2241 4398 > barplot(as.matrix(t(gpn))) > barplot(as.matrix(t(gpn)),beside=T) barplot(as.matrix(t(gpn)),beside=T,legend=T, col = c("lightblue", "mistyrose", "lightcyan", "lavender") ) > gpp <- gpn[,2:4]/gpn[,1] barplot(as.matrix(t(gpp)),beside=T,legend=T, col = c( "mistyrose", "lightcyan", "lavender") ) > gp[1:10,] Group nclust ngenes nexons same_near same near_other near_only far_only tgenes 1 dappul_DGIL_SNO 1944 2814 12179 2219 8220 487 767 726 37172 2 dappul_JGI 1849 2298 12869 2315 8025 1083 825 962 17350 3 dappul_NCBI_GNO 1996 2779 14523 3393 9151 629 911 715 27084 4 dromel_BATZ_CON 1650 2314 9829 1143 7422 574 726 201 14132 5 dromel_DGIL_SNO 1816 2799 10202 2098 6793 582 793 180 22244