SpeciesTreeDiscordanceTest.drw

##
#  Script to run the species tree discordance test  
#
#  in: resFn, projs, lock
Set(printgc=false):
SetRandSeed();

WriteLock(lock, 'initialize Species Tree discordance test compuation');
CreateDayMatrices();
NR_TRIALS := 50000;
MAX_NR_TREES := 3000;

prj := [seq(GetPKeyData(z),z=projs)]:
dataset := prj[-1,'ref'];
hashOfAllProjPIDs := hash_sha2(projs):

if not assigned(confidence) then
    confidence := NULL;
elif type(confidence, posint) then
    confidence := '_conf'.confidence;
else error('unexpected type for confidence parameter');
fi:

if not member(treebuilder, {'LSTree','BIONJ'}) then 
    error('unexpected treebuilder method: '.treebuilder);
fi;


# load species tree as nested list of list, e.g.
# sptree := [[['HUMAN','PANTR'],['MOUSE','RATNO'],'CANFA'],'CHICK','XENTR'];
# i.e. the tree is non-binary
ReadProgram(datdir.'speciestree_'.problem.'_'.dataset.string(confidence).'.drw'); 

ShuffleTreeOrder := proc(t)
    cpy := copy(t):
    for i to length(cpy) do if type(cpy[i],list) then
        cpy[i] := procname(cpy[i]);
    fi od:
    return( Shuffle(cpy) ):
end:

SelectAcceptedBinarySubTree := proc(t)
    res := NULL;
    if type(t,{string,symbol}) then res := t;
    else 
        tL := procname(t[1]);
        tR := procname(t[2]);
        res := tL,tR;
    fi:
    return(res);
end:

GetAcceptedSpeciesTopology := proc(t;(h0=0):numeric)
    res := NULL;
    if type(t,{string,symbol}) then res := Leaf(t,h0);
    else 
        tL := procname(t[1],h0+1);
        tR := procname(t[2],h0+1);
        res := Tree(tL,h0,tR);
    fi:
    return(res);
end:

TreeToConsecExprSeq := proc(t)
    if type(t,Leaf) then return(t['Label']);
    else
        return( procname(t['Left']), procname(t['Right']) );
    fi:
end:


GetSeqOfOrthologs := proc( curProt, remSp )
    global get_orth_calls, fail_fast;
    get_orth_calls := get_orth_calls + 1;

    #printf('calls: %d; -->%d, %A\n', get_orth_calls, curProt, remSp):
    res := NULL;
    vps := ParseLongList(SearchTag('VP',Entry(curProt)));
    nextSpeciesRange := GenomeRange(remSp[1], dataset);
    idx1 := SearchOrderedArray(nextSpeciesRange[1]-1, vps)+1;
    idx2 := SearchOrderedArray(nextSpeciesRange[2], vps);
    cand_vps := Shuffle(vps[idx1..idx2]);
    if length(cand_vps)>20 then cand_vps := cand_vps[1..20] fi:

    for vp in cand_vps do if fail_fast[vp]='0' and member(vp,iSet) then
        if vp>=nextSpeciesRange[1] and vp<=nextSpeciesRange[2] then
            if length(remSp)=1 then res := vp;
            else 
                res := remember(procname( vp, remSp[2..-1]));
            fi:
            if res<>NULL then break else fail_fast[vp] := '1'; fi:
        fi:
    fi od:
    if res<>NULL then 
        res := curProt, res;
    else
        fail_fast[curProt] := '1';
    fi:
    #printf(' <--%d, %A: %A\n', curProt, remSp, If(res<>NULL, res, 'n/a')):
    return( res );
end:


SelectNextFam := proc()
    global fail_fast, get_orth_calls:
    t0 := time();
    maxTry := 10;
    for x to maxTry do 
        t := ShuffleTreeOrder(sptree):
        accTree := GetAcceptedSpeciesTopology(t);
        specList := [TreeToConsecExprSeq(accTree)];
        if length(specList)>3 then break fi:
    od:
    if x>maxTry then return(NULL) fi:

    startRange := GenomeRange(specList[1], dataset);
    fnd := false;
    for xx to 50 do 
        start := Rand(startRange);
        if member(start,iSet) then fnd := true; break fi:
    od:
    if not fnd then 
        warning('no start protein in '.specList[1].' found in iSet');
        return(NULL) 
    fi:
    t1 := time();
    fail_fast := CreateString(SeqDB[TotEntries], '0'):
    get_orth_calls := 0;
    orths := GetSeqOfOrthologs(start, specList[2..-1] );
    if orths<>NULL then res := [accTree, [orths]];
    else res := NULL fi:
    #print('%d calls to get_orthologs. sampling completed: %A\n', get_orth_calls, nxt<>NULL);
    return( res );
end:

FindCases := proc(iSet, prjNrs:list)
    global DB;
    nP := length(prjNrs);
    cases := CreateArray(1..nP,[]):
    NrOrthologs := CreateArray(1..nP):
    protSet := intersect(iSet, GetProtsInSPtree()):

    for i2pNr to nP do
        pNr := prjNrs[i2pNr];
        assert( prj[pNr,'directMap'] );
        DB := dbs[pNr];
        
        msg := 'sampling possible trees for '.prj[pNr,'Title'];
        WriteLock(lock, msg);
        cand := []: sampling_times := Stat('Stats of sampling times for '.prj[pNr,'Title']);
        tLast := tBegin := time();
        for trial to NR_TRIALS do
            tCur := time();
            nxt := SelectNextFam();
            if nxt <> NULL then cand := append(cand, nxt) fi:
            sampling_times + (time() - tCur);
            if time()-tLast > 30 then
                print(sampling_times);
                WriteLock(lock, sprintf('%s\n %d of %d samplings done; %d samplings are complete.\n'.
                                        '  sampling for %.1fmin (%.3fsec/sampling)\n\n', 
                                        msg, trial, NR_TRIALS, length({op(cand)}), (time()-tBegin)/60,
                                        (time()-tBegin)/trial));
                tLast := time();
            fi:
        od:
        cand := [op({op(cand)})]:
        cases[i2pNr] := cand;
        NrOrthologs[i2pNr] := CntOrthologs( protSet ):
    od:
   
    Logger( sprintf('nr of cases in total: %d',sum(length(z),z=cases)), 'INFO');
    return(cases, NrOrthologs):
end:

GetProtsInSPtree := proc()
    genomesInTree := indets(sptree, string):
    protSet := {seq(seq(eNr, eNr=GenomeRange(g,dataset)),g=genomesInTree)};
    return( protSet );
end:

CntOrthologs := proc(iSet)
    # precondition: project has directMap
    cnt := 0;
    for eNr to DB[TotEntries] do
        if not member(eNr,iSet) then next fi:
        vps := ParseLongList(SearchTag('VP',Entry(eNr)));
        for vp in vps do if member(vp, iSet) then cnt := cnt+1 fi od:
    od:
    return( cnt/2 );
end:

InferDistTree := proc(D, V, labs)
    if treebuilder='LSTree' then
        tree := LeastSquaresTree(D,V,labs);
    elif treebuilder='BIONJ' then
        treeRes := BioNJ(D, labs);
        tree := treeRes['Tree'];
    else 
        error('invalid tree reconstruction method: '.treebuilder) 
    fi:
    return(tree);
end:

ComputeCongruenceStat := proc(prjNrs:list(posint), cases:list, totNrTrees:{0,posint})
    global DB, Ndone:
    RFstat := [seq( Stat(prj[i,'Title']), i=prjNrs)];
    IdentStat := [seq( Stat(prj[i,'Title']), i=prjNrs)];
    rawDat := CreateArray(1..length(prjNrs),[]):
    
    DB := SeqDB: 
    cErr := Counter('# errors in PhylogeneticTree');
    if not type(Ndone, integer) then Ndone:=0 fi:

    tLast := 0;
    for pNr to length(prjNrs) do
        for cas in cases[pNr] do
            seqs := labs := spec := [];
            speciestree := cas[1];
            for i to length(cas[2]) do 
                eNr := cas[2,i];
                seqs := append(seqs, Sequence(Entry(eNr)));
                labs := append(labs, eNr);
                spec := append(spec, GenomeCodeFromENr(eNr,dataset));
            od:
            nseq := length(seqs):
            msa := traperror(MafftMSA(seqs)):
            if msa=lasterror then cErr+1; next fi:
            D := CreateArray(1..nseq,1..nseq):
            V := CreateArray(1..nseq,1..nseq):
            for i to nseq do for j from i+1 to nseq do
                 dps := EstimatePam(msa[AlignedSeqs,i],msa[AlignedSeqs,j],DMS):
                 D[i,j] := D[j,i] := dps[2];
                 V[i,j] := V[j,i] := dps[3];
            od od:
            genetree := traperror( InferDistTree(D,V,labs) ):
            if genetree=lasterror then cErr+1; next fi;
            
            pruned_st := PruneTree(copy(speciestree), spec);
            for l in Leaves(pruned_st) do 
                x := SearchArray(l['Label'], spec);
                l['Label'] := labs[x];
            od:

            rfdist := RobinsonFoulds([genetree, pruned_st])[1,2];
            RFstat[pNr] + rfdist;
            IdentStat[pNr] + If(rfdist>0,1,0);
            rawDat[pNr] := append(rawDat[pNr], [genetree, pruned_st, rfdist]);
            Ndone := Ndone + 1;
            if time()-tLast>30 then
               t := Ndone/totNrTrees;
               msg := sprintf( 'Finished %d / %d (%.1f%%) trees in %.1f min. '.
                   'Estimated remaining time: %.1fmin\n', 
                   Ndone, totNrTrees, 100*t, (time()-t0)/60, (1-t)/t*(time()-t0)/60 );
               WriteLock(lock,msg);
               tLast := time();
            fi;
        od;
    od:
    Logger( sprint('STDTest result for %s:\n  %A\n  %A', 
        ConcatStrings([seq(prj[i,'Title'],i=prjNrs)],' - '), op(RFstat), cErr), 'INFO');
    return( [RFstat, rawDat, IdentStat] );
end:

StoreResult := proc(fn:string, RFstat, IdentStat, NrOrthologs, NrTrees, RawDataFn, PIDs)
    OpenWriting(fn):
    printf('RFstat := %A:\n', RFstat):
    printf('PIDs := %A:\n', PIDs );
    printf('RawDataFn := %A:\n', RawDataFn);
    printf('IdentStat := %A:\n', IdentStat);
    printf('NrOrthologs := %A:\n', NrOrthologs);
    printf('NrTrees := %A:\n', NrTrees);
    printf('#problem := %A: confidence := %A: treebuilder := %A:\n', 
           problem, confidence, treebuilder);
    OpenWriting(previous):
    CallSystem('gzip -9f '.fn);
end:

StoreRawData := proc(fn:string, RawData, projects)
    OpenWriting(fn);
    printf('# Project <TAB> gene tree <TAB> species tree <TAB> RF distance\n');
    for pNr to length(projects) do 
        pName := projects[pNr,'Title'];
        for z in RawData[pNr] do
            for k to 2 do for l in Leaves(z[k]) do 
                oE := l['Label']; 
                if not type(oE,integer) then next fi:
                l['Label'] := ENr2XRef(oE, dataset); 
            od od:
            printf('%s\t%s\t%s\t%f\n', pName, Tree_Newick(z[1]), 
                Tree_Newick(z[2]), z[3]);
        od:
    od:
    OpenWriting(previous);
    CallSystem('gzip -9f '.fn);
end:

# load the project databases
nP := length(prj); dbs := CreateArray(1..nP):
for i to nP do dbs[i] := DB := ReadDb( prj[i,'fnBase'].'.db' ); od:

# load index database
ddir := eval(symbol(lowercase(dataset).'DBpath'));
SeqDB := ReadDb( ddir.'/ServerSeqs.db');


if length(meth)>=6 and meth[1..6]='REFSET' then
    cases := CreateArray(1..nP,[]):
    RFstats := CreateArray(1..nP):
    IdentStats := CreateArray(1..nP):
    NrOrthologList := CreateArray(1..nP):
    NrTreesList := CreateArray(1..nP):
    RawData := CreateArray(1..nP):
    prjNrs := [seq(i,i=1..nP)];
    resFnParts := CreateArray(1..nP):
    resFnBase := resFn[1..(SearchAllString('.',resFn)[-2])];
    rawDataFn := CreateArray(1..nP):
    iSet := intersect():
    for i to nP do 
        hash_of_cur_proj := hash_sha2([problem, prj[i,'PID'], meth, confidence, treebuilder]);
        resFnParts[i] := sprintf('%s%a.drw', resFnBase, hash_of_cur_proj); 
        rawDataFn[i] := sprintf('STDTest/STD_%s.%a.tsv.gz', meth, hash_of_cur_proj); 
        if FileExists(resFnParts[i].'.gz') then next fi:
        
        casesAndCnts := FindCases( intersect(), [i] );
        cases[i] := casesAndCnts[1,1]; NrOrthologList[i] := casesAndCnts[2,1];
        NrTreesList[i] := length(cases[i]);
        cases[i] := If(NrTreesList[i]>MAX_NR_TREES, Shuffle(cases[i])[1..MAX_NR_TREES], cases[i]):
    od:
    t0 := time();
    totNrTrees := sum(length(z),z=cases):
    for i to nP do
        if FileExists(resFnParts[i].'.gz') then
            RFstat := RawDataFn := 0;
            ReadProgram(resFnParts[i]): 
            assert(type(RFstat,Stat)); 
            assert(type(IdentStat,Stat)); 
            assert(type(RawDataFn,string));
            RFstats[i] := RFstat;
            IdentStats[i] := IdentStat;
            rawDataFn[i] := RawDataFn;
            NrOrthologList[i] := NrOrthologs;
            NrTreesList[i] := NrTrees;
        else
            compRes := ComputeCongruenceStat([prjNrs[i]], [cases[i]], 
                totNrTrees):
            # compRes[1]: RFstat ; compRes[2]: raw data to be stored away
            RFstats[i] := compRes[1,1]; rawdata := compRes[2,1]; IdentStats[i] := compRes[3,1];
            StoreRawData(rawdir.rawDataFn[i,1..-4], [rawdata], [prj[i]]);
            StoreResult(resFnParts[i], RFstats[i], IdentStats[i], NrOrthologList[i], NrTreesList[i], rawDataFn[i], 
                [prj[i,'PID']]);
        fi:
    od:
    StoreResult(resFn, RFstats, IdentStats, NrOrthologList, NrTreesList, rawDataFn, [seq(prj[i,'PID'],i=1..nP)]);
else error('unknown strategy method:'.meth); 
fi: 

CallSystem('rm -f '.lock);
done;