OMA.drw

# ***** BEGIN LICENSE BLOCK *****
# Version: MPL 2.0
#
# The contents of this file are subject to the Mozilla Public License Version
# 2.0 (the "License"); you may not use this file except in compliance with
# the License. You may obtain a copy of the License at
# http://www.mozilla.org/MPL/2.0/
#
# Software distributed under the License is distributed on an "AS IS" basis,
# WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
# for the specific language governing rights and limitations under the
# License.
#
# The Original Code is OMA standalone.
#
# The Initial Developer of the Original Code is CBRG Research Group; 
# ETH Zurich; Switzerland.
# Portions created by the Initial Developer are Copyright (C) 2005-2015
# the Initial Developer. All Rights Reserved.
#
# Contributor(s):
#   Christophe Dessimoz <cdessimoz@ucl.ac.uk>
#   Adrian Altenhoff <adrian.altenhoff@inf.ethz.ch>
#   Stefan Zoller <stefan.zoller@inf.ethz.ch>
#   Adrian Schneider <adrian.schneider@inf.ethz.ch>
#   Alexander Roth <alexander.roth@inf.ethz.ch>
#   Gaston Gonnet <gonnet@inf.ethz.ch>
#
# ***** END LICENSE BLOCK *****
# Set(quiet):
##############################################################################
#       ___           ___           ___
#      /\  \         /\__\         /\  \    
#     /::\  \       /::|  |       /::\  \   
#    /:/\:\  \     /:|:|  |      /:/\:\  \  
#   /:/  \:\  \   /:/|:|__|__   /::\~\:\  \ 
#  /:/__/ \:\__\ /:/ |::::\__\ /:/\:\ \:\__\
#  \:\  \ /:/  / \/__/~~/:/  / \/__\:\/:/  /
#   \:\  /:/  /        /:/  /       \::/  / 
#    \:\/:/  /        /:/  /        /:/  /  
#     \::/  /        /:/  /        /:/  /   
#      \/__/         \/__/         \/__/    
#
#       Stand-alone version
#          (c) CBRG (ETH Zurich) 2004-2015
#          (c) Dessimoz Lab (UCL and UNIL) 2013-2015
##############################################################################
# Usage
#
#   This script identifies orthologs among genomes using the OMA procedure.
#   (http://omabrowser.org).
#
#   If you received this script per standalone installer for OMA, please have
#   also a look at the included README.oma file in the installation folder
#   (default location is /usr/local/oma).
#
#   - Run this script from the same directory.
#   - Script parameters are set in parameters.drw in the same directory
#   - Each genome should have its own fasta file (filename ending with .fa)
#     and be placed in a directory called DB/
#
##############################################################################

Set(gc=1e8):
SetRandSeed():
GetStats := true;

if assigned(AuxDataPath) then 
    datadirname := AuxDataPath;
    if datadirname[1..2] = '~/' then 
        datadirname := getenv('HOME').datadirname[2..-1]
    fi:
else
    datadirname := getenv('HOME').'/.cache/oma';
fi:

# Read libraries
if assigned(libdir) then
    if libdir[-1]<>'/' then libdir := libdir .'/' fi:
    ReadProgram( libdir . 'darwinit'):
else
    ReadProgram('lib/darwinit'):
fi:

# Convert FASTA databases to darwin databases
print('Starting database conversion and checks...');

if not assigned(Datasets) then
    G := TimedCallSystem('find -L DB -name "*.fa" -print')[2];
    G := sort(SearchDelim('\n',G));
else
    G := Datasets;
fi:

# get nr of total parallel process and number of this job.
pinf := traperror(DetectParallelInfo()):
if pinf=lasterror then
    printf('ERROR: %a\n', pinf);
    exit(1);
fi:
print(pinf);

# check that OutputFolder contains no spaces
for c in [' ','\n','\t','\c','\r'] do if SearchString(c, OutputFolder) >= 0 then
    error('found invalid character in ''OutputFolder'' argument (no spaces allowed).');
fi od:
OutputFolder := NormalizePath(OutputFolder);

if InputDataType='DNA' then
    ReadProgram(If(assigned(libdir),libdir,'lib') . '/DNA_matrix');
    CreateDayMatrices(DNA_Counts);
    DM := DayMatrix(50):
elif not assigned(InputDataType) or InputDataType='AA' then
    CreateDayMatrices();
    DM := DayMatrix(224):
else error('unexpected value for ''InputDataType'': '.str(InputDataType));
fi:

# this assinges DBs, GS, genomes, isContig and DB
ConvertGenomes(G, pinf);

NG := length(genomes):
if NG < 2 then
    printf('\n\n'
          .'%d genome found... make sure that your genomes (in fasta format)\n'
          .'are located in a directory called "DB", itself located in\n'
          .'the present directory.\n\n', NG);
    exit(1);
fi;

ValidateParameters():

if only_run_dbconv=true then 
    printf('\n\n*******\n');
    printf('*** You specified to stop after the database conversion step (i.e. you ***\n');
    printf('*** set the "-c" flag). Database conversion successfully finished.     ***\n');
    done
fi:


##############################################################################
# All-against-all
NumberFormat := '%.10g':
JustScoreBound := MinScore;
Set(gc=1e7):

NumberOfAlignments := proc(name1:string,name2:string)
    if name1=name2 then
        nrAlignments := GS[name1,TotEntries]*(GS[name1,TotEntries]-1)/2;
    else
        nrAlignments := GS[name1,TotEntries]*GS[name2,TotEntries];
    fi:
    return(nrAlignments);
end:

NumberOfChunks := proc(name1, name2)
    return(ceil(NumberOfAlignments(name1, name2)/AlignBatchSize));
end:

UpperTriangInd2Tuple := proc(first, last, N)
    tot := last-first+1;
    # find starting tuple
    c := first; m := 1;
    while c-N+m > 0 do
        c := c-N+m;
        m := m+1;
        if m > N then error('starting value out of range'); fi;
    od:
    out := CreateArray(1..tot);
    # generate all pairs
    for i to tot do:
        if m > N then error('end value out of range'); fi;
        out[i] := [m, c+m];
        c := c+1;
        if c+m > N then
            m := m+1;
            c := 1;
        fi;
    od;
    return(out)
end:

FullInd2Tuple := proc(first,last, N)
    tot := last-first+1;
    # find starting tuple
    m := floor((first-1)/N);
    c := first - m*N;

    out := CreateArray(1..tot);
    # generate all pairs
    for i to tot do:
        out[i] := [m+1, c];
        c := c+1;
        if c > N then
            m := m+1;
            c := 1;
        fi;
    od;
    return(out)
end:

RangeOfChunk := proc(name1,name2,chunk)
    totChunks := NumberOfChunks(name1, name2);
    if not type(chunk,posint) or chunk>totChunks then
        error('invalid chunk number: '.chunk)
    fi:
    all := NumberOfAlignments(name1, name2);
    first := floor((chunk-1)*all/totChunks)+1;
    last := floor(chunk*all/totChunks);
    if name1=name2 then
        # retrieve them in one go
        return(UpperTriangInd2Tuple(first,last,GS[name1,TotEntries]));
    else
        return( FullInd2Tuple(first,last,GS[name2,TotEntries]));
    fi:
end:



IsMyJob := proc(name1, name2, chunk)
    h := sha2(sprintf('%a',[name1,name2,chunk]))[1..8]:
    nr := sscanf(h, '%x')[1];
    return( evalb(mod(nr, pinf['NrParallelProc']) = pinf['ProcNr']-1) );
end:

CheckRemainingTimeAndQuitIfOver := proc(fn_chunk, matches, skip_until)
    if stop_at_time - UTCTime() < 30 then
        printf('\n\n%s\nReached timelimit. Quitting now...\n', CreateString(34,'*')):
        OpenWriting(fn_chunk.'.ckpt');
        printf('# checkpoint created by job %d (PID: %d, host %s) at %s\n',
            pinf['ProcNr'], getpid(), hostname(), date());
        printf('_matches := %A:\n_skip_until := %A:\n', matches, skip_until);
        OpenWriting(previous):
        exit(99):
    fi:
end:

if UseExperimentalHomologousClusters=true then
    ReadProgram(If(not assigned(libdir),'lib',libdir).'/HomologousClustering'):
else
# run standard all-against-all comparisons
tot := [seq(seq(NumberOfChunks(genomes[g1],genomes[g2]),g2=g1..NG), g1=1..NG)];
chunk_cnt := chunk_processed := genomepair_cnt := tot_align_time := 0:
for g1 to NG do
    for g2 from g1 to NG do
        genomepair_cnt := genomepair_cnt + 1;
        name1 := genomes[g1];
        name2 := genomes[g2];

        # Phase 0
        db1 := DBs[g1];
        db2 := DBs[g2];
        n1 := db1[TotEntries]:
        n2 := db2[TotEntries];
        if n1 > n2 or n1=n2 and name1 > name2 then
            t := db1;  db1 := db2;  db2 := t;
            t := name1;  name1 := name2;  name2 := t;
            t := n1;  n1 := n2;  n2 := t;
        fi:

        # skip if 1-chunk archive already present (import from OMA AxA)
        fn1chunk := sprintf('Cache/AllAll/%s/%s.gz', name1, name2);
        if FileExists(fn1chunk) then
            Logger(sprintf('Pair %s/%s already computed (probably imported). Skipped.\n',
                           name1, name2 ), mem=false);
            next;
        fi;

        Seqs1 := CreateArray(1..n1):
        Seqs2 := CreateArray(1..n2):
        DB := db1:
        for i1 to n1 do Seqs1[i1] := SearchTag('SEQ',Entry(i1)) od:
        DB := db2:
        for i2 to n2 do Seqs2[i2] := SearchTag('SEQ',Entry(i2)) od:

        # from release 0.99z.3 on, we store files in subdirectories per genome *pair*
        dir := sprintf('Cache/AllAll/%s/%s', name1, name2);
        if not FileExists(dir) then
            CallSystem('mkdir -p '.dir);
        fi:

        # Phase 1
        nrChunks := NumberOfChunks(name1, name2);
        for part to nrChunks do
            chunk_cnt := chunk_cnt+1;
            if not IsMyJob(name1,name2,part) then next fi:

            # skip if cached or already computed
            fn := sprintf('Cache/AllAll/%s/%s/part_%d-%d', name1, name2, part, nrChunks);
            if FileExists(fn.'.gz') then
                Logger(sprintf('Pair %s/%s part %d of %d already computed. Skipped.\n',
                    name1, name2, part, nrChunks), mem=false);
                next;
            elif FileExists(fn) then
                if not FileExists(fn.'.ckpt') then
                    Logger(sprintf('Pair %s/%s part %d of %d seems to be currently computed\n',
                        name1, name2, part, nrChunks), mem=false);
                    next;
                else
                    # load checkpoint
                    _matches := 0; _skip_until := 0:
                    err := traperror(ReadProgram(fn.'.ckpt'));
                    if err=lasterror or _matches=0 or _skip_until=0 then
                         # something bad happend. remove checkpoint file and start over with this chunk
                         res := []; skip_until := 0:
                         Logger(sprintf('Failed to load checkpoint file for %s/%s part %d of %d. Redo chunk\n',
                              name1, name2, part, nrChunks), mem=true);
                    else
                         res := _matches: skip_until := _skip_until:
                         Logger(sprintf('Loaded checkpoint file for %s/%s part %d of %d. Continue with this chunk\n',
                              name1, name2, part, nrChunks), mem=true);
                         Logger(sprintf(' loaded %d matches from checkpoint. will continue with pair %A',
                              length(res), skip_until),  2);
                    fi:
                    CallSystem('rm -f '.fn.'.ckpt');
                fi:
            else
                if FileExists(fn.'.ckpt') then
                    Logger(sprintf('Failed to load checkpoint file for %s/%s part %d of %d. No raw output. Redo chunk\n',
                        name1, name2, part, nrChunks), mem=true);
                    CallSystem('rm -f '.fn.'.ckpt');
                fi:
                res := []; skip_until := 0:
            fi:

            fn_oldstyle := sprintf('Cache/AllAll/%s/%s_%d-%d.gz', name1, name2, part, nrChunks);
            if FileExists(fn_oldstyle) then
                CallSystem('mv '.fn_oldstyle.' '.fn.'.gz');
                warning('using AllAll files from old release. moved file '.fn_oldstyle.'.gz into subfolder');
                next;
            fi:

            Logger(sprintf('Computing %s vs %s (Part %d of %d)',
                name1, name2, part, nrChunks));
            OpenWriting(fn);
            printf( '# AllAll of %s vs %s (Part %d of %d); pid %d; %s\n',
                name1, name2, part, nrChunks, getpid(), hostname() );
            OpenWriting(previous);

            t0 := UTCTime():
            dat0 := time():
            c2 := 0;
            algnInRange := min(AlignBatchSize,NumberOfAlignments(name1,name2));

            # iterate over all the protein pairs in this chunk and compute alignment
            for proteinPair in RangeOfChunk(name1, name2, part) do
                i1 := proteinPair[1];
                i2 := proteinPair[2];
                c2 := c2+1;
                # if we continue from a checkpoint, we skip until we reach the
                # proteinPair where we left off.
                if skip_until<>0 then
                      if proteinPair[2] = proteinPair[1] + 1 then
                          Logger(sprintf('  pair is %A. skip_until %A', proteinPair, skip_until), 4);
                      fi;
                      if skip_until = proteinPair then
                          skip_until := 0:
                          Logger('reached pair in Chunk where AllAll continues...', 2);
                      else next fi:
                fi:

                len1 := length(Seqs1[i1]);
                len2 := length(Seqs2[i2]);
                Logger(sprintf('aligning %s/%d vs %s/%d (seq lens: %d/%d)', 
                               name1, i1, name2, i2, len1, len2), 3);
                if mod(c2,109)=0 then
                    CheckRemainingTimeAndQuitIfOver(fn, res, proteinPair);
                fi:
                if mod(c2,50000)=0 then
                    Logger(sprintf('  %.2f%% complete, time left for this '.
                        'part=%.2fh, %.1f%% of AllAll done.', 100*c2/algnInRange,
                        (algnInRange/c2-1)*(UTCTime()-t0)/3600,
                        100*( (chunk_cnt-1+c2/algnInRange)/sum(tot) )));
                fi:
                dps := DynProgScore(Seqs1[i1],Seqs2[i2],DM,JustScoreShort);
                Logger(sprintf('fast sign check done. score: %.5g', dps[1]), 4);
                if dps[1] >= .75*MinScore then
                    a := Align(Seqs1[i1],Seqs2[i2],DMS);
                    Logger(sprintf('optimal alignment done. score: %f dist %f', a[Score], a[PamDistance]), 4);
                    if a[Score] < MinScore or min(a[Length1],a[Length2]) <
                                max(30,0.3*min(len1,len2)) then next fi;
                    i := GetOffset(a[Seq1])-GetOffset(Seqs1[i1]);
                    ran1 := i+1 .. i+a[Length1];
                    i := GetOffset(a[Seq2])-GetOffset(Seqs2[i2]);
                    ran2 := i+1 .. i+a[Length2];
                    res := append(res,[i1,i2,a[Score],a[PamDistance],
                                ran1,ran2,a[PamVariance]] );
                fi;
                if length(res)>1e5 then
                    WriteMatches(fn, res);
                    res := []:
                fi:
            od:
            WriteMatches(fn,res):

            drt := UTCTime()-t0:
            dat := time()-dat0:
            chunk_processed := chunk_processed + 1;
            tot_align_time := tot_align_time + drt;
            OpenAppending(fn);
            printf('# %s: %ds/%ds active = %.1f%%  %s\n',hostname(),round(dat),
                round(drt), 100*dat/drt,If(dat/drt<0.7,'SLOW!',''));
            OpenAppending(previous);
            # Compress file
            CallSystem('gzip -9f '.fn);
            msg := [sprintf('Finished %s vs %s (Part %d of %d)', name1, name2, part, nrChunks)];
            if c2 > 0 then
                msg := append(msg, sprintf('  %.1f%% of AllAll done, estimated remaining time: %.1fh\n', 
                    100*chunk_cnt/sum(tot), tot_align_time/chunk_processed * (sum(tot[genomepair_cnt..-1])+nrChunks-part)/3600));
            fi:
            Logger(msg, mem=false);
        od:
    od;
od:
fi:

## Parallel version:
## At this point, check whether all files are ready. If not
## exit (another process will pick from there on!)
f := TimedCallSystem('find Cache/AllAll -type f | grep -v ".gz" | grep -v -i ".DS_store" | grep -v ".ckpt"')[2]:
if length(f) > 0 then
    print('** At least 1 process appears to be still computing the all-vs-all.');
    print('** The following file(s) is (are) not yet completed:');
    prints(f);
    print('** If no other process is running, delete these files and restart.');
    quit;
fi:

if only_run_allall=true then
    print('\n\n*** All all-vs-all jobs successfully terminated.     ***');
    print('*** terminating after AllAll phase due to "-s" flag. ***');
    print('*** if you see this message at the end of one job,   ***'); 
    print('*** this means that all jobs successfully finished.  ***\n');
    quit;
fi:

##############################################################################
# Candidate Pairs
Set(gc=1e8):

print('\n\n****\n\nReading the all-against-all files...\n');

ns := CreateArray(1..NG):
Lengths := CreateArray(1..NG):
for i to NG do
    if GS[genomes[i]]='unassigned' then
        lprint(genomes[i],'missing GenomeSummary information');  quit fi;
    t := GS[genomes[i],EntryLengths];
    ns[i] := length(t);
    Lengths[i] := t
od:
t1 := Counter('RefinedMatches structures read'):
t2 := Counter('Matches read'):
t3 := Counter('files read'):
t4 := Counter('Matches of genomes against themselves'):
t5 := Counter('Matches above parameters');
t6 := Counter('empty all-all'):
t7 := Counter('Number of matches discarded below MinSeqLen');

# Read the BestMatch from all x all files #
BestMatch := CreateArray(1..NG,1..NG):
for i to NG do for j from i to NG do
    BestMatch[i,j] := CreateArray(1..ns[i],[]);
    if j>i then BestMatch[j,i] := CreateArray(1..ns[j],[]) fi;
    if ns[i] < ns[j] or ns[i]=ns[j] and genomes[i] < genomes[j] then
        g1 := i;  g2 := j else g2 := i;  g1 := j fi;
    BestMatch1 := BestMatch[g1,g2];  BestMatch2 := BestMatch[g2,g1];
    t3+1;
    Exclude1 := Exclude2 := {}:
    db1 := DBs[g1]; db2 := DBs[g2]:
    # if there exist a 1-chunk archive, process it with higher priority
    fn1chunk := sprintf('Cache/AllAll/%s/%s.gz',genomes[g1],genomes[g2]);
    if FileExists(fn1chunk) then
        fnSepHash := sprintf('Cache/AllAll/%s/%s.sha2.gz', genomes[g1], genomes[g2]);
        if FileExists(fnSepHash) then ReadProgram(fnSepHash) fi:
        err := traperror(ReadProgram(fn1chunk));
        if err=lasterror then error('corrupted AllAll file: '.fn1chunk.' Please remove and restart') fi;
    else
        nrChunks := NumberOfChunks(genomes[g1],genomes[g2]);
        for part to nrChunks do
            fn := sprintf('Cache/AllAll/%s/%s/part_%d-%d',genomes[g1],genomes[g2],
                part, nrChunks);
            err := traperror(ReadProgram( fn ));
            if err=lasterror then
               if lasterror='syntax error' then 
                   error('corrupted AllAll file: '.fn1chunk.' Please remove and restart');
               else error(lasterror) fi:
            fi:
        od:
    fi;
    # no Pairs should be cut for being below the Scoretol, as those
    # may serve to break a stable pair
    sij := sum(length(BestMatch[i,j,k]), k=1..length(BestMatch[i,j]));
    sji := sum(length(BestMatch[j,i,k]), k=1..length(BestMatch[j,i]));
    t5+max(sij,sji);
    if sij=0 or sji=0 then
        t6+1;
        # relax in case of mitochondion genome or synthetic evolution
        warning('empty all-all for '.genomes[g1].' '. genomes[g2] )
    fi;
    for k to length(BestMatch[i,j]) do
        if BestMatch[i,j,k] <> [] then
        BestMatch[i,j,k] := compress(BestMatch[i,j,k]) fi od;
    for k to length(BestMatch[j,i]) do
        if BestMatch[j,i,k] <> [] then
        BestMatch[j,i,k] := compress(BestMatch[j,i,k]) fi od;
    printf( '# %s vs %s done, (%d/%d)\n', genomes[g1],genomes[g2],
        sij, sji );
od od:
printf('# Matches loaded. Mem: %.3fGB\n', Set(BytesAlloc)/2^30);


##############################################################################
# Verification of stable pairs and clique building or esprit

cput_mid := time(); time_mid := UTCTime();
NumberFormat := '%.8g':
params := [MinScore, LengthTol, StablePairTol, VerifiedPairTol, UseEsprit]:
ParamSet := [];

if UseEsprit then
    print('\n\n****\n\nRunning ESPRIT for the detection of split-genes\n');
    Esprit([MinProbContig, MinBestScore, DistConfLevel,
            true, MaxContigOverlap]):
    WriteEspritResults();
else
    if UseOnlyOneSplicingVariant=true then
        FilterBestSplicingMatches();
    fi:

    # check whether orthologs have already been computed
    orthoCacheFile := 'Cache/ortholog_milestone.drw.gz':
    alreadyComputed := false;
    if FileExists(orthoCacheFile) then
        alreadyComputed := LoadOrthologMilestone(orthoCacheFile):
    fi:

    if not alreadyComputed then
        print('\n\n****\n\nIdentification of Stable-Pairs (Putative Orthologs)\n');
        Ort := OrthologousMatrix( op(params) ):
        StoreOrthologMilestone( orthoCacheFile );
    fi:

    if not DoHierarchicalGroups=false or (
           DoGroupFunctionPrediction=true and CladeDefinition='default') then
        if SpeciesTree='estimate' then;
            speciesTree := BuildSpeciesTree( Ort );
            if OutgroupSpecies <> 'none' then
                speciesTree := RerootWithOutgroups(speciesTree, OutgroupSpecies);
            fi:
        elif SpeciesTree='lineage' then
            speciesTree := NULL:
        elif type(SpeciesTree,string) then
            speciesTree := ParseNewickTree(SpeciesTree,
                'defaultBranchLength'=1):
        fi:
    fi:

    if DoHierarchicalGroups='top-down' then
        # call hierarchical group procedure
        HOG := traperror(HierarchicalGroups(ReachabilityCutoff, speciesTree));
        if HOG=lasterror then
            warning('computation of HierarchicalGroups failed: '.string(lasterror));
        fi:
    fi:

    CallSystem('rm -rf '.OutputFolder.'; mkdir -p '.OutputFolder);

    # # Retrieve all IDs
    IDtable := CreateArray(1..NG);
    for i to NG do
        DB := DBs[i];
        tmp := CreateArray(1..ns[i]);
        for j to DB[TotEntries] do
            e := Entry(j):
            tmp[j] := SearchTag('FASTAHEADER', e); 
            if tmp[j]='' then tmp[j] := SearchTag('ID', e); fi:
        od;
        IDtable[i] := tmp;
    od:

    # store used splicing variant
    main_isoform := CreateString(sum(ns), '1'):
    Goff := CreateArray(1..NG+1): 
    for i to NG do Goff[i+1] := Goff[i] + ns[i] od:
    assert(length(main_isoform) = Goff[NG+1]);

    if UseOnlyOneSplicingVariant=true then
        used_variants := []:
        for i to NG do
            g := genomes[i];
            splicemap := LoadSpliceMap(g);
            for alt in splicemap do
                NM := [seq(sum(length(decompress(BestMatch[i,j,k])),j=1..NG),k=alt)];
                best := SearchArray(max(NM), NM);
                used_variants := append(used_variants, [g, IDtable[i,alt[best]]]);
                for j to length(alt) do if j<>best then 
                    main_isoform[Goff[i]+alt[j]] := '0';
                fi od:
            od:
        od:
        if length(used_variants)>0 then
            OpenWriting(OutputFolder.'/used_splicing_variants.txt');
            for z in used_variants do
                printf('%s\t%s\n', op(z));
            od:
            OpenWriting(previous);
        fi:
    fi:


    # # Create File MapSeqNum-ID:
    OpenWriting(OutputFolder.'/Map-SeqNum-ID.txt');
    print('\#Format: genome<tab>sequence number<tab>id<tab>is_main_isoform');
    for i to NG do for j to ns[i] do
        printf('%s\t%d\t%s\t%s\n', genomes[i], j, IDtable[i,j], If(main_isoform[Goff[i]+j]='1', 'True', 'False'));
    od od:
    OpenWriting(previous);

    if DoHomologyFunctionPrediction=true then
        PredictViaHomology();
    fi:

    if DoGroupFunctionPrediction=true then
        can_do_function := true;
        if SpeciesTree='lineage' then
            warning('Annotation propagation is not implemented for parameter SpeciesTree := ''lineage'''):
            can_do_function := false;
        elif not CheckExported(GroupFunctionCutoff) then
            printf('Too few OMA exported genomes are included to do group function prediction.\n'):
            printf('To run either include more exported genomes or override by setting \n'):
            printf('GroupFunctionCutoff lower in the parameters file.\n'):
            warning('Too few OMA exported genomes for function prediction.'):
            can_do_function := false;
        fi:

        if can_do_function then 
            if CladeDefinition<>'default' then
                new_annots := PropagateFunctionByGroup(Ort, CladeDefinition, GroupFunctionCutoff):
            else
                new_annots := PropagateFunctionByGroup(Ort, speciesTree, GroupFunctionCutoff):
            fi:
        fi:
    fi:

    ############################################################################
    # Store results
    Set('printgc'=false);
    printf('\n\n%s\n\nStoring results\n', CreateString(80,'*'));

    if lasterror <> 'lasterror' then
        warning('An error occured during inference: '.string(lasterror));
        # quit;
    fi:

    OGtable := CreateArray(1..NG):
    for i to NG do OGtable[i] := CreateArray(1..ns[i]); od:
    for r to length(Ort) do
        for k to NG do
            t := Ort[r,k];
            if t <> 0 then
                OGtable[k,t] := r;
            fi;
        od;
    od:

    # Store the pairwise orthologs
    if not WriteOutput_PairwiseOrthologs=false or DoHierarchicalGroups='bottom-up' then
        CallSystem('mkdir '.OutputFolder.'/PairwiseOrthologs');
        for i to NG do for j from i+1 to NG do
            if ns[i] < ns[j] or ns[i]=ns[j] and genomes[i] < genomes[j] then
                g1 := i;  g2 := j else g2 := i;  g1 := j fi;
            OpenWriting(OutputFolder.'/PairwiseOrthologs/'.genomes[g1].'-'.genomes[g2].'.txt');
            printf('\# Format: Protein 1<tab>Protein 2<tab>Protein ID1<tab>ProteinID2'
                 .'<tab>Orthology type<tab>OMA group (if any)\n');
            printf('\# Every pair is listed only once, and in no particular order.\n');
            printf('\# The map between sequence number and ID are given\n');
            printf('\# in the file "Map-SeqNum-ID.map\n');

            for k to ns[g1] do
                vp := decompress(VPairs[g1,g2,k]);
                typ2 := If(length(vp)>1,'many','1');
                for l in vp do
                    # if only 1 in target -> :1, else :many
                    # l
                    typ1 := If(length(decompress(VPairs[g2,g1,l]))>1,'many','1');
                    printf('%d\t%d\t%s\t%s\t%s:%s\t%s\n',
                        k,l,
                        IDtable[g1,k],
                        IDtable[g2,l],
                        typ1, typ2,
                        If(OGtable[g1,k]=OGtable[g2,l] and OGtable[g1,k] <> 0,
                            string(OGtable[g1,k]),
                            '')
                    );
                od;
            od:
            OpenWriting(previous);
        od od:
    fi:

    if not WriteOutput_Paralogs=false and ParalogTol>0 then
        CallSystem('mkdir '.OutputFolder.'/PairwiseParalogs');
        for i to NG do for j from i+1 to NG do
            if ns[i] < ns[j] or ns[i]=ns[j] and genomes[i] < genomes[j] then
                g1 := i;  g2 := j else g2 := i;  g1 := j fi;

            OpenWriting(OutputFolder.'/PairwiseParalogs/'.genomes[g1].'-'.genomes[g2].'.txt');
            printf('#nr1\tnr2\tid1\tid2\tScore\tPamDist\n');
            for x1 to ns[g1] do 
                pp := decompress(ParalogPairs[g1,g2,x1]);
                for y2 in pp do 
                    printf('%d\t%d\t%s\t%s\t%.2f\t%.2f\n',
                        x1, y2['Entry'], IDtable[g1,x1], IDtable[g2, y2['Entry']],
                        y2['Score100']/100, y2['PamDist10000']/10000);
                od:
            od:
            OpenWriting(previous);
        od;od:
    fi:
    # store the pairwise orthologs in orthoxml format
    if not WriteOutput_OrthologousPairs_orthoxml=false then
        WriteOrthoXML( VPairs, OutputFolder.'/OrthologousPairs.orthoxml');
    fi:

    # Store the groups of orthologs
    WriteMatrix(Ort, OutputFolder.'/OrthologousMatrix.txt');
    WriteGroups(Ort, OutputFolder.'/OrthologousGroups.txt');
    # store orthologous groups in orthoxml format
    WriteOrthoXML(Ort, OutputFolder.'/OrthologousGroups.orthoxml',
        StableIdsForGroups=true);
    if not WriteOutput_PhyleticProfileOG = false then
        WriteProfileFromOGs(Ort, OutputFolder.'/PhyleticProfileOMAGroups.txt');
    fi:

    if type(QuasiCliquesCutoff, numeric) and QuasiCliquesCutoff<1 then
        mergedOMAGroups := MergeOMAGroups(Ort, VPairs, QuasiCliquesCutoff);
        WriteOrthoXML( mergedOMAGroups, OutputFolder.'/MergedOMAGroups.orthoxml',
            StableIdsForGroups=true);
        WriteMatrix( mergedOMAGroups, OutputFolder.'/MergedOMAMatrix.txt');
        WriteGroups( mergedOMAGroups, OutputFolder.'/MergedOMAGroups.txt');
        if not WriteOutput_OrthologousGroupsFasta=false then
            WriteFastaGroups(OutputFolder.'/MergedOMAGroupsFasta', mergedOMAGroups);
        fi:
        if not WriteOutput_PhyleticProfileOG = false then
            WriteProfileFromOGs(mergedOMAGroups, OutputFolder.'/PhyleticProfileMergedOMAGroups.txt');
        fi:
    fi:

    # store all orthologous groups in fasta format
    if not WriteOutput_OrthologousGroupsFasta=false then
        WriteFastaGroups(OutputFolder.'/OrthologousGroupsFasta', Ort):
    fi:

    # store hierarchical groups if computed
    if member(DoHierarchicalGroups, {'top-down','bottom-up'}) then
        if SpeciesTree='estimate' then 
            OpenWriting(OutputFolder.'/EstimatedSpeciesTree.nwk');
            prints( Tree_Newick(speciesTree) .';' );
            OpenWriting(previous):
            WritePhyloXML(speciesTree, OutputFolder.'/EstimatedSpeciesTree.phyloxml');
        fi:
        if DoHierarchicalGroups='bottom-up' then
            HOG := GetHOGsBottomUp(speciesTree):
        elif DoHierarchicalGroups='top-down' and type(HOG,list) then
            WriteOrthoXML( HOG, OutputFolder.'/HierarchicalGroups.orthoxml',
                StableIdsForGroups=true);
        fi:

        if not WriteOutput_HOGFasta=false then
            WriteFastaGroups(OutputFolder.'/HOGFasta', HOG);
        fi:
        if not WriteOutput_PhyleticProfileHOG = false then
            WriteProfileFromHOGs(HOG, OutputFolder.'/PhyleticProfileHOGs.txt');
        fi:
    fi:

    # Write out function predictions
    if DoGroupFunctionPrediction=true and can_do_function then
        WriteGAF(new_annots, OutputFolder.'/gene_function.gaf', GroupFunctionCutoff):
    fi:

    ###
    # Final comments to user
    printf('\n\n%s\n', CreateString(80, '*'));
    if SpeciesTree='estimate' and member(DoHierarchicalGroups, {'top-down','bottom-up'})
            and OutgroupSpecies='none' then
        printf('NOTE: The HOGs OMA computed depend strongly on the correct species tree.\n');
        printf('      OMA estimated a species tree based on some marker genes, and stored\n');
        printf('      it in %s.\n', OutputFolder.'/EstimatedSpeciesTree.nwk');
        printf('      Please verify that it is in agreement with your expectations, especially\n');
        printf('      that it is rooted correctly. If not, please modify it and put the\n');
        printf('      corrected tree in the parameter file and rerun OMA. You can visualize and\n');
        printf('      reroot the tree for example with http://phylo.io\n');
    fi:
fi:

printf('\n\n%s\n\nDone!!\n\n', CreateString(80, '*'));