File Coverage

blib/lib/obogaf/parser.pm

Criterion	Covered	Total	%
statement	192	192	100.0
branch	83	86	96.5
condition	3	3	100.0
subroutine	11	11	100.0
pod	6	6	100.0
total	295	298	98.9

package obogaf::parser;

require 5.006;
our $VERSION= '1.271';
$VERSION= eval $VERSION;

use strict;
use warnings;
use Graph;
use IO::File;
use PerlIO::gzip;

sub build_edges{
    my ($obofile)= @_;
    my ($namespace, $idname, $isname, $pofname, $source, $destination, $pof, $res);
    if($obofile=~/.obo$/){ open FH, "<", "$obofile" or die "cannot open $obofile. $!.\n"; } else { die "cannot open $obofile. The extension must be obo.\n"; }
    while(){
        chomp;
        next if $_=~/^\s*$/;
        if($_=~/^namespace:\s+(\D+)/){
            $namespace=$1;
        }elsif($_=~/^name:\s+(.+)/){
            $idname=$1;
        }elsif($_=~/^id:\s+(\D+\d+)/){
            $destination=$1;
        }elsif($_=~/^is_a:\s+(\D+\d+)/){
            $source=$1;
            ($isname)= ($_=~/!\s+(.+)/);
            if(defined $namespace){
                $res .= "$namespace\t$source\t$destination\t$isname\t$idname\tis-a\n";
            }else{
                $res .= "$source\t$destination\t$isname\t$idname\tis-a\n";
            }
        }elsif($_=~/^relationship: part_of\s+(\D+\d+)/){
            $pof=$1;
            ($pofname)= ($_=~/!\s+(.+)/);
            if(defined $namespace){
                $res .= "$namespace\t$pof\t$destination\t$pofname\t$idname\tpart-of\n";
            }else{
                $res .= "$pof\t$destination\t$pofname\t$idname\tpart-of\n";
            }
        }
    }
    close FH;
    return \$res;
}

sub build_subonto{
    my ($edgesfile, $namespace)= @_;
    my ($res, %checker);
    open FH, "<", $edgesfile or die "cannot open $edgesfile. $!.\n";
    while(){
        next if $_=~/^[!,#]|^\s*$/;
        my @vals= split(/\t/, $_);
        $checker{$vals[0]}=1;
        if($vals[0] eq $namespace){ $res .= join("\t", @vals[1..$#vals]); }
    }
    close FH;
    unless(exists($checker{$namespace})){die "$edgesfile does not include $namespace or $namespace is not in the first column of $edgesfile.\n";}
    return \$res;
}

sub make_stat{
    my ($edgesfile, $parentIndex, $childIndex)= @_;
    my (%indeg, %outdeg, %deg, $ed, $nd, $mindeg, $maxdeg, $medeg, $avgdeg, $den, $scc, $resdeg, $stat, $res);
    ## create graph
    my $g= Graph->new(directed => 1);
    open FH, "<", $edgesfile or die "cannot open $edgesfile. $!.\n";
    while(){
        chomp;
        my @vals= split(/\t/,$_);
        $g->add_edge($vals[$parentIndex], $vals[$childIndex]);
    }
    close FH;
    ## compute indegree/outdegree/degree
    my @V= $g->vertices;
    foreach my $nd (@V){
        my $i= $g->in_degree($nd);
        my $o= $g->out_degree($nd);
        my $d= $i+$o;
        $indeg{$nd}=$i;
        $outdeg{$nd}=$o;
        $deg{$nd}=$d;
    }
    foreach my $node (sort{$deg{$b}<=>$deg{$a} or ($a cmp $b)} keys %deg){ $resdeg .= "$node\t$deg{$node}\t$indeg{$node}\t$outdeg{$node}\n"; }
    ## compute: median/max/min degree
    my @sortdeg= sort{$a<=>$b} values (%deg);
    my $len= $#sortdeg+1;
    my $mid = int $len/2;
    if($len % 2){ $medeg = $sortdeg[$mid]; }else{ $medeg = ( $sortdeg[$mid-1] + $sortdeg[$mid] ) / 2; }
    $medeg= sprintf("%.4f", $medeg);
    $mindeg= $sortdeg[0];
    $maxdeg= $sortdeg[$#sortdeg];
    ## compute number of nodes and edges
    $ed= $g->edges;
    $nd= $g->vertices;
    ## compute average degree and density
    $avgdeg= $ed/$nd;
    $den= $ed / ( $nd * ($nd -1) );
    $avgdeg= sprintf("%.4f", $avgdeg);
    $den= sprintf("%.4e", $den);
    ## return stat
    $stat .= "nodes: $nd\nedges: $ed\nmax degree: $maxdeg\nmin degree: $mindeg\n";
    $stat .= "median degree: $medeg\naverage degree: $avgdeg\ndensity: $den\n";
    $res= "#oboterm  degree  indegree  outdegree\n".$resdeg."\n"."~summary stat~\n".$stat;
    return $res;
}

sub get_parents_or_children_list{
    my ($edgesfile, $parentIndex, $childIndex, $chdORpar)= @_;
    my (%nodelist);
    if($chdORpar ne "parents" && $chdORpar ne "children"){ die "$chdORpar can be 'parents' or 'children'.\n";}
    open FH, "<", $edgesfile or die "cannot open $edgesfile. $!.\n";
    while(){
        chomp;
        my @vals= split(/\t/,$_);
        if($chdORpar eq "parents"){
            $nodelist{$vals[$childIndex]} .= $vals[$parentIndex]."|";
        }else{
            $nodelist{$vals[$parentIndex]} .= $vals[$childIndex]."|";
        }
    }
    close FH;
    foreach my $term (keys %nodelist){ chop $nodelist{$term}; }
    return \%nodelist;
}

sub gene2biofun{
    my ($annfile, $geneIndex, $classIndex)= @_;
    my (%gene2biofun, @genes, @biofun, $stat)= ();
    my ($sample, $oboterm)= (0)x2;
    if ($annfile =~ /.gz$/){
        open FH, "<:gzip", $annfile or die "cannot open $annfile. $!.\n";
    }else{
        open FH, "<", "$annfile" or die "cannot open $annfile. $!.\n";
    }
    while(){
        next if $_=~/^[!,#]|^\s*$/;
        chomp;
        my @vals=split(/\t/,$_);
        push(@genes, $vals[$geneIndex]);
        push(@biofun, $vals[$classIndex]);
        $gene2biofun{$vals[$geneIndex]} .= $vals[$classIndex]."|";
    }
    close FH;
    foreach my $gene (keys %gene2biofun){ chop $gene2biofun{$gene}; }
    my %seen=();
    my @uniqgenes= grep{!$seen{$_}++} @genes;
    $sample= scalar(@uniqgenes);
    undef %seen;
    my @uniqpbiofun= grep{!$seen{$_}++} @biofun;
    $oboterm= scalar(@uniqpbiofun);
    $stat .= "genes: $sample\nontology terms: $oboterm\n";
    return (\%gene2biofun, \$stat);
}

sub map_OBOterm_between_release{
    my ($obofile, $annfile, $classIndex)= @_;
    my (%altid, %oldclass, %old2new, $header, $id, $fln, $pair, $stat, $pstat);
    my ($alt, $classes, $seen, $unseen)= (0)x4;
    ## step 0: pairing altid_2_id (key: alt_id)
    if($obofile=~/.obo$/){ open FH, "<", "$obofile" or die "cannot open $obofile. $!.\n"; } else { die "cannot open $obofile. The extension must be obo.\n"; }
    while (){
        chomp;
        next if $_=~/^\s*$/;
        if($_=~/^id:\s+(\D+\d+)/){ $id=$1; }
        if($_=~/^alt_id:\s+(\D+\d+)/){ $altid{$1}=$id; }
    }
    close FH;
    $alt= keys(%altid);
    # step 1: storing old ontology terms in a hash
    if ($annfile =~ /.gz$/){
        open FH, "<:gzip", $annfile or die "cannot open $annfile. $!.\n";
    }else{
        open FH, "<", "$annfile" or die "cannot open $annfile. $!.\n";
    }
    while(){
        chomp;
        if($_=~/^[!,#]|^\s*$/){ $header .= "$_\n"; }
        next if $_=~/^[!,#]|^\s*$/;
        my @vals=split(/\t/,$_);
        $oldclass{$vals[$classIndex]}=$vals[$classIndex];
    }
    close FH;
    $classes= keys(%oldclass);
    ## step 2: mapping old GO terms to the new one using *alt_id* as key
    my $tmp= "";
    foreach my $k (sort{$a cmp $b} keys(%altid)){
        if($oldclass{$k}){
            $old2new{$k}=$altid{$oldclass{$k}};  ## pairing
            $seen++;
            $tmp= $altid{$oldclass{$k}};
        }else{
            $tmp= "unseen";
            $unseen++;
        }
        if($tmp ne "unseen"){
            $pair .= "$k\t$altid{$oldclass{$k}}\n";
        }
    }
    ## step 3: substitute ALT-ID with the updated ID, then the annotation file is returned.
    if ($annfile =~ /.gz$/){
        open FH, "<:gzip", $annfile or die "cannot open $annfile. $!.\n";
    }else{
        open FH, "<", "$annfile" or die "cannot open $annfile. $!.\n";
    }
    while(){
        chomp;
        next if $_=~/^[!,#]|^\s*$/;
        my @vals= split(/\t/, $_);
        my $oboterm= $vals[$classIndex];
        if($old2new{$oboterm}){
            $oboterm= $old2new{$oboterm};
            $_=~ s/$vals[$classIndex]/$oboterm/g;
            $fln .= "$_\n";
        }else{
            $fln .= "$_\n";
        }
    }
    close FH;
    if(defined $header){$fln = $header.$fln;}
    ## print mapping stat
    $stat .= "Tot. ontology terms:\t$classes\nTot. altID:\t$alt\nTot. altID seen:\t$seen\nTot. altID unseen:\t$unseen\n";
    unless(not defined $pair){
        $pstat .= "#alt-id  id\n$pair\n$stat";
        return (\$fln, \$pstat);
    }
    return (\$fln, \$stat);
}

1;

__END__