Synthetic Evolution Exercise - Solution

CreateAdjustedParamFile := proc( origParams:string,  invSites:nonnegative )
    fn := sprintf('/tmp/param.%.2f.drw', invSites);
    param := copy(origParams):
    varName := 'motifFreq': lenVarName := length(varName):

    # we search the variable and replace its value
    off := SearchString(varName, param):
    for k1 from off to length(param) while param[k1]<>'=' and param[k1-1]<>':' do od:
    for k2 from k1+1 to length(param) while param[k2]<>':' and param[k2]<>';' do od:
    param := param[1..k1+1].string(invSites).param[k2..-1]:

    OpenWriting( fn ):
    printf('%s\n', param):
    OpenWriting( previous ):
    return( fn ):
end:

RunSim := proc( origParams:string, invSites:nonnegative )
    # we want to have a float value between 0 and 1 representing the percentage 
    # of invariant sites
    if invSites >= 1 then error('0 <= invSites < 1, was:'.string(invSites) ) fi;
    
    # generate synthetic sequences
    paramFile := CreateAdjustedParamFile( origParams, invSites );
    printf('File is %s\n', paramFile):
    eval(parse(ReadRawFile(paramFile))):
    TimedCallSystem( 'rm -rf ' . wdir );
    TimedCallSystem( 'bin/alfsim ' . paramFile );

    # initialize variables for evalutation
    lsTrees := pTrees := []:

    # get list of fasta files
    tmpFiles := TimedCallSystem( 'ls ' . wdir . '/' . mname . '/MSA/*.fa'):
    files := []:
    if tmpFiles[1] = 0 then
        files := [seq(i[1..-2], i=SplitLines(tmpFiles[2]))]:
    fi:

    # read the true alf tree
    ReadProgram( wdir . '/' . mname . '/RealTree.drw'):

    # for every sequence, reconstruct its LS and parsimony tree.
    # adapt the leaf labels to correspond to the true alf tree.
    for fastafile in files do
        msa := ReadFasta(fastafile):
        dt := PhylogeneticTree( msa[1], msa[2], DISTANCE):
        for leaf in Leaves(dt) do leaf[1] := SearchDelim('/', leaf[1])[1] od:
        lsTrees := append(lsTrees, dt):
        pt := PhylogeneticTree( msa[1], msa[2], PARSIMONY);
        for leaf in Leaves(pt) do leaf[1] := SearchDelim('/', leaf[1])[1] od:
        pTrees := append(pTrees, pt):
    od:

    return( [ invSites,
        sum( If(z=0,1,0), z=RobinsonFoulds([RealTree, op(lsTrees)])[1,2..-1] ),
        sum( If(z=0,1,0), z=RobinsonFoulds([RealTree, op(pTrees)])[1,2..-1] )
    ]);
end:

DrawResult := proc(d:matrix(numeric), legend:list(string) )
    colors := [[1,0,0],[0,0,1]]:
    cmds := []:
    for j from 2 to length(d[1]) do
        cmds := append(cmds,
           seq( LINE(d[i,1],d[i,j],d[i+1,1],d[i+1,j], 'color'=colors[j-1]), i=1..length(d)-1));
        cmds := append( cmds, LTEXT(0,40-5*j, legend[j-1],'color'=colors[j-1]) ):
    od:
    DrawPlot(cmds,axis):
end:

# starting point of program:
# read parameterfile into a string
params := ReadRawFile('synEvol_parameters.drw'):

# iterate over several values for invariante sites parameters an
# evaluate performance of reconstructed trees.
performance := []:
for invSites from .1 to .9 by .2 do
    performance := append(performance, RunSim( params, invSites)):
od:

# the final result:
print(performance):

DrawResult(performance,['DISTANCE','PARSIMONY']):
View():

After running the program, you will obtain a result similar to this one.

We can see that the parsimony tree reconstruction method seems to be more robust against this kind of model violation. But this is not generally valid!