Learning to use RPy2 (3)

So, I'd been thinking that I'd do the Bioconductor example in Python. The first group of commands in R is:

library('ALL')
data('ALL')
library('limma')

In Python we do:

>>> import rpy2.robjects as robjects
>>> from rpy2.robjects.packages import importr
>>> ALL = importr('ALL')
>>> limma = importr('limma')

and to get the data, we can (cheat a bit and) do

>>> robjects.r('data("ALL")')
<StrVector - Python:0x10eb15a28 / R:0x7f8371183f38>
['ALL']
>>> data = robjects.globalenv['ALL']

The second group of R commands is:

eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')]
f <- factor(as.character(eset$mol.biol))
design <- model.matrix(~f)
fit <- eBayes(lmFit(eset,design))
sel <- p.adjust(fit$p.value[,2]) < 0.05
esetSel <- eset[sel,]

I figured out how to take apart the data object above, and implement the first line. That was enough to convince me that it's not worth it.

I realized now that the use case is to get the underlying data into and back out of Python. That is, eventually we'll want to learn about how to construct an ExpressionSet object with data starting in Python data structures. But first:

>>> featureNames = robjects.r['featureNames']
>>> fn = featureNames(data)
>>> len(fn)
12625
>>> fn[:5]
<StrVector - Python:0x10eb19ef0 / R:0x7f8371eba940>
['1000_at', '1001_at', '1002_f..., '1003_s..., '1004_at']
>>> fn.rx(1)
<StrVector - Python:0x10eb19440 / R:0x7f8371f85208>
['1000_at']
>>> tuple(fn[:4])
('1000_at', '1001_at', '1002_f_at', '1003_s_at')

That's what we'll want.

The expression data:

>>> exprs = robjects.r['exprs']
>>> ed = exprs(data)
>>> ed
<Matrix - Python:0x10eb10710 / R:0x1137e6000>
[7.597323, 5.046194, 3.900466, ..., 3.095670, 3.342961, 3.842535]
>>> dim = robjects.r['dim']
>>> dim(ed)
<IntVector - Python:0x10eb19cb0 / R:0x7f8372023548>
[ 12625, 128]
>>> t = tuple(ed)
>>> t[:2]
(7.597322981163869, 5.046194285620063)
>>> len(t)
1616000
>>> import numpy as np
>>> A = np.array(list(t))
>>> A.shape = tuple(dim(ed))
>>> A.shape
(12625, 128)
>>> A[:5,:3]
array([[ 7.59732298, 5.04619429, 3.90046642],
[ 5.20211564, 5.81351983, 5.81160593],
[ 6.14364259, 7.1470674 , 6.77438247],
[ 6.76609498, 4.61191719, 4.91830898],
[ 7.12201406, 7.61231255, 6.82803053]])

The phenotypic data including mol.biol:

>>> phenoData = robjects.r['phenoData']
>>> pd = phenoData(data)
>>> mb = pd.rx2('mol.biol')
Traceback (most recent call last):
File "", line 1, in
AttributeError: 'RS4' object has no attribute 'rx2'

UPDATE: this way

>>> slotNames = r['slotNames']
>>> phenoData = r['phenoData']
>>> pd = phenoData(ALL)
>>> slotNames(pd)
<StrVector - Python:0x1023f77a0 / R:0x103140528>
['varMetada..., 'data', 'dimLabels', '.__classV...]
>>>
>>> df = pd.do_slot("data")
>>> df.rx2('sex')
<FactorVector - Python:0x1023e2170 / R:0x10268da40>
[ 2, 2, 1, ..., 2, 2, NA_integer_]
>>> mb = df.rx2('mol.biol')
>>> levels(mb)
<StrVector - Python:0x1023e2ea8 / R:0x103742d88>
['ALL1..., 'BCR/..., 'E2A/..., 'NEG', 'NUP-..., 'p15/...]
>>>

pData works:

>>> pData = robjects.r['pData']
>>> pd = pData(data)
>>> mb = pd.rx2('mol.biol')
>>> mb
<FactorVector - Python:0x1154c5cf8 / R:0x7f8374e72010>
[ 2, 4, 2, ..., 4, 4, 4]
>>> len(mb)
128
>>> mb.levels
<StrVector - Python:0x1154c5f38 / R:0x7f8371a5c660>
['ALL1..., 'BCR/..., 'E2A/..., 'NEG', 'NUP-..., 'p15/...]
>>> L = list(mb)
>>> L[:5]
[2, 4, 2, 1, 4]
>

Notice how the factor changed into ints (see the Rpy2 docs).

Now, to select the entries we want:

>>> i = mb.levels.index('ALL1/AF4') + 1
>>> j = mb.levels.index('BCR/ABL') + 1
>>> L = [True if e in [i,j] else False for e in L]
>>> sum(L)
47


>>> from rpy2.robjects import BoolVector
>>> bv = BoolVector(L)
>>> robjects.r('f <- function(e,bv) { e[,bv] }')
<SignatureTranslatedFunction - Python:0x1154cf2d8 / R:0x7f8371fa4400>
>>> f = robjects.r['f']
>>> eset = f(data,bv)
>>> dim(eset)
<IntVector - Python:0x1154cf638 / R:0x7f8375052c68>
[ 12625, 47]
>>> eset.rclass
<rpy2.rinterface.SexpVector - Python:0x10eb17ee8 / R:0x7f83750f5698>

And that's what we set out to do.

Constructing an R object from data in Python will have to wait.

Another look at the selected data:

>>> t = tuple(exprs(eset))
>>> A = np.array(list(t))
>>> A.shape = tuple(dim(eset))
>>> A.shape
(12625, 47)
>>> B = A[:100,]

Might as well plot something while we're here:

>>> fig = plt.figure()
>>> plt.hot()
>>> plt.pcolormesh(B)
<matplotlib.collections.QuadMesh object at 0x10d727610>
>>> plt.colorbar()
<matplotlib.colorbar.Colorbar instance at 0x10d707560>
>>> plt.savefig('example.png')

That is what is at the top of the post.

Learning to use RPy2 (2)

We're working on RPy2. Before I continue with Bioconductor, we should work through more of the examples in the docs.

I rashly asked a question on Stack Overflow this morning (after being stumped for several hours), then stumbled across the answer. The question is, how to load the data from ALL. In R we would just do:

data('ALL')

and then the objects would be available to us through their names. It turns out that in Python it is almost that simple. I had actually seen the answer, which is a kind of "shortcut" in which we just execute R code directly, and then grab a variable through the name assigned in R. (I used this device several times in the previous post).

(And I should note that there is apparently a more elegant solution in the bioconductor extensions to rpy2.

>>> import rpy2.robjects as robjects
>>> from rpy2.robjects.packages import importr
>>> base = importr('base')
>>> ALL = importr('ALL')
>>> data = robjects.r('data(ALL)')
>>> data.rclass
<rpy2.rinterface.SexpVector - Python:0x258b190 / R:0xdf86c8>
>>> data = robjects.globalenv['ALL']
>>> data
<RS4 - Python:0x2591490 / R:0x29a2134>
>>> data.rclass
<rpy2.rinterface.SexpVector - Python:0x258b3b0 / R:0xdf85c8>
>>> featureNames = robjects.r('featureNames')
>>> featureNames(data)
<StrVector - Python:0x23e4f08 / R:0x304fc00>
['1000..., '1001..., '1002..., ..., 'AFFX..., 'AFFX..., 'AFFX...]
>>> exprs = robjects.r['exprs']
>>> e = exprs(data)
>>> e
<Matrix - Python:0x23b2da0 / R:0x84d8000>
[7.597323, 5.046194, 3.900466, ..., 3.095670, 3.342961, 3.842535]
>>>

Anyway, we want not to make any more silly mistakes. Another elementary thing I happened upon in the docs (in Vectors) is

Access to R-style extracting/subsetting is granted though the two delegators rx and rx2, representing the R functions [ and [[ respectively.

It seems to return a strange little beast:

>>> from rpy2.robjects import r
>>> m = r.matrix(range(4),nrow=2)
>>> value = m.rx(3)
>>> value
<ListVector - Python:0x1023db5a8 / R:0x102cf9c28>
[IntVector]
<no name>:
<IntVector - Python:0x1023db368 / R:0x102cf99e8>
[ 2]
>>> type(value )
<class 'rpy2.robjects.vectors.ListVector'>
>>> value[0]
<IntVector - Python:0x1023db518 / R:0x102cf99e8>
[ 2]
>>> value[0][0]
2

We copy the code to make an R list, directly out of the docs (see ?list):

>>> r('''e1 <- new.env(); e1$a <- 10; e1$b <- 20; L = as.list(e1)''')
<ListVector - Python:0x1023da4d0 / R:0x103306140>
[FloatVector, FloatVector]
b:
<FloatVector - Python:0x1023da680 / R:0x102d1bd58>
[20.000000]
a: <class 'rpy2.robjects.vectors.FloatVector'>
<FloatVector - Python:0x1023da3b0 / R:0x102d1bd88>
[10.000000]
>>> L = r['L']
>>> L.rx2('a')
<FloatVector - Python:0x1023da998 / R:0x102d1bd88>
[10.000000]


---

Here is a somewhat fancier example involving a linear model (from the Intro). Look in the docs for the R code. In Python:

>>> from rpy2.robjects import FloatVector
>>> from rpy2.robjects.packages import importr
>>> from rpy2.robjects import globalenv as g
>>>
>>> stats = importr('stats')
>>> base = importr('base')
>>>
>>> ctl = FloatVector([4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14])
>>> trt = FloatVector([4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69])
>>> weight = ctl + trt
>>>
>>> group = base.gl(2, 10, 20, labels = ["Ctl","Trt"])
>>> g["weight"] = weight
>>> g["group"] = group
>>> lm_D9 = stats.lm("weight ~ group")
>>> print(stats.anova(lm_D9))
Analysis of Variance Table

Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
group 1 0.6882 0.68820 1.4191 0.249
Residuals 18 8.7292 0.48496
>>>
>>> print(lm_D9.names)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "contrasts" "xlevels" "call" "terms"
[13] "model"

>>> print names(lm_D9)
Traceback (most recent call last):
File "<stdin>", line 1, in
NameError: name 'names' is not defined
>>>


The next example in this section fails with two errors:

ImportError: cannot import name NA_real.

If we just use zeroes, it fails with

in __setitem__
globalenv_ri)
TypeError: All keywords must be strings (or None).

I'm just going to leave it for now.

---

PCA example

>>> from rpy2.robjects.packages import importr
>>> from rpy2.robjects import globalenv as g
>>> from rpy2.robjects import r
>>>
>>> base = importr('base')
>>> stats = importr('stats')
>>> graphics = importr('graphics')
>>>
>>> # each of these x.y could be replaced by `r`
... m = base.matrix(r.rnorm(100), ncol=5)
>>> pca = stats.princomp(m)
>>> graphics.plot(pca, main="Eigen values")
rpy2.rinterface.NULL

The graphic is at the top of the post. I didn't both to print this properly. See the first post on the docs here. But the idea is that although these functions come from the imported modules they are present in the global R namespace (if that's what they call it). In general, it's better to document your code by using the qualified name. (I know, I know, I used r here extensively! Do as I say, not as I do).

---

kwargs example

>>> from rpy2.robjects import r
>>> r('''f <- function(x='foo',y='blue') {
... if (x=='foo') print(x)
... else print (y) }''')
<SignatureTranslatedFunction - Python:0x1023d76c8 / R:0x100ef6b58>
>>> f = r['f']
>>> kwargs = {'x':"bar", 'y':"red"}
>>> f(**kwargs)
[1] "red"
<StrVector - Python:0x1023db170 / R:0x102d0fc28>
['red']

I guess that's enough for one post.

Random Quotes, get your Random Quotes

A nice collection of quotes here. Perhaps my favorite:

"A month in the laboratory can often save an hour in the library."

-- F. H. Westheimer

Advice to a neuroscientist.

Money quote for folks at my place:

Third: learn how to do your own data analysis. Know statistics well. Know at least some basic programming/scripting in Python, R, Matlab, etc. This will be of immense value in helping you get your research done efficiently and correctly, without needing to rely on other people's code (and time and commitment). This will become more important as our field becomes more data driven.

--Bradley Voytek (here)
h/t Tyler Cowen

My prior on the probability that they'll notice? 0.01 :)

Learning to use the RPy2 module (1)

Before we tackle turning the Bioconductor example into Python code, I need to review some basic usage of RPy2 as given in the docs. I posted a bit about this before, but my exploration wasn't systematic, and also the posts contain issues which I solved and described in later posts.

This exercises the first half of what's on the linked page. Python code followed output in bold.

import rpy2.robjects as robjects
from rpy2.robjects.packages import importr

r = robjects.r
g = robjects.globalenv
print type(r)

<class 'rpy2.robjects.R'>



# from the R base package
# any other name comes from .globalEnv
pi = r['pi']
print pi[0]
print type(pi)

3.14159265359
<class 'rpy2.robjects.vectors.FloatVector'>



# r is callable with code to be evaluated
piplus = r('piplus = pi + 1')
print type(piplus)
print piplus
print piplus[0]

<class 'rpy2.robjects.vectors.FloatVector'>
[1] 4.141593

4.14159265359
[1] 4.141593



# explicitly get from .globalEnv
piplus_from_g = g['piplus']
print piplus_from_g
print piplus == piplus_from_g
print piplus[0] == piplus_from_g[0]

False
True



# define an R function and call it
r('''
f <- function(start=1,stop = 5) {
n = 0
for (j in start:stop) { n = n + j }
print (n)
}
''')

# it's available in .globalEnv
f = g['f']
f()


[1] 15



# but also from the running R process
f = r['f']
f(4)

[1] 9



# interpolating an R object into code
letters = robjects.r['letters']
s = letters[1:6].r_repr()
rcode = 'paste(%s, collapse="-")' %(s)
res = robjects.r(rcode)
print(res)

[1] "b-c-d-e-f"



# more on calling R functions
rsum = r['sum']
print rsum(robjects.IntVector([1,2,3]))[0]

6



# with keyword
rsort = r['sort']
iv = robjects.IntVector([3,1,2])
res = rsort(iv, decreasing=True)
print res.r_repr()

c(3L, 2L, 1L)



# plotting
gd = importr('grDevices')
ofn = '/Users/telliott_admin/Desktop/plot.pdf'
gd.pdf(ofn)

x = robjects.IntVector(range(10))
y = r.rnorm(10)
r.layout(r.matrix(robjects.IntVector([1,2,3,2]), nrow=2, ncol=2))
r.plot(r.runif(10), y, xlab="runif", ylab="foo/bar", col="red")
gd.dev_off()



Finally, there are R's special operators like %in% and %*%. I didn't find this yet, so try wrapping it:

r('''
is_in <- function(value, container) {
value %in% container
}
''')

is_in = g['is_in']
L = robjects.StrVector(list('abcde'))
print is_in('a',L)
print is_in('f',L)

[1] TRUE

[1] FALSE





r('''m_mult <- function(m1,m0) { m1 %*% m0 }''')
m_mult = g['m_mult']
r('''mtoi <- function(m) { as.integer(m) }''')
mtoi = g['mtoi']

# an old friend
m0 = r.matrix(robjects.IntVector([0,1,1,1]), nrow=2)
m = m0
for i in range(10):
print i + 3,
m = m_mult(m,m0)
print mtoi(m)

3 [1] 1 1 1 2

4 [1] 1 2 2 3

5 [1] 2 3 3 5

6 [1] 3 5 5 8

7 [1] 5 8 8 13

8 [1] 8 13 13 21

9 [1] 13 21 21 34

10 [1] 21 34 34 55

11 [1] 34 55 55 89

12 [1] 55 89 89 144



Explanation of the example code

This is a brief post to explain the R (Bioconductor) code that we used to draw a heat map for gene expression data (here and here). Again, first we load the example data:

library('ALL')
data('ALL')
library('limma')

The second group of commands is:

eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')]
f <- factor(as.character(eset$mol.biol))
design <- model.matrix(~f)
fit <- eBayes(lmFit(eset,design))
sel <- p.adjust(fit$p.value[,2]) < 0.05
esetSel <- eset[sel,]

The first of these constructs a "selector." The values in ALL$mol.biol are factors; two of the possibilities are in the vector c('BCR/ABL','ALL1/AF4'). The result is a vector of booleans. To give a simpler example:

> vowels = c('a','e','i','o','u')
> letters[1:6]
[1] "a" "b" "c" "d" "e" "f"
> letters %in% vowels
[1] TRUE FALSE FALSE FALSE TRUE FALSE FALSE
[8] FALSE TRUE FALSE FALSE FALSE FALSE FALSE
[15] TRUE FALSE FALSE FALSE FALSE FALSE TRUE
[22] FALSE FALSE FALSE FALSE FALSE

The vector of True/False values from ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4' is used to select the appropriate columns from the data, but is acting on the complex ExpressionSet object, as explained last time. The result is reduced in size.

> dim(eset)
Features Samples
12625 47

The second line is just a bit of coaxing. Because the mol.biol descriptions for the smaller eset were derived from a "factor" with 6 possibilities, we remove the other 4 before moving on to model construction:

> f = eset$mol.biol
> class(f)
[1] "factor"
> length(f)
[1] 47
> levels(f)
[1] "ALL1/AF4" "BCR/ABL" "E2A/PBX1" "NEG"
[5] "NUP-98" "p15/p16"
> f = factor(as.character(f))
> class(f)
[1] "factor"
> length(f)
[1] 47
> levels(f)
[1] "ALL1/AF4" "BCR/ABL"


---

The middle two lines is where the magic happens. I'm going to skip that for the moment. You can read about it in (Gentleman et al, PMID 15461798). The result is in the variable `fit`. One of its sub-objects is a list of p-values:

> class(fit$p.value)
[1] "matrix"
> fit$p.value[1:3,]
(Intercept) fBCR/ABL
1000_at 2.713661e-58 0.05437095
1001_at 3.001944e-43 0.24423133
1002_f_at 2.273600e-46 0.10091357
> colnames(fit$p.value)
[1] "(Intercept)" "fBCR/ABL"
> pv = fit$p.value[,2]
> class(pv)
[1] "numeric"
> pv[1:4]
1000_at 1001_at 1002_f_at 1003_s_at
0.05437095 0.24423133 0.10091357 0.30303775

p.adjust is a correction for the large number of comparisons. You can read about it by doing ?p.adjust. I wasn't clear on which is the default method, but it's not Bonferroni.

> r1 = p.adjust(pv)
> r2 = p.adjust(pv,method='bonferroni')
> all(r1 == r2)
[1] FALSE

So we just construct a selector for those rows with p < 0.05, and again, whittle down the ExpressionSet to something manageable:

> sel <- p.adjust(pv) < 0.05
> esetSel <- eset[sel,]
> dim(esetSel)
Features Samples
165 47


---

In the third part, we draw the heatmap.

color.map <- function(s) {
if (s=='ALL1/AF4') 'red' else 'blue'}
patient.colors <- unlist(lapply(esetSel$mol.biol,color.map))
heatmap(exprs(esetSel),
col=topo.colors(100),
ColSideColors=patient.colors)

In the first two lines we construct a function called color.map, and then use it to construct a vector of patient.colors. This is for the colored bars at the top of the plot.

> color.map <- function(s) {
+ if (s=='ALL1/AF4') 'red' else 'blue'}
>

We use lapply to feed values to our function, but the result is a list. So we use unlist to extract what we want.

> values = esetSel$mol.biol
> values[1:4]
[1] BCR/ABL BCR/ABL ALL1/AF4 BCR/ABL
6 Levels: ALL1/AF4 BCR/ABL E2A/PBX1 ... p15/p16
> colors = lapply(values,color.map)
> class(colors)
[1] "list"
> colors[1:2]
[[1]]
[1] "blue"

[[2]]
[1] "blue"

> patient.colors = unlist(colors)
> patient.colors[1:4]
[1] "blue" "blue" "red" "blue"
> class(patient.colors)
[1] "character"

We draw the heatmap with the command:

heatmap(exprs(esetSel)

All the rest is eye candy. In particular, we use the topo.colors. Without that, we get the default.



With the topo.colors, we would get colors as in the graphic at the top of the post.

heatmap(exprs(esetSel),col=topo.colors(100))

The last argument to heatmap gives us the bars on the columns to identify the sample by translocation.

We've explained everything but the actual statistical model fitting. That's for another time.

Dissecting an ExpressionSet object

Last time, we were in R using some Bioconductor code. If we do this:

library('ALL')
data('ALL')

then the variable ALL is available

> ALL
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 128 samples
element names: exprs
protocolData: none
phenoData
sampleNames: 01005 01010 ... LAL4 (128 total)
varLabels: cod diagnosis ... date last seen (21 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 14684422 16243790
Annotation: hgu95av2

(note: your output may vary slightly depending on the R and Bioconductor version.)

Obviously the object that ALL refers to is complicated. Its components are called "slots" and they are accessible with the syntax: object_name@slot_name. For example:

pData = ALL@phenoData

You can find out which slots are available by parsing the above output or just doing:

> slotNames(ALL)
[1] "assayData" "phenoData"
[3] "featureData" "experimentData"
[5] "annotation" "protocolData"
[7] ".__classVersion__"


---

In turn the slots are also complex objects. Technically, a phenoData object is:

> class(ALL@phenoData)
[1] "AnnotatedDataFrame"
attr(,"package")
[1] "Biobase"

More important, its sub-objects are available with the syntax slot_name$sub-object_name. But first we must discover which names are available:

> pData = ALL@phenoData
> varLabels(pData)
[1] "cod" "diagnosis" "sex"
[4] "age" "BT" "remission"
[7] "CR" "date.cr" "t(4;11)"
[10] "t(9;22)" "cyto.normal" "citog"
[13] "mol.biol" "fusion protein" "mdr"
[16] "kinet" "ccr" "relapse"
[19] "transplant" "f.u" "date last seen"

So, for example:

> pData$sex[1:5]
[1] M M F M M
Levels: F M

> pData$mol.biol[1:5]
[1] BCR/ABL NEG BCR/ABL ALL1/AF4 NEG
Levels: ALL1/AF4 BCR/ABL E2A/PBX1 NEG NUP-98 p15/p16

> levels(pData$mol.biol)
[1] "ALL1/AF4" "BCR/ABL" "E2A/PBX1" "NEG"
[5] "NUP-98" "p15/p16"

At this level, the objects are reasonably simple.

> class(pData$mol.biol)
[1] "factor"
> class(pData$sex)
[1] "factor"


---

Don't forget the other method names:

> sampleNames(pData)[1:5]
[1] "01005" "01010" "03002" "04006" "04007"

> featureNames(ALL)[1:8]
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at"
[5] "1004_at" "1005_at" "1006_at" "1007_s_at"

> head(varMetadata(pData))
labelDescription
cod Patient ID
diagnosis Date of diagnosis
sex Gender of the patient
age Age of the patient at entry
BT does the patient have B-cell or T-cell ALL
remission Complete remission(CR), refractory(REF) or NA. Derived from CR

And, I should point out that (without actually knowing it at the time) I shadowed a real function in the namespace by assigning to pData above.

If we reboot R and do:

library('ALL')
data('ALL')

we can use a different method for access:

> p = pData(ALL)
> class(p)
[1] "data.frame"
> names(p)
[1] "cod" "diagnosis" "sex"
[4] "age" "BT" "remission"
[7] "CR" "date.cr" "t(4;11)"
[10] "t(9;22)" "cyto.normal" "citog"
[13] "mol.biol" "fusion protein" "mdr"
[16] "kinet" "ccr" "relapse"
[19] "transplant" "f.u" "date last seen"
> p$sex[1:5]
[1] M M F M M
Levels: F M


---

The very first slot holds the gene expression data.

> aData = ALL@assayData
> class(aData)
[1] "environment"

I'm not really sure what that class is. And, unfortunately, it is not easy to find which sub-object names are available for this slot:

> names(aData)
NULL
> elements(aData)
Error: could not find function "elements"
> elementNames(aData)
Error: could not find function "elementNames"

except we can probably guess that exprs is one (based on the output of the call `ALL` above:

> e = aData$exprs
> class(e)
[1] "matrix"
> e[1:3,1:3]
01005 01010 03002
1000_at 7.597323 7.479445 7.567593
1001_at 5.046194 4.932537 4.799294
1002_f_at 3.900466 4.208155 3.886169
> dim(e)
[1] 12625 128
> rownames(e)[1:8]
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at"
[5] "1004_at" "1005_at" "1006_at" "1007_s_at"
> colnames(e)[1:3]
[1] "01005" "01010" "03002"

-------------------------------------------------------
So far, so good.

Now, obviously we can select rows and columns of a matrix like ALL@assayData$exprs, which we've assigned to the variable e:

> e[1:2,1:3]
01005 01010 03002
1000_at 7.597323 7.479445 7.567593
1001_at 5.046194 4.932537 4.799294

What the Bioconductor guys have done is make this assignment valid for the ExpressionSet object itself:

> ALL[1:2,1:3]
ExpressionSet (storageMode: lockedEnvironment)
assayData: 2 features, 3 samples
element names: exprs
protocolData: none
phenoData
sampleNames: 01005 01010 03002
varLabels: cod diagnosis ... date last seen (21 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 14684422 16243790
Annotation: hgu95av2

That's either very cool or very confusing depending on where you sit.

---

Finally, the annotation. This requires another package:

ann = ALL@annotation
> ann
[1] "hgu95av2"

library('annaffy')
probeids <- featureNames(ALL)
symbols <- aafSymbol(probeids, 'hgu95av2')

> symbols <- aafSymbol(probeids, 'hgu95av2')
[1] "You are missing hgu95av2 looking to see if it is available."
Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu/bioc/bin/macosx/leopard/contrib/2.11
Package hgu95av2 is available for download, would you like to install? [y/n]
y
Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu/bioc/bin/macosx/leopard/contrib/2.11
trying URL 'http://bioconductor.org/packages/2.6/data/annotation/bin/macosx/leopard/contrib/2.11/hgu95av2_2.2.0.tgz'
Content type 'application/x-gzip' length 9340297 bytes (8.9 Mb)
opened URL
==================================================
downloaded 8.9 Mb


The downloaded packages are in
/var/folders/3g/3gHxNStTGZaWkGngOL39-++++TI/-Tmp-//RtmptUUJr1/downloaded_packages
Loading required package: hgu95av2

******* Deprecation warning *******:

The package 'hgu95av2' is deprecated and will not be supported in the future.

Instead we strongly reccomend that you should start using the 'hgu95av2.db' package.

We hope you will enjoy these new packages. If you still have questions, you can search
(http://dir.gmane.org/gmane.science.biology.informatics.conductor) or ask the mailing list
(http://bioconductor.org/docs/mailList.html).

To repeat:

library('annaffy')
probeids <- featureNames(ALL)
symbols <- aafSymbol(probeids, 'hgu95av2')
> class(symbols)
[1] "aafList"
> genenames = getText(symbols)
> genenames[1:3]
[1] "MAPK3" "TIE1" "CYP2C19"
> probeids[1:3]
[1] "1000_at" "1001_at" "1002_f_at"

We'll write this all down on a couple post-it notes and put it on the monitor for reference.

A first heat map using Bioconductor

I don't recall what the exact Google search target was that turned up a post about making a heatmap of gene expression data, but it really changed my life. It lead me to these two pages describing how to do it in R (here) and Python (here), and that motivated me to get into R, as much as I've been able to in spite of all I don't like about it.

I thought I'd work on this a little bit again; it's been 3 years or so since I gave a couple of talks about this, so I'm rusty, but it'll be fun. I would recommend that you go through the original posts by Peter Cock first. I'll try to add something of my own.

Of course, it's all based on the Bioconductor project in R, and there is a very nice write-up in the original paper (Gentleman et al, PMID 15461798). I even bought the book (our library sucks), and it's pretty good too, especially if you are serious about the field.

To begin with, to motivate things, we will recreate the graphic. We need a couple of R packages from Bioconductor to do this. You could get them from the Package Installer, but I was having trouble and so I tried this:

source("http://bioconductor.org/biocLite.R")
biocLite('Biobase')
biocLite('ALL')
biocLite('limma')
biocLite('annaffy')

Note: all these lines would actually look like this in R:

> source("http://bioconductor.org/biocLite.R")

For most of the commands here, I've left out the prompt so you can more easily copy and paste them. And I'm not showing the info that's printed as the packages are downloaded and installed.

We're going to just quickly do the plots, then later I'll explain more about it. First we need to load two modules and some data:

library('ALL')
data('ALL')
library('limma')

The first command will output:

Loading required package: Biobase

Welcome to Bioconductor

Vignettes contain introductory material. To view, type
'openVignette()'. To cite Bioconductor, see
'citation("Biobase")' and for packages 'citation(pkgname)'.

Next we do the analysis (I know this won't make much sense yet, but we're just making sure we get the graphic without an error):

eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')]
f <- factor(as.character(eset$mol.biol))
design <- model.matrix(~f)
fit <- eBayes(lmFit(eset,design))
sel <- p.adjust(fit$p.value[,2]) < 0.05
esetSel <- eset[sel,]

And in the third part we plot the data:

color.map <- function(s) {
if (s=='ALL1/AF4') 'red' else 'blue'}
patient.colors <- unlist(lapply(esetSel$mol.biol,color.map))
heatmap(exprs(esetSel),
col=topo.colors(100),
ColSideColors=patient.colors)

The graphic is at the top of the post.

An overview of what we've done. We filtered for all samples with one of two chromosomal translocations. In the plots, the columns from the 'ALL1/AF4' samples are red, those with 'BCR/ABL' are blue. We now have 47 samples.

> dim(eset)
Features Samples
12625 47

Each sample has 12625 "features"---gene expression measurements. We've done some fancy statistics and then filtered for those with p < 0.05 (correcting for the large number of measurements). In the end, we have 165 genes of interest to plot.

> dim(esetSel)
Features Samples
165 47

By default, heatmap does clustering and draws dendrograms, but you can turn this off. Let's turn off the column clustering only, and since we haven't substituted the real gene names, let's turn that off too:

heatmap(exprs(esetSel),
col=topo.colors(100),
Colv=NA,
ColSideColors=patient.colors,
labRow=NA)




More later on this. Our goal is to understand more about "ExpressionSet" objects and the kinds of manipulation that are possible, the second is to see how much of this we can do from Python (using rpy2), so that the code is transparent and makes intuitive sense.

Zenoss

I found this question on SF about network monitoring for OS X. One of answers described an open source solution called Zenoss.

From the looks of things it's very powerful. Here is a slightly old post about using it on OS X.

I downloaded and installed Zenoss. Unfortunately, it made my system feel unresponsive, even when the App was not running. (There are a bunch of process that run at startup, and eat up CPU).



Luckily, there is an uninstaller. I couldn't figure out how to get there in the Finder, but I just did:

open -a /usr/local/zenoss/uninstall.app/

It's probably no surprise that the uninstaller took more than 5 min to run, but it seems at least to have done its job and now everything is nice and snappy again.

High standards

The nice thing about standards is that there are so many to choose from. And if you really don't like all the standards you just have to wait another year until the one arises you are looking for.

-- A. Tanenbaum, "Introduction to Computer Networks", quoted here

Trying Ubuntu Linux (6)

Continuing with my Ubuntu Linux (11.0.4) server running in a VirtualBox VM under Mac OS X Lion, and making its capabilities available to other machines on my LAN. In two previous posts I described my search for a solution (here), and the actual solution (here), which used port forwarding with VirtualBox networking in NAT mode.

There is at least one other possible solution, which is to do the networking in bridged mode. To quote the docs:

With bridged networking, VirtualBox uses a device driver on your host system that filters data from your physical network adapter. This driver is therefore called a "net filter" driver. This allows VirtualBox to intercept data from the physical network and inject data into it, effectively creating a new network interface in software. When a guest is using such a new software interface, it looks to the host system as though the guest were physically connected to the interface using a network cable: the host can send data to the guest through that interface and receive data from it. This means that you can set up routing or bridging between the guest and the rest of your network.

To begin with I did something (that may not be necessary), I turned off the NAT rule from before (in OS X):

VBoxManage modifyvm Ubuntu --natpf1 delete "server"

Then I set networking to bridged in the VirtualBox Settings:



And restarted Apache (in Ubuntu):

sudo /etc/init.d/apache2 restart

ifconfig gives me the IP address that we have obtained from the DHCP server in Airport Extreme:



If I use the AirPort Utility I can see the corresponding MAC and IP addresses.



Now if I go to another machine on the network and point Safari at that IP, it works:



[ UPDATE: In the screenshot we specified port 8080, and that works because of what we did last time. But you can also do the ipadress only and it works too. ]

Pretty code

[ UPDATE2: In an earlier version of this post, I missed one HTML tag in the code below, as a result the example wouldn't run as described when copy/pasted. Sorry. ]
I've been looking into solutions for posting code on the blog. Since the beginning, I've used <table> tags like this:
<div style="overflow-x: scroll ">
<table bgcolor="#ffffb0" border="0" width="100%" padding="4">
<tbody><tr><td><pre style=" hidden;font-family:monaco;">
my code here
</pre></table></div>


With results looking like a bit like this, but larger :)


I borrowed the technique from Peter Norvig. He's moved on, but I hadn't.

In the last few days Blogger changed something that inserts newlines into such tables. I found another guy who noticed, but no solution. (Perhaps if I used the "new" editor?).

So, I googled around for syntax highlighting. Here is a listing of a bunch of methods. The one I decided to try is google-code-prettify. Like all the methods, it uses javascript and some css to do its work. Instructions are in the readme. The footprint in the html is pretty small, and I like the minimalist style.

To test it, I downloaded the distribution (see the link), and then put a file with the following html in the same directory as the source:


<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE html>
<html ... >
.
<head ... >
.
<link href="prettify.css" type="text/css" rel="stylesheet" />
<script type="text/javascript" src="prettify.js"></script>

</head>
<body onload="prettyPrint()">

<pre class="prettyprint" lang-python>
def f(n, s='abc'):
    print 'Hello world!'
    for i in range(n):
        # this is a loop
        print s
</pre>

<pre class="prettyprint">
#include <math.h>
#include "Python.h"

int process (const char *c, double f) {
    int i = 0;
    if (c[0] == 'c') {
        i = floor(f);
    }
    return i;
}
</pre>


</body>
</html>


and then drag-n-dropped it on Safari:



One small problem is this note:

Include the script and stylesheets in your document (you will need to make sure the css and js file are on your server, and adjust the paths in the script and link tag)

Server, what server?

Seriously, I'd be happy to do it, but I don't want to pay Comcast $$ for the privilege.

However, there is a very nice answer on Stack Overflow (here) that shows how to embed the relatively simple css instructions into the html document. Then you just load the code from


http://google-code-prettify.googlecode.com/svn/trunk/src/prettify.js
http://google-code-prettify.googlecode.com/svn/trunk/src/lang-css.js


I tested it in a simple html document and it worked. Now I'll have to try to incorporate all this into a Blogger template.

[ UPDATE: I see it choked on the <math.h>. I'll have to remember that. ]

Trying Ubuntu Linux (5)

The project I had yesterday (here) was to set up my server (Apache) in Linux under VirtualBox so that it is visible from the host (OS X Lion). I solved that eventually by sticking with NAT mode and issuing this command in Terminal (with Ubuntu powered down):

VBoxManage modifyvm Ubuntu --natpf1 "server,tcp,,8080,,8080"

This is exactly what the docs say to do. The extra twist was to modify /etc/apache2/ports.conf to add:

NameVirtualHost *:8080
Listen 8080

and to modify /etc/apache2/sites-available/default. I basically duplicated the whole <VirtualHost entry and changed the second part to read:

<VirtualHost *:8080>
..

And now it works. From the host in Safari or using curl in Terminal:

http://127.0.0.1:8080/cgi-bin/test.py

SERVER_SOFTWARE Apache/2.2.17 (Ubuntu)
SCRIPT_NAME /cgi-bin/test.py
SERVER_SIGNATURE
Apache/2.2.17 (Ubuntu) Server at localhost Port 8080
REQUEST_METHOD GET
SERVER_PROTOCOL HTTP/1.1
QUERY_STRING
..
HTTP_USER_AGENT Mozilla/5.0 (Macintosh; Intel Mac OS X 10_7) AppleWebKit/534.48.3 (KHTML, like Gecko) Version/5.1 Safari/534.48.3
HTTP_CONNECTION keep-alive
SERVER_NAME localhost
..
SERVER_PORT 8080
..
SERVER_ADMIN webmaster@localhost
HTTP_HOST localhost:8080
REQUEST_URI /cgi-bin/test.py
HTTP_ACCEPT text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8
GATEWAY_INTERFACE CGI/1.1
..
HTTP_ACCEPT_LANGUAGE en-us
HTTP_ACCEPT_ENCODING gzip, deflate


The other thing I've done today is to install EMBOSS and PyCogent.

sudo apt-get install emboss

How easy is that! The rebase files are at

http://rebase.neb.com/rebase/rebase.files.html

We want # 5 and # 31 as described. chmod for the EMBOSS directories so we can write there:

sudo chmod 777 -R /usr/share/EMBOSS
rebaseextract

tell it where the files are, and it'll do its thing.

head -5 /usr/share/EMBOSS/data/REBASE/embossre.enz > cat

# REBASE enzyme patterns for EMBOSS
#
# Format:
# namepatternlenncutsbluntc1c2c3c4
#

put a sequence file DA19.txt on the Desktop

remap -sequence DA19.txt

DA19


SgeI
| PcsI
SgeI | SgeI
Cac8I | NspI | | SgeI HindIII
| PcsI Cac8I | SgeI | | | SgeI | Cac8I
\ \ \ \ \ \ \ \ \ \ \
GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGGAGCGAATGAAAGCTTGCT

..

[ UPDATE: Something wrong with the formatting here. The HindIII site, AAGCTT, is at the far right. ]

Get Cython:

sudo apt-get install cython

Get Pycogent from http://sourceforge.net/projects/pycogent/:

tar -xzvf ~/Downloads/PyCogent-1.5.1.tgz
python setup.py build
sudo python setup.py install

I tested Pocogent by downloading sequences and making a phylogenetic tree as described here (I had to install a MUSCLE binary first). I pasted the source into four scripts.

Not exactly as before, but it looks reasonable.

Trying Ubuntu Linux (4)

As discussed in previous posts (here, here, here), I have Ubuntu Linux 11.0.4 running under VirtualBox 4.1.0 on my Powerbook with the host running OS X Lion 10.7 (in a second partition, actually). I set up and tested Apache Server 2.2 on the guest. Now I'd like to access the guest server from the host machine, and later the other machines on my Wi-Fi network. I haven't yet succeeded, so yesterday I asked a question on StackOverflow, which got moved by someone to a related site, "serverfault" (here).

I got a couple of good suggestions, but didn't yet solve my problem. However, there are a bunch of "related questions" in the sidebar that look promising. Time to be (even more) systematic.

The VirtualBox manual (here) has a section: Configuring port forwarding with NAT.

So it's pretty clear that we could stay with NAT. In that case we need to do "port forwarding," although other people have suggested using "bridged mode." In the NAT method, we simply instruct the VM (as we would a router) that packets arriving on certain ports should be "forwarded to the guest, on the same or a different port." The manual says to do this from the command line in the host (the first line is their template, the second my implementation):

VBoxManage modifyvm "VM name" --natpf1 "guestssh,tcp,,2222,,22"
VBoxManage modifyvm Ubuntu --natpf1 "apache,tcp,,8888,,80"

There are six fields in the last argument: a name, which is "purely descriptive", the protocol for forwarding, and then two sets of an ip address and port. So, this example "forwards all TCP traffic arriving on the localhost interface (127.0.0.1) via port 2222 to port 22 in the guest":

VBoxManage modifyvm "VM name" --natpf1 "guestssh,tcp,127.0.0.1,2222,,22"

From the docs

Forwarding host ports < 1024 impossible:
On Unix-based hosts (e.g. Linux, Solaris, Mac OS X) it is not possible to bind to ports below 1024 from applications that are not run by root. As a result, if you try to configure such a port forwarding, the VM will refuse to start.

So that's why we used 8888. Under Port Forwarding we forward (in the second example) from the host port 8888 to the guest port 80.

I used the GUI to set this.




I believe we should not have to do the command line version. Just to be safe, I remove the forwarding rule I set when playing around, by doing this from the command line:

VBoxManage modifyvm Ubuntu --natpf1 delete "apache"

And then repeated setting up forwarding in the GUI (Settings > Network). To be sure that everybody has got the word, we:

Quit Ubuntu
Restart VirtualBox
Restart Ubuntu
Restart Apache with:

sudo /etc/init.d/apache2 restart

But it doesn't work. From the guest, I point Firefox at localhost or 127.0.0.1 and I can see the index page or run my scripts. Note that this works even though the Apache restart gave this message:

apache2: Could not reliably determine the server's fully qualified domain name, using 127.0.1.1 for ServerName

But from Safari on the host with the same ip address I get "Safari can't connect to the server" and if I add the port (127.0.0.1:8888) it just hangs. 127.0.1.1 also doesn't work.

Now, what could be wrong? From serverfault (here):

Iain says to try bridged mode.

I tested that early on but it didn't make any difference. Also, the docs clearly say that NAT should work. I'll have to test some more, but first..

Eric Fortis says:

"then from the host PC browser access the IP of the Ubuntu virtual machine, instead of 127.0.0.1. ifconfig will show you your IP."

te@VB:~$ ifconfig
eth0 Link encap:Ethernet HWaddr 08:00:27:88:33:a4
inet addr:10.0.2.15 Bcast:10.0.2.255 Mask:255.255.255.0
inet6 addr: fe80::a00:27ff:fe88:33a4/64 Scope:Link
UP BROADCAST RUNNING MULTICAST MTU:1500 Metric:1
RX packets:122 errors:0 dropped:0 overruns:0 frame:0
TX packets:182 errors:0 dropped:0 overruns:0 carrier:0
collisions:0 txqueuelen:1000
RX bytes:43296 (43.2 KB) TX bytes:23263 (23.2 KB)

lo Link encap:Local Loopback
inet addr:127.0.0.1 Mask:255.0.0.0
inet6 addr: ::1/128 Scope:Host
UP LOOPBACK RUNNING MTU:16436 Metric:1
RX packets:18 errors:0 dropped:0 overruns:0 frame:0
TX packets:18 errors:0 dropped:0 overruns:0 carrier:0
collisions:0 txqueuelen:0
RX bytes:1642 (1.6 KB) TX bytes:1642 (1.6 KB)

If I'm reading this correctly, the second entry is for "local loopback", that is, it doesn't go through the virtual network card. And the first one uses eth0. So that might make a difference (Ethernet v. Wi-Fi)

Also, the ip address is 10.0.2.15, which I think should not be visible beyond the "router"---the VM.

The third comment is from anthonysomerset:

is apache in your vhost configured to listen on the correct ports? this will work in bridged or nat mode, also check any firewall rules in the guest and that network access works out of the guest as well.

So:
- this will work in bridged or NAT
- make sure the Listen directive is correct
- check the firewall
- check that network access works out of the guest

The guest is set up to listen on port 80. That is the standard port set up in Apache. I wrote a Python script to filter comments (filter.py) and do:

python filter.py /etc/apache2/apache2.conf

Include mods-enabled/*.load
Include mods-enabled/*.conf
Include httpd.conf
Include ports.conf
..
Include conf.d/
Include sites-enabled/

So there are lots of places to look for possible conflicts! Luckily, httpd.conf is empty. ports.conf:

..
NameVirtualHost *:80
Listen 80
..

I'm not quite sure about the first one, but the Listen part is correct.

Firewall rules are here (at least, some are):

te@VB:/etc/apache2/sites-available$ python ~/Desktop/filter.py default
<VirtualHost *:80>
ServerAdmin webmaster@localhost
DocumentRoot /var/www
<Directory />
Options FollowSymLinks
AllowOverride None
</Directory>
<Directory /var/www/>
Options Indexes FollowSymLinks MultiViews
AllowOverride None
Order allow,deny
Allow from all
# Allow from 127.0.0.0/255.0.0.0 ::1/128
</Directory>
..

Skipping the script stuff, this has been edited from my previous rule restricting access to go back to Allow from all.

About checking network access from the guest: if he means accessing a server on the host machine, I haven't tried that yet. But of course, Firefox works.

So, out of all this, the only thing I can see is that I don't know what this stuff is about.

NameVirtualHost *:80
<VirtualHost *:80>

I tried switching from Wi-Fi to Ethernet in the host, but that didn't help.

I'm going to have to look through all the other posts related to this topic. But it seems like it's getting to be too hard, after all, I used VirtualBox to make things easy.

Keeping up appearances

Something's changed in blogger that messes with my code examples. The last normal one is from Sat August 6, and the first unusual one is from last Thu August 11. It looks the same as ever in Preview mode, but after publishing, all the code text has acquired an extra newline and is double-spaced.

Code is enclosed within tags:

<div style="overflow-x: scroll ">
<table bgcolor="#ffffb0" border="0"
width="100%" padding="4">
<tbody>
<tr>
<d>
<pre style=" hidden;font-family:monaco;">
..
</pre>
</table>
</div>

I liked this method because the output is monospaced, and the background colored according to the type of code (Python, R, output, etc).

I guess I'll have to look into CSS. If anyone has a good set of instructions to point me to, I'd be grateful.

Python C Extension

Following Andrew Dalke's example, this is perhaps the simplest possible C extension for Python. For explanation of what's happening, see his page.

[ Something's gone wacko with this post---it's single spaced in preview, double after posting ]

x.c

#include <math.h> #include "Python.h" int process (const char *c, double f) { int i = 0; if (c[0] == 'c') { i = floor(f); } return i; } PyDoc_STRVAR(x__doc__, "x module for process'ing stuff"); PyDoc_STRVAR(process__doc__, "c,f -> do something with f if c == 'c'"); static PyObject * py_process(PyObject *self, PyObject *args) { double f = 0; char *c; int i; if (!PyArg_ParseTuple(args, "sd:process", &c, &f)) return NULL; i = process(c,f); return PyInt_FromLong((long) i); } static PyMethodDef x_methods[] = { {"process", py_process, METH_VARARGS, process__doc__}, {NULL, NULL} /* sentinel */ }; PyMODINIT_FUNC initx(void) { Py_InitModule3("x", x_methods, x__doc__); }

setup.py

from distutils.core import setup, Extension setup(name="x", version="0.0", ext_modules = [Extension("x", ["x.c"])])

> python setup.py build running build running build_ext building 'x' extension creating build creating build/temp.macosx-10.6-universal-2.6 gcc-4.2 -fno-strict-aliasing -fno-common -dynamic -DNDEBUG -g -fwrapv -Os -Wall -Wstrict-prototypes -DENABLE_DTRACE -arch i386 -arch ppc -arch x86_64 -pipe -I/System/Library/Frameworks/Python.framework/Versions/2.6/include/python2.6 -c x.c -o build/temp.macosx-10.6-universal-2.6/x.o creating build/lib.macosx-10.6-universal-2.6 gcc-4.2 -Wl,-F. -bundle -undefined dynamic_lookup -arch i386 -arch ppc -arch x86_64 build/temp.macosx-10.6-universal-2.6/x.o -o build/lib.macosx-10.6-universal-2.6/x.so

You can control which architecture is build via setup.py with:

export ARCHFLAGS="-arch x86_64"

> ln -s build/lib.macosx-10.6-universal-2.6/x.so x.so > > python -c "import math; import x; print x.process('c',math.pi)" 3 > python -c "import math; import x; print x.process('a',math.pi)" 0 > python Python 2.6.1 (r261:67515, Jun 24 2010, 21:47:49) [GCC 4.2.1 (Apple Inc. build 5646)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> import x >>> x.__doc__ "x module for process'ing stuff" >>> x.process.__doc__ "c,f -> do something with f if c == 'c'" >>>

Looks good.

Note: on my 64-bit machines I'm getting a warning:

x.c: In function ‘process’: x.c:7: warning: implicit conversion shortens 64-bit value into a 32-bit value

Converting i to long doesn't help. I'll try to look into that.

Trying Ubuntu Linux (3)

I have an idea that someday I might try to teach students how to set up a web site (for example, EMBOSS on a local network would be nice). Python would be a perfect fit and it should be fun and extensible. I would probably want a neutral platform to do that (since most would be trained on Windows and I know nothing except OS X), so I looked into adding Apache Server to my Ubuntu 11.0.4 Desktop install (posts here and here).

There is a nice guide on the web from Ubuntu (here), and some simple instructions are here, and the Apache docs are here. I could have started by grabbing the Server version of Ubuntu, but instead I did:

sudo apt-get install apache2

As easy as that. Start and stop commands:

sudo /etc/init.d/apache2 start sudo /etc/init.d/apache2 stop sudo /etc/init.d/apache2 restart

To test it, just use Firefox and point it at localhost or 127.0.0.1. The index page is at /var/www/index.html (the "document root" is /var/www). If you want to get a little fancier, you can grab PHP:

sudo apt-get install php5 libapache2-mod-php5 sudo /etc/init.d/apache2 restart

And put this in /var/www/test.php:

<?php phpinfo(); ?>

Go to localhost/test.php and it should print a bunch of details about your setup.

I spent a lot of time working on access control (Apache docs here). I should probably have read the VirtualBox guide first (here):

A virtual machine with NAT enabled acts much like a real computer that connects to the Internet through a router. The "router", in this case, is the VirtualBox networking engine, which maps traffic from and to the virtual machine transparently. The disadvantage of NAT mode is that, much like a private network behind a router, the virtual machine is invisible and unreachable from the outside internet; you cannot run a server this way unless you set up port forwarding (described below).

Actually, we will want to do this later. A lot of time was wasted because I followed the "simple" instructions (e.g. Apache's) and couldn't figure out why stuff didn't work.

It turns out that although /etc/apache2/apache2.conf and /etc/apache2/httpd.conf have relevant settings there are a bunch of other files (in directories under /etc/apache2/). If you want to change access to the document root or scripts, you have to modify /etc/apache2/sites-available/default. Guess I should have read the first page (here) of the Ubuntu Server docs. It's right there!

Here is part of the file in the original form:

ScriptAlias /cgi-bin/ /usr/lib/cgi-bin/ <Directory "/usr/lib/cgi-bin"> AllowOverride None Options +ExecCGI -MultiViews +SymLinksIfOwnerMatch Order allow,deny Allow from all </Directory>

I'm not sure what all of this does, but using ScriptAlias means that (Apache docs again):

The ScriptAlias directive tells Apache that a particular directory is set aside for CGI programs. Apache will assume that every file in this directory is a CGI program, and will attempt to execute it, when that particular resource is requested by a client.

I just modified /usr/lib/cgi-bin/ to be /home/te/cgi-bin/.

And, since we will ultimately care about access permissions, you should read (at least) this page of the docs, which explains that

Order allow,deny

is explained as:

Allow,Deny
First, all Allow directives are evaluated; at least one must match, or the request is rejected. Next, all Deny directives are evaluated. If any matches, the request is rejected. Last, any requests which do not match an Allow or a Deny directive are denied by default.
Deny,Allow
First, all Deny directives are evaluated; if any match, the request is denied unless it also matches an Allow directive. Any requests which do not match any Allow or Deny directives are permitted.

I copied from further down in the file (about access to /usr/share/docs) to the ScriptAlias directive above:

Order deny,allow Deny from all Allow from 127.0.0.0/255.0.0.0 ::1/128

The subnet masking stuff will need an explanation of its own.

So now, I can put the first script from the second part of this post (that prints data from os.environ in cgi-bin under my home directory, and it works.

Trying Ubuntu Linux (2)

I'm installing software on Ubuntu 11.0.4 (in a Virtual Box VM running under OS X Lion) and my standard test is matplotlib. To begin with, Ubuntu 11.0 has python but does not have numpy. So, let's start there.

I found build instructions here, but didn't actually follow them. Instead, I took their advice and used apt-get. While I was still thinking I would build everything I did:

sudo apt-get install gfortran sudo apt-get install python-dev

And then I went after BLAS, and LAPACK. Based on the library name on this page, I did:

sudo apt-get install libblas3gf

And using the library name from this page, I did:

sudo apt-get install liblapack3gf

Then, I just did:

sudo apt-get install python-numpy

And, amazingly enough, it works:

>>> import numpy >>>

Now:

sudo apt-get install zlib-bin sudo apt-get install libpng3 sudo apt-get install libfreetype6 sudo apt-get install python-matplotlib Failed to fetch http://us.archive.ubuntu.com/ubuntu/pool/main/t/tcl8.5/tcl8.5_8.5.9-2_amd64.deb Something wicked happened resolving 'us.archive.ubuntu.com:http' (-5 - No address associated with hostname) E: Unable to fetch some archives, maybe run apt-get update or try with --fix-missing?

Oops..

sudo apt-get update

And now it gets tcl8.5

What else to do for matplotlib?
Nothing!

We'll try the example from here, and we get a png in the right place, which I just double click and get the graphic at the top. Here is the whole Desktop:



I really can't take any credit for this. I am very impressed speechless with the Ubuntu dev guys.

[ UPDATE: Did the same apt-get for scipy and the test at the bottom of this post works too! ]

Trying Ubuntu Linux (1)

I thought I would play around with Linux a bit, trying to improve my Unix skills. I'm running it in a second partition on my Macbook that has OS X Lion installed (post here) and runs a Virtual Box VM.

I first heard about and tried Oracle's Virtual Box in connection with QIIME (post here).

I grabbed a VirtualBox version 4.1 binary installer for Mac OS X here.

Then I downloaded Ubuntu Linux 11.04 (now why is it called Natty Narwhal?) from here. There is another download page here but I think it has the same choices.

I ended up trying both Desktop "CD"s, because I was a little confused about the version I needed. There appears to be only 32-bit for i386 (ubuntu-11.04-desktop-i386.iso) but 64-bit for AMD (ubuntu-11.04-desktop-amd64.iso). However, the chip running the host OS is not relevant. The QIIME Linux is labeled AMD64 and is 64-bit and can run on my 64-bit machines under a VM. I suppose it would matter if I wanted to run Linux directly in the partition instead of in Virtual Box.

I have to confess, I've installed several different Linux versions (and then removed them) and the one described here a number of times, because I ran into problems with the install that I couldn't figure out how to back out of. But experimentation is relatively cheap, since it's just a VM. I don't have any worries about leftover crap on my disk. If worst comes to worst I can just erase the whole partition. Although it takes patience for the install steps.

---

Stock 64-bit Ubuntu Linux

So, I ran the Virtual Box installer, set the Name as Ubuntu, Operating System as Linux and the Version as Ubuntu 64. Then I set up the VM with all standard options (except 1 GB for RAM). Just click on through.

To actually install Linux using a .iso file, use the "First Run Wizard" by double-clicking the new machine in the list. For Select Installation Media: I navigated to the saved file. On one try I somehow got past the Wizard without properly installing and had to delete and re-configure the VM.

It takes quite a while to boot the first time, copying all the files over to the VM, I guess. Choose Erase disk and install Ubuntu (the disk size should be only the size of the VM disk). There are a bunch of dialogs to go through. After the install and reboot we get a message that Unity won't run on our hardware and so we're going to get Classic. Oh well.

Reboot..

Now, there are two things that seem really important to have but are not present in this stock install:

• a Shared Folder for file sharing between host and guest OS
• a Clipboard for copy/paste between host and guest OS

To have these features you must install what are called the "Guest Additions," which are additions to the Linux kernel. Some instructions are on this page:

After you have installed the Linux Guest system, there are a few additional packages needed for the Guest Additions. These packages enables you to create kernel modules. For Debian and Debian based distro's like Ubuntu, you need the following packages:
dkms
build-essential
linux-headers-generic

Find the Terminal (why not just drag it to the menu bar) and do:

sudo apt-get install dkms build-essential linux-headers-generic

You'll have to type in the command, since copy/paste doesn't work yet :)
When it's done:

Reboot to make the changes effective..

Next, "mount" the Guest Additions "iso": in the VB menu (in OS X) under Devices, choose Install Guest Additions. It shows up as a CD on the Desktop. Then from the dialog choose the default that automatically runs the shell script. Or hit Open AutoRun Prompt > OK and then in Terminal do:

media/VBOXADDITIONS{version}/VBoxLinuxAdditions.run

After it's done, reboot again..

At this point, copy/paste between host/guest should work. You'll have to click twice, once to switch focus to the VM and once to switch focus again to whatever application is receiving the paste.

There are two different key combos for Linux and neither is CMD. Not sure if this is standard Linux or it is just since the left CMD key is a special key for the VM. To copy and paste within Terminal on Linux use CTL-SHIFT-x and CTL-SHIFT-v rather than CMD (for Text Edit just use CTL-x and CTL-v).

Now we can try to get shared folders. In the Virtual Box window, from the Devices menu choose Shared Folders. Add a folder (navigate the host OS to the desired path). Also, make a folder under Linux that we'll use to share with, say /home/telliott/Desktop/share

In the dialog, choose Automount (and Make Permanent, see below). Even after the specification of a shared folder, you still have to mount it within Linux. The Qiime guys do:

sudo mount -a

Unfortunately, it doesn't seem to work for me. (Not in the Qiime install either). So then I tried this:

sudo mount -t vboxsf -o uid=1000,gid=1000 shared ~/Desktop/share

The first argument after vboxsf above is the name of the folder on the host OS (as defined by us in the Devices menu above). The second argument is the name of the sharepoint --- the path to a folder on the Guest.

It didn't work at first, because it turned out, I had failed to check Automount above. And, if you want to make a change to the Shared Folder setting, it seems that you have to reboot to be able to actually do it, otherwise there's a message about the resource being in use.

Also, I never did figure out how to change Make Permanent. Any change or attempt to delete in the Shared Folders dialog (even after a reboot) gives an error message: VERR_PERMISSION_DENIED. There's a note in the link about putting the mount command in

/etc/rc.local

but I don't see how that's going to work since it needs to run as root? Luckily, the VM allows me to suspend Linux in its current state. And now, thinking that I need to re-run that command every time upon booting, I try it and get:

sudo mount -t vboxsf -o uid=1000,gid=1000 shared ~/Desktop/share
/sbin/mount.vboxsf: mounting failed with the error: Invalid argument

After re-boot, the same mount fails! I am starting to get angry with this thing.

I gave up and set up Dropbox access from within the VM. It works great.

Now, to build some software.

PyMOL (1)

If you're interested in Python and bioinformatics you've undoubtedly heard of Pymol.

If you're at a place that cares about molecular biology (my institution does not), you may even have a site license. I have a copy that I got for educational use about 3 years ago, and still have some examples that I made for a class I taught. That version has a Cocoa GUI.

There is also an open source version that is up on SourceForge (here). I thought it would be a good challenge to try to build the open source Pymol.

[ Update: I think the Cocoa GUI one is MacPymol here). I just noticed this on the legacy page. ]

The dependencies are listed in the README as:

PyMOL has the following external dependencies: 1. OpenGL 2. GLUT library for OpenGL (or freeglut) 3. Python (v 2.1 or better) compiled with threads support 4. libpng (can compiled without it, but image saves won't work) 5. Tcl/Tk (for the external GUI) 6. Numerical Python (technically optional - now using numpy instead) 7. Python megawidgets (Pmw -- only required for the "external" GUI) 8. FreeType2 (if you want nice fonts)

Most of these are already present in Lion, so I didn't think it would be too difficult. What I failed to realize is that the Tcl/Tk in Lion is apparently messed up (brief SO discussion here).

I downloaded pymol-v1.4.1.tar.bz2 from the main page at SourceForge, and did the standard distutils thing:

python setup.py build

The first error was:

layer0/ShaderMgr.c:173: error: ‘GLEW_OK’ undeclared (first use in this function)

The top Google hit shows it's related to OpenGL and glut.
If present, there should be a file glew.h (and likely glut.h).

> find /usr -iname *glut* /usr/X11/include/GL/glut.h /usr/X11/include/GL/glutf90.h /usr/X11/lib/libglut.3.7.dylib /usr/X11/lib/libglut.3.dylib /usr/X11/lib/libglut.dylib /usr/X11/lib/pkgconfig/glut.pc > find /usr -iname *glew* > find /usr -iname *GLEW* >

No glew.

brew install glew

works and we have

$ find /usr/local/lib -iname *glew* /usr/local/lib/libGLEW.1.5.8.dylib /usr/local/lib/libGLEW.1.5.dylib /usr/local/lib/libGLEW.a /usr/local/lib/libGLEW.dylib /usr/local/lib/pkgconfig/glew.pc

and we find "GLEW_OK defined on line 15273 of glew.h. Unfortunately, the build fails with a link error at the end:

ld: warning: ignoring file /usr/local/lib/libGLEW.dylib, file was built for unsupported file format which is not the architecture being linked (i386) ld: duplicate symbol __CShaderMgr in build/temp.macosx-10.7-intel-2.7/layer2/RepCylBond.o and build/temp.macosx-10.7-intel-2.7/layer2/RepCartoon.o for architecture i386 > file /usr/local/lib/libGLEW.dylib /usr/local/lib/libGLEW.dylib: Mach-O 64-bit dynamically linked shared library x86_64

So it seems that Homebrew has built x86_64, which can't work with what pymol is building (i386). What I should do is rebuild libGLEW.dylib with the correct architecture. I haven't accomplished that yet. Instead, I noticed a page that discusses this very same problem. One answer is in a different thread (here):

Kasper, not sure if this will help you, but I think the GLEW stuff was
added recently in V1.4.x and it seems to be still somewhat experimental.
May be to get you started you could switch to v1.3r1 from the repository
and try to embedd that version. V1.3.x should not contain the new GL
stuff.

So let's try 1.3. I grabbed pymol-1.3r2-src.tar.bz2 from here and it got past the previous problem but then dies with:

/usr/X11R6/include/ft2build.h:56:38: error: freetype/config/ftheader.h: No such file or directory

It's looking in /usr/X11R6 and can't find freetype which is in /usr/X11 (my post here).

So.. I modify the setup.py in pymol:

inc_dirs=["ov/src", "layer0","layer1","layer2", "layer3","layer4","layer5", "/usr/X11R6/include", # my additions "/usr/X11/include", "/usr/X11/include/freetype2", "/usr/local/include/GL", EXT+"/include", EXT+"/include/GL", EXT+"/include/freetype2", "modules/cealign/src", "modules/cealign/src/tnt", ] libs=[] pyogl_libs = [] lib_dirs=[] def_macros=[("_PYMOL_MODULE",None), ("_PYMOL_LIBPNG",None), ("_PYMOL_FREETYPE",None), ] ext_comp_args=[] ext_link_args=[ "-L/usr/X11R6/lib", "-lGL", "-lXxf86vm", "-L"+EXT+"/lib", "-lpng", "-lglut", "-lfreetype", # my additions "-lGLEW" ]

And now it builds.

But.. after doing ./pymol the GUI hangs. I think it's Python (IDLE?), which I never use.

However, I can get pymol to run from the command line with no GUI:

> pymol -c PyMOL(TM) Molecular Graphics System, Version 1.3. .. Command mode. No graphics front end. Detected 2 CPU cores. Enabled multithreaded rendering. PyMOL: normal program termination.

I ran a short script and made a figure, it looks fine.

So then, we get wise words from Ned Deily (here), and see (here)

- You should not rely on the Apple-suppled Pythons if you want to use IDLE.

I thought about putting Tcl/Tk 8.4 in a different place, etc.

But finally I said (f*** it) and I reformatted my drive to add a new partition, and then used my USB flash drive with the Lion Installer (post here), and then I installed MacPorts and did

$ sudo port install pymol ---> Computing dependencies for pymol ---> Dependencies to be installed: freetype bzip2 zlib glew libpng mesa makedepend pkgconfig glib2 autoconf help2man gettext expat libiconv gperf ncurses ncursesw p5-locale-gettext perl5 perl5.12 perl5 perl5 m4 automake libtool xorg-xproto py27-libxml2 libxml2 libxml2 python27 db46 gdbm openssl python_select readline sqlite3 python27 xorg-dri2proto xorg-glproto xorg-libXfixes xorg-fixesproto xorg-libX11 xorg-bigreqsproto xorg-inputproto xorg-kbproto xorg-libXau xorg-libXdmcp xorg-libxcb xorg-libpthread-stubs xorg-xcb-proto xorg-util-macros xorg-xcmiscproto xorg-xextproto xorg-xf86bigfontproto xorg-xtrans xorg-libXi xorg-libXext xorg-libXmu xorg-libXt xorg-libsm xorg-libice py26-numeric python26 py26-pmw py26-tkinter tk Xft2 fontconfig xrender xorg-renderproto tcl xorg-libXScrnSaver xorg-scrnsaverproto xdpyinfo xorg-libXinerama xorg-xineramaproto xorg-libXtst xorg-recordproto xorg-libXxf86vm xorg-xf86vidmodeproto xorg-libdmx xorg-dmxproto

Incredible! Everything but the kitchen sink. Perhaps even more incredible, it worked. An hour or two later, we finish with:

---> Attempting to fetch pymol-1.4_2.darwin_11.x86_64.tbz2 from http://packages.macports.org/pymol ---> Fetching pymol ---> Verifying checksum(s) for pymol ---> Extracting pymol ---> Applying patches to pymol ---> Configuring pymol ---> Building pymol ---> Staging pymol into destroot ---> Installing pymol @1.4_2 ---> Activating pymol @1.4_2 ---> Cleaning pymol $

I got 1LMB.pdb from the PDB, and I ran this script with pymol script.pml:

load ~/Desktop/1LMB.PDB hide everything select D1, chain 1 select D2, chain 2 select R1, chain 3 select R2, chain 4 deselect all show sticks, D1 show sticks, D2 show cartoon, R1 show cartoon, R2 color red, R1 zoom active, -10 move x, 4 bg_color white ray 1200,1200 png ~/Desktop/x.png

and I got what's at the top of the post. Pretty cool! More later on this topic.