Using Python to run simulations calling ASReml

I was having trouble to run some simulations using asreml-r when, pressed on time, I decided to use a combination of Python and plain-vanilla ASReml. This page documents a small subset of what I was doing, replacing random realizations of variance components, predicting breeding values and storing them for later analysis using R (or whatever one may want to use for the last step).

I used pretty much the standard Python distribution with the addition of numpy (a separate install) because I was playing with the idea of doing everything in Python. Random is used for the generation of normally distributed random numbers and subprocess to make calls to the operating system (in this case calling ASReml).

import numpy, random, subprocess

fname = '/Users/lap44/Dropbox/Pitfalls/simu.as'
sname = '/Users/lap44/Dropbox/simu.sln'
nsim = 1000
nbv = 3736

# Matrix to hold data
bvmat = numpy.zeros((nsim, nbv))

template = u'''
Analysis of two trials for density
 TrialNO  2
 TrialID  2 !A
 Nws      1000
 Estab    
 Plot_stp  !I
 CP Set_id !I
 Age 
 REP !I
 SETS !A
 TREE !I
 PLOT !I
 N_PLOT 
 TreeID !P 
 FCLN !I
 MCLN !I
 FAM  !I
 DEN 
 DENSVNEW
dens09pairped.txt !SKIP 1 !ALPHA !MAKE
dens09pairdat.txt !SKIP 1 !DOPART $A !BLUP 1

!part 10
DEN ~ mu TrialNO !r TreeID addvarT
2 1 0
425 0 IDEN !S2=re1varT
2562 0 IDEN !S2=re2varT
'''

def update_template(varcomps, template):
    for key in varcomps.keys():
        template = template.replace(key, str(varcomps[key]))

    return template

def save_asreml_code(filename, text):
    f = open(filename, 'w')
    f.write(text)
    f.close()
    return

def read_bv(filename, term):
    f = open(filename, 'r')
    bv = []
    for line in f:
        record = line.split()
        if record[0] == term:
            bv = bv + [float(record[2])]

    return bv

if __name__ == '__main__':

    for sim in range(nsim):
        varcomps = {u'addvarT' : random.gauss(330.2101, 45.33277),
                    u're1varT' : random.gauss(509.8517, 65.28963),
                    u're2varT' : random.gauss(169.0345, 24.93024)}

        asreml_code = update_template(varcomps, template)
        save_asreml_code(fname, asreml_code)

        # Runs asreml, reads and stores results
        subprocess.check_call('asreml ' + fname + ' 10', shell = True)
        bvs = numpy.array(read_bv(sname, 'TreeID'))
        bvmat[sim,] = bvs

    bvmean = numpy.mean(bvmat, axis=0)

    numpy.savetxt('/Users/lap44/Dropbox/Pitfalls/bvmean.txt', bvmean, fmt='%.6e')

There were a few more things going on, but this skeleton should give you an idea of how to do the basics.

Sketchbook

Using Python to run simulations calling ASReml