File Coverage

blib/lib/obogaf/parser.pm

Criterion	Covered	Total	%
statement	228	228	100.0
branch	96	102	94.1
condition	5	6	83.3
subroutine	13	13	100.0
pod	7	7	100.0
total	349	356	98.0

package obogaf::parser;
require 5.006;
our $VERSION= '1.373';
$VERSION= eval $VERSION;

require Exporter;
our @ISA= qw(Exporter);
our %EXPORT_TAGS=(
    all => [qw(
        &build_edges
        &build_subonto
        &make_stat
        &get_parents_or_children_list
        &obo_filter
        &gene2biofun
        &map_OBOterm_between_release
    )],
);
our @EXPORT_OK= (@{$EXPORT_TAGS{'all'}});

use strict;
use warnings;
use Graph;
use List::MoreUtils qw(uniq);
use IO::File;
use PerlIO::gzip;

sub build_edges{
    my ($obofile)= @_;
    my ($namespace, $idname, $isname, $pofname, $source, $destination, $pof, $res);
    if($obofile=~/.obo$/i){ open FH, "<", "$obofile" or die "cannot open $obofile. $!.\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, %outlist);
    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"){
            push (@{$nodelist{$vals[$childIndex]}}, $vals[$parentIndex]);
        }else{
            push (@{$nodelist{$vals[$parentIndex]}}, $vals[$childIndex]);
        }
    }
    close FH;
    foreach my $term (keys %nodelist){
        $outlist{$term} = join('|', sort{($a=~/(\d+)/)[0] <=> ($b=~/(\d+)/)[0]} uniq @{$nodelist{$term}});
    }
    return \%outlist;
}

sub obo_filter{
    my ($obofile, $termsfile)= @_;
    my (@oboterms, $header, $counter, $res, $obo);
    ## obo terms of interest
    open FH, "<", "$termsfile" or die "cannot open $termsfile. $!.\n";
    while(){
        chomp;
        next if $_=~/^\s*$/;
        push(@oboterms, $_);
    }
    close FH;
    my @unique= do { my %seen; grep { !$seen{$_}++ } @oboterms};
    @oboterms=@unique;
    ## store obofile header
    if($obofile=~/.obo$/i){ open FH, "<", "$obofile" or die "cannot open $obofile. $!.\n"; }
    while(){
        next unless 1 .. /\[Term\]/;
        next if /^\[Term\]/;
        $header .= $_;
    }
    ## extract from obo file the terms of interest
    foreach my $oboterm (@oboterms){
        $counter++;
        open FH, $obofile;
        while(){
            if(/^id:\s+($oboterm)/ .. (/^$/ || eof FH)){ ## check for end of file
                $res .= $_;
            }
        }
        close FH;
        print($counter, "\t", $oboterm, "\tdone\n");
    }
    unless(defined $res){ die "none obo terms listed in $termsfile file was found in $obofile file. $!.\n" }
    ## print obofile filtered
    $obo = $header.$res."\n";
    $obo=~ s/\b(id:\s+(\D+\d+))/\[Term\]\n$1/g; ## note: word boundary (\b) is necessary to match only "id:"" and not "alt_id:" in the string $obo
    return \$obo;
}

sub gene2biofun{
    my ($annfile, $geneIndex, $classIndex)= @_;
    my (%gene2biofun, %biofun, @genes, @biofun, $stat)= ();
    my ($sample, $oboterm)= (0)x2;
    if($annfile=~/.gz$/i){ 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]);
        push(@{$biofun{$vals[$geneIndex]}}, $vals[$classIndex]);
    }
    close FH;
    foreach my $gene (keys %biofun){
        $gene2biofun{$gene} = join('|', sort{($a=~/(\d+)/)[0] <=> ($b=~/(\d+)/)[0]} uniq @{$biofun{$gene}});
    }
    my @uniqgenes= uniq @genes;
    my @uniqpbiofun= uniq @biofun;
    $sample= scalar(@uniqgenes);
    $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$/i){ open FH, "<", "$obofile" or die "cannot open $obofile. $!.\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$/i){ 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$/i){ 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__