#!/usr/bin/perl
# swlist2gff.perl
=item notes
# attach to gff attrib, swiss annots (from list table not full sw records) and tandy ids
cat tandyfourprot.pa?.list dpulex1_tandy6gfour.genes dpulex1_ncbigno_tandyfour2.gff \
| perl sw2gff.pl > dpulex1_ncbigno_tandyfour2ann1.gff
# view short notes:
grep Note= dpulex1_ncbigno_tandyfour2ann.gff | perl -ne\
'chomp; m/Note=([^;]+)/;$n=$1; $n=~s,.*/,,; \
$n=~s/,[A-Z0-9]{2,6}_.*$//; \
$n=~s/(\w),? [0-9][A-Z0-9\.,-]+$/$1/; \
$n=~s/(Cytochrome P450).*/\u$1/i; $n=~s/(Histone).*/\u$1/i;\
$n=~s/,\s*isoform \w+//; \
$n=~s/^(ENS\w+\d+|CG\d+[\w-]+),\s*(\w)/$2/; $n=~s/^.*product://; \
$n=~s/,$//; print "$n\n";' | sort | uniq -c | sort -k1,1nr | more
# add here, elsewhere EST-treatment class info
CG31522-PD, isoform D,GH22993p
=cut
use strict;
use Getopt::Long;
# use warnings;
my $debug=0;
my $act= "summary"; # summary, stitch
my $html=1;
my $url_td="http://server2.eugenes.org/cgi-bin/gbrowse/dpulex_jgi060905/";
my $optok= GetOptions(
"action=s", \$act,
"html!", \$html,
"debug!", \$debug,
);
stitch_swiss() if($act =~ /stitch/); # cat tandyfourprot.pa?.list dpulex1_tandy6gfour.genes dpulex1_ncbigno_tandyfour2.gff |
summary_notes() if($act =~ /summary/); # grep Note= dpulex1_ncbigno_tandyfour2ann.gff
# add html table : links to gmap; drop Gene_IDs?
sub summary_notes {
my (%notes,%tandy);
while(<>){
next unless(/Note=/);
s/\s+$//; # my @v=split"\t";
my($id)=m/ID=([^;]+)/;
my($tandy)=m/tandy=([^;]+)/; #? better than ID for output?
my($scaf)=m/^(\S+)/; $scaf =~ s/scaffold_/s/; # really need to add to tandy ids
m/Note=([^;]+)/; my $n=$1;
$n=~s,.*/,,;
$n=~s/,[A-Z0-9]{2,6}_.*$//;
$n=~s/(\w),? [0-9][A-Z0-9\.,-]+$/$1/;
$n=~s/,?\s+isoform\s+\w+//g;
$n=~s/.*(Cytochrome P450).*/\u$1/i;
$n=~s/.*(Hemoglobin).*/\u$1/i;
$n=~s/Predicted protein/Hypothetical protein/i;
$n=~s/Novel protein/Hypothetical protein/i;
$n=~s/^(ENS\w+\d+|CG\d+[\w-]+),\s*(\w)/$2/;
$n=~s/^.*product://;
$n=~s/,$//;
#print "$n\n";
# do here: sort | uniq -c | sort -k1,1nr
$notes{$id}= $n;
my @tandy= map{ $scaf.$_ } split",",$tandy;
push( @{$tandy{$id}},@tandy) if($tandy);
}
my (%count,%nids);
my @ids= sort{ $notes{$a} cmp $notes{$b} } keys %notes;
foreach my $id (@ids) {
my $n= $notes{$id};
$count{$n}++;
push( @{$nids{$n}},$id);
}
my $ntot=0;
if($html) {
print "
Daphnia Sample Duplicate Genes
Daphnia pulex genome: Selected Tandem Duplicate Genes
Tandem duplicate gene clusters, selected from those with homology
to other organisms. Descriptions here suggested from homology
may not be the most accurate.
Duplicate genes were located with 'tandy'
duplicate finding software.
Don Gilbert, June 2007, gilbertd\@indiana.edu
| N_dupl |
Homolog Description |
Tandem IDs |
";
} else {
my $nl= sprintf "%-45s", "Description";
print join("\t","N_dupl",$nl,"Tandem_IDs"),"\n"; #"Gene_IDs",
}
my @cnotes= sort{$count{$b} <=> $count{$a} or $a cmp $b} keys %count;
foreach my $cn (@cnotes) {
my @nids= sort @{ $nids{$cn} };
my $n= scalar(@nids); $ntot += $n;
my @tds= map{ @{$tandy{$_}} } @nids;
my %tds= map{ s/g\d+//; $_,1; } @tds; # drop gene num
@tds= sort{
my ($sa)= $a =~ m/s(\d+)/;
my ($sb)= $b =~ m/s(\d+)/;
$sa<=>$sb or $a cmp $b } keys %tds;
my $etc=""; if(@tds>8) { @tds= splice(@tds,0,8); $etc=" .."; }
my $tds= join(",",@tds).$etc; # shorten long list for view?
# way to many ids to show all, but print here?
my $maxn= ($n>2) ? 2 : $n-1;
my $nido= join(",",@nids[0..$maxn]);
$nido.=" .." if($maxn<$n-1);
if($html) {
# ? http://server2.eugenes.org/cgi-bin/gbrowse/dpulex_jgi060905/?ref=scaffold_1;name=td_c173
my $flag="";
$flag=' bgcolor="pink"' if($cn=~/Hemoglobin/);
my $tds= join ",", map {
my ($sc)= $_ =~ m/s(\d+)/;
my $tdv= $_;
#$tdv =~ s/s\d+//; ## FIXME: tandy GFF needs scaffoldID in tandy id
$tdv= "td_$tdv";
"$_"
##"$_"
} @tds;
print "
$n |
$cn |
$tds $etc |
";
} else {
$cn= sprintf "%-45s", substr($cn,0,45);
print join("\t",$n,$cn,$tds),"\n"; #$nido,
}
}
if($html) {
print "
";
} else {
print "# Total listed dupl. gene set=$ntot\n";
}
}
#Putative core 3-alpha-L-fucosyltransferase
#.........1.........2.........3.........4...
sub stitch_swiss {
my (%td, %ptab, );
while(<>){
s/\s+$//; my @v=split"\t";
# tandy gff: match gene predict ID to tandy ID
if(/^scaffold_\d+\standy/){ # intandy
my ($tdid)= $v[-1] =~ m/ID=td_(\w+)/;
while(/(NCBI_GNO_\d+)/g){ my $id=$1; $td{$id}.=$tdid."," unless($td{$id} && $td{$id}=~m/$tdid,/); }
}
# uniprot id match table
# ID ACC ENTNAME LENGTH ORG TAXGRP NF90 NF50 Matched Query ID(s)
#271 HSP90_CAEEL Q18688 Heat shock protein 90; Abnormal dauer formation protein 21 702 Caenorhabditis elegans Euk/Animal UniRef90_Q18688 UniRef50_P07900 NP_506626,
elsif(/^\d/){ # inptab
my @id= grep/\w/, split(/[,;]/, $v[-1]);
$v[2]=~s/\s*;\s*/,/g; $v[3]=~s/\s*;\s*/,/g;
my $pv=join("/",@v[1,2,5,3]);
map{ $ptab{$_}= $pv; }@id;
}
# main gene predict gff
elsif(/^scaffold/){ # (/^scaffold_\d+\sNCBI_GNO\smRNA/)
my $ptab=""; my $at=$v[-1];
# replace all those Fly CGnnn with no annot for Mouse alternate description
if($at=~/protein_hit=([^;\s]+)/) {
my $ph=$1;
my($gl,$gi,$db,$dx)=split(/[|]/,$ph); $dx=~s/\.\d+$//;
#$ptab .= $ptab{$dx}."," if $ptab{$dx};
$ptab = $ptab{$dx} if $ptab{$dx};
}
unless(0) {
$at =~ m/Dbxref=([^;\s]+)/; my @dx= split",", $1; my %pv=();
foreach (@dx){ s/\.\d+$//;
if(my $pv= $ptab{$_}) {
# alternate ,Homo sapiens, if no Mus ?
if($ptab && $ptab=~m,Drosophila, && $ptab =~ m,/CG\d+, && $pv=~m,Mus musculus,){
$ptab = $pv; } # prefer Mouse over useless Fly descrip=cgid
elsif($ptab){ $pv=~s,/.*,,; $pv{$pv}++; }
else{ $ptab = $pv; } #$ptab .= "$pv,";
}
}
# $ptab.= join",",sort keys %pv if(%pv); # drop this?
}
my($id)=m/(NCBI_GNO_\d+)/; if(my $td=$td{$id}){ $v[-1].=";tandy=$td"; }
$v[-1].=";Note=$ptab" if($ptab);
print join("\t",@v),"\n";
}
}
}