search

chapmanb / bcbb

Incubator for useful bioinformatics code, primarily in Python and R
http://bcbio.wordpress.com
603 stars 243 forks source link

globalenv versus globalEnv in count_diffexp (rpy2 call) #100

Closed yarden closed 9 years ago

yarden commented 9 years ago

Hi Brad,

In your nice count_diffexp.py script (https://github.com/chapmanb/bcbb/blob/5a773e3212bc7d20776cdecc12ad75c64ab2a694/stats/count_diffexp.py), you use rpy2 to call edgeR. The script references robjects.globalEnv but it looks like at least in my rpy2 version it's robjects.globalenv -- do you know what explains the difference?

>> robjects.globalenv
<Environment - Python:0x1037fcb00 / R:0x1008e8340>
>> robjects.globalEnv
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-8-c083f4c27be8> in <module>()
----> 1 robjects.globalEnv
>> rpy2.__version__
'2.5.6'

Thanks, Yarden

yarden commented 9 years ago

Nevermind, I realized this was written for a very old edgeR version (and possibly old rpy2). I rewrote a wrapper to edgeR with rpy2 that looks something like this:

    robjects.r('''
        library(edgeR)
    ''')
    base = importr("base")
    # Assume two groups (pairwise comparison)
    groups = robjects.r('''1:2''')
    read_delim = robjects.r["read.delim"]
    counts_file_params = {"sep": delimiter,
                          "row.names": gene_id_col}
    counts = read_delim(counts_fname, **counts_file_params)
    print "counts: ", counts
    # can include lib.sizes here if needed
    params = {'group' : groups}
    y = robjects.r.DGEList(counts, **params)
    y = robjects.r.calcNormFactors(y)
    y = robjects.r.estimateCommonDisp(y)
    et = robjects.r.exactTest(y, dispersion=dispersion)
    tags = robjects.r.topTags(et)
    tags_df = tags[0]
    result = {"exactTest": et,
              "tags": tags_df,
              "y": y}
    return result
chapmanb commented 9 years ago

Yarden -- sorry I missed this the other day. You're right, that code is crazy old and certainly not up to date with the best thing to do now. Thanks so much for the updated code for anyone who's looking at this. If you want to submit a pull request I'm happy to pull in a new script to not confuse anyone. Sorry again about the old code and thanks again.