HomologyTest.drw

# 
# A script to measure missing and fp homolog relations of 
# multidomain proteins based on the method and data of 
# Dannie Durand (Song et al, PLOS Comp Biol, 2008)
#
#                             Adrian Altenhoff, Jun 8, 2011  
#
#  in: resFn, projs, problems, measure, lock, 
Set(printgc=false); Set(gc=5e6):
WriteLock(lock,'initialize test data');

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

# 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:
if not type(dbs, list(database)) then 
   OpenWriting(resFn);
   printf('error(''could not open databases'');\n');
   OpenWriting(previous);
   CallSystem('rm -f '.lock);
   quit
fi:

ReadProgram(datdir.'homologyTest_'.dataset.'.drw'); # assigns pos/neg
nProb := length(problems):
WriteLock(lock, 'computing preformance...');

#ddir := eval(symbol(lowercase(dataset).'DBpath')):
#SeqDB := ReadDb(ddir.'ServerSeqs.db'):
ReadProgram(datdir.'homology.mapback'):
InitLogger := proc()
    global LoggedPairs:
    LoggedPairs := table([],[]):
end:
LogPairs := proc(pairs,typ)
    global LoggedPairs:
    t := []:
    for z in pairs do
        k1 := SearchOrderedArray( z[1], oENr2jacob[2] ); 
        if k1<0 or k1>length(oENr2jacob[2]) or oENr2jacob[2,k1]<>z[1] then 
            Logger( sprintf('LogPairs: should not happen z[1]: %a',z[1]), 'ERROR');
        else id1 := oENr2jacob[1,k1]; fi:
        k2:= SearchOrderedArray( z[2], oENr2jacob[2] ); 
        if k2<0 or k2>length(oENr2jacob[2]) or oENr2jacob[2,k2]<>z[2] then 
            Logger( sprintf('LogPairs: should not happen z[2]: %a',z[2]), 'ERROR');
        else id2 := oENr2jacob[1,k2]; fi:
        t := append( t, [id1, id2] );
    od:
    LoggedPairs[typ] := append(LoggedPairs[typ], op(t)):
    DB := curDB:
end:
GetLoggedData := proc(typ)
    return( LoggedPairs[typ] ):
end:

# pairs is list of length NrOfProteins containing a set of 
# homologous/non-homologous partners
ComputeRates := proc(pos, neg, iSet, directMap:boolean)
    fp := tp := covPos := covNeg := 0; 
    InitLogger():
    for i to DB[TotEntries] do 
        e := Entry(i);
        oE1 := If(directMap, i, parse(SearchTag('OE',e))):
        if not member(oE1, iSet) then next fi:
        
        if length(pos[oE1])>0 or length(neg[oE1])>0 then
            vps := {op(ParseLongList(SearchTag('VP',e)))}:
            p := length( intersect(vps, pos[oE1]) );
            n := length( intersect(vps, neg[oE1]) );
            fp := fp+n; tp := tp+p; covPos := covPos+p; covNeg := covNeg+n;
            if n>0 then LogPairs( [seq([oE1,z],z=intersect(vps, neg[oE1]))],'FP') fi:
            if p>0 then LogPairs( [seq([oE1,z],z=intersect(vps, pos[oE1]))],'TP') fi:
        fi:
    od:
    FPdata := GetLoggedData('FP'):
    TPdata := GetLoggedData('TP'):
    return( [[tp,0], EstimateAndVarOfFraction(fp, covPos+covNeg), FPdata, TPdata] );
end:

EstimateAndVarOfFraction := proc( tp:integer, tot:integer )
    if tot=0 then return([0,-1]);
    else 
        p := tp/tot;
        return( [p, p*(1-p)/tot] );
    fi:
end:

ComputePerformance := proc(iSet, prjNrs:list(posint))
    global DB;
    np := length(prjNrs);
    fpr := CreateArray(1..nProb,1..np,0);
    tpr := CreateArray(1..nProb,1..np,0);
    FPdata := CreateArray(1..nProb,1..np,0);
    TPdata := CreateArray(1..nProb,1..np,0);
    for prNr to nProb do for pNr to nP do
        DB := dbs[prjNrs[pNr]];
        directMap := prj[prjNrs[pNr],'directMap'];
        
        t := ComputeRates( pos[problems[prNr]], neg[problems[prNr]], iSet, directMap);
        fpr[prNr,pNr] := t[1];
        tpr[prNr,pNr] := t[2];
        FPdata[prNr,pNr] := t[3];
        TPdata[prNr,pNr] := t[4];
    od od:
    return( [fpr, tpr, FPdata, TPdata] );
end:

StoreResult := proc(fn:string, FPR, TPR, PIDs, measure, problems, meth, FPdata, TPdata)
    OpenWriting(fn):
    printf('# created on %s\n', date());
    printf('FPR := %A:\nTPR := %A:\n', FPR, TPR):
    printf('#PIDs := %A:\n', PIDs );
    printf('#problems:=%A:\nmeasure:=%A:\n#meth := %A:\n', problems, measure, meth);
    printf('#FPdata := %A:\n#TPdata := %A:\n', FPdata, TPdata);
    OpenWriting(previous):
    CallSystem('gzip -9f '.fn);
end:
          
if meth='INTERSECT' or (length(meth)>=6 and meth[1..6]='REFSET') then
    iSet  := GetProjISets(prj):
    prjNrs:= [seq(i,i=1..nP)];
    Res := ComputePerformance(iSet, prjNrs):
    StoreResult(resFn, Res[1], Res[2], [seq(prj[i,'PID'], i=prjNrs)], 
        problems, measure, meth, Res[3], Res[4]);
elif meth='PAIRWISE' then
    maxP := nP-1;
    SimsRes := CreateArray(1..maxP):
    prjNrs := [seq([i,nP],i=1..maxP)];
    resFnBase := resFn[1..(SearchAllString('.',resFn)[-2])];
    for curP to maxP do 
        curPrjPIDs := [seq(prj[j,'PID'],j=prjNrs[curP])]:
        resFnPart := sprintf('%s%a.drw', resFnBase, hash_sha2([curPrjPIDs, problems, 
            measure,'INTERSECT']) );
        if FileExists(resFnPart.'.gz') then 
            Res[curP] := 0; ReadProgram(resFnPart): Res[curP] := [FPR,TPR];
        else
            iSet := GetProjISets( [seq(prj[j],j=prjNrs[curP])] ):
            Res[curP] := ComputePerformance(iSet, prjNrs[curP]):
            StoreResult(resFnPart, Res[curP,1], Res[curP,2], 
                [seq(prj[j,'PID'], j=prjNrs[curP])], problems, measure,
                'INTERSECT',Res[curP,3], Res[curP,4]);
        fi:
    od:
    StoreResult(resFn, SimsRes, [seq([prj[i,'PID'],prj[nP,'PID']],i=1..maxP)], 
        problems, measure, meth);
#elif length(meth)>=6 and meth[1..6]='REFSET' then
#    Res := CreateArray(1..nP):
#    resFnBase := resFn[1..(SearchAllString('.',resFn)[-2])];
#    maxP := 1:
#    for curP to nP do 
#        resFnPart := sprintf('%s%a.drw', resFnBase, hash_sha2([prj[curP,'PID'],problems, 
#            measure,meth]) );
#        if FileExists(resFnPart.'.gz') then 
#            Res[curP] := 0; ReadProgram(resFnPart): Res[curP] := [FPR, TPR];
#        else
#            Res[curP] := ComputePerformance(intersect(), [curP])[1]:
#            StoreResult(resFnPart, Res[curP,1], Res[curP,2], [prj[curP,'PID']],
#                problems, measure, meth);
#        fi:
#    od:
#    StoreResult(resFn, Res[1], Res[2], [seq(prj[i,'PID'], i=1..nP)], problems,
#        measure, meth);
#
else error('unknown strategy method:'.meth); 
fi:


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