Bioonductor: summary

[ Update: Sorry for the dup. I tried using the new formatting on this one, but found a problem, which I hope is fixed now.]

I promise this is the last one. Here we load the sample data from ExpressionSet. I show (once again) how to obtain all the relevant stuff: expression data, sample names, feature names, and gene names.

One of the nice things about ExpressionSet objects is the very fancy indexing one can do with them. A great example is given on p. 166 of the Bioconductor book.

I made an stab at implementing range selectors here but it's almost trivial. Also, in this one, the Python code does not contain any informative whitespace, so I've used the <code> tag.
The script is first and then the output. If you can't figure something out, just ask.

script.py

import rpy2.robjects as robjects from rpy2.robjects.packages import importr import rpy2.rinterface as ri import bioc.biobase as biobase r = robjects.r g = robjects.globalenv #gd = importr('grDevices') dim = r['dim'] r('data("sample.ExpressionSet")') sample_ExpressionSet = r['sample.ExpressionSet'] eset = sample_ExpressionSet #print tuple(ri.globalenv) print eset print dim(eset) ed = eset.exprs print ed.rx(1,1) # no way to do es[x,y] from Python? # just wrap it r(''' do_sel <-function(obj,sel,bycol=FALSE) { if (bycol) { return(obj[,sel]) } else { return(obj[sel,]) } } ''') do_sel = r['do_sel'] # both work: v = robjects.IntVector([2,3,6]) v = robjects.BoolVector([True,False]*13) eset_sub = do_sel(eset,v,bycol=True) eset_sub = do_sel(eset_sub,v) print dim(eset_sub) # no explicit featureData, just from data rows rownames = r['rownames'] fd = rownames(ed) print fd[:2] # get the gene names importr('annaffy') r(''' get.genenames <- function(probeids) { symbols <- aafSymbol(probeids, 'hgu95av2') getText(symbols) } ''') get_genenames = r['get.genenames'] gn = get_genenames(fd) # slightly fancy, no dups, original order from collections import OrderedDict z = zip(gn,[None]*len(gn)) D = OrderedDict(z) L = D.keys() for i in range(4): print L[i*5:i*5+5] print pd = eset.do_slot('phenoData') data = pd.do_slot('data') v = robjects.IntVector(range(1,4)) print do_sel(data,v)

> python script.py 
ExpressionSet (storageMode: lockedEnvironment) 
assayData: 500 features, 26 samples  
element names: exprs, se.exprs  
protocolData: none 
phenoData 
sampleNames: A B ... Z (26 total) 
varLabels: sex type score 
varMetadata: labelDescription 
featureData: none 
experimentData: use 'experimentData(object)' 
Annotation: hgu95av2  

Features  Samples  
500       26  

[1] 192.742 

Features  Samples  
3        3  

[1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" 


******* 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).   


Loading required package: hgu95av2 
['', 'STAT1', 'GAPDH', 'ACTB', 'TFRC'] 
['C15orf31', 'ACTL6B', 'GLRA1', 'KCNB2', 'MGAT5'] 
['BMP3', 'RPL14', 'ZCWPW2', 'KITLG', 'GPR12'] 
['TRA@', 'ATP8A2', 'C1orf68', 'FAM20B', 'SLC34A1'] 

sex    type score 
A Female Control  0.75 
B   Male    Case  0.40 
C   Male Control  0.73 

>  

Bioconductor: multiple plots

Actually, I have a bit more to do with Bioconductor ExpressionSets, as long as I'm thinking about it. When I first played with the ALL example, I tried other pairs besides ('ALL1/AF4','BCR/ABL')---I'd have to look now, but they are probably in the original paper too. Anyway, in the first talk I gave about this stuff, I made plots of all the pairs. Some are really nice. So that's what I'd like to show now.

The difference is that then, I wrote the code in R, and I did not like the experience.

In Python, we do the standard imports. We factor out the R code into several functions, then load it from a text file in Python, and execute the code in R-space. In Python, we pick pairs to compare and call a plotting function in R.

The R code is a little hackish (I'm thinking of color.map), but this is a good start.

R code:

color.map <- function(g) { colors = c('red','green','blue', 'cyan','magenta','maroon') if (g == 'ALL1/AF4') i=1 if (g == 'BCR/ABL') i=2 if (g == 'E2A/PBX1') i=3 if (g == 'NEG') i=4 if (g == 'NUP-98') i=5 if (g == 'p15/p16') i=6 colors[i] } library('limma') library('annaffy') get.genenames <- function(eset) { probeids <- featureNames(eset) symbols <- aafSymbol(probeids, 'hgu95av2') getText(symbols) } compare.factors <- function(eset,L) { eset.sub = eset[,eset$mol.biol %in% L] f <- factor(as.character(eset.sub$mol.biol)) design <- model.matrix(~f) fit <- eBayes(lmFit(eset.sub,design)) pv = fit$p.value[,2] sel <- p.adjust(pv) < 0.001 eset.sel <- eset.sub[sel,] } plot.eset <-function(eset,...) { patient.colors <- unlist(lapply( eset$mol.biol,color.map)) genenames = get.genenames(eset) heatmap(exprs(eset), col=topo.colors(100), ColSideColors=patient.colors, labRow=genenames) }

The R code is in a file Rcode.txt. Python code:

import rpy2.robjects as robjects from rpy2.robjects.packages import importr import bioc.biobase as biobase r = robjects.r g = robjects.globalenv gd = importr('grDevices') importr('ALL') r('data("ALL")') ALL = g['ALL'] eset = biobase.ExpressionSet(ALL) #-------------------------------------------- FH = open('Rcode.txt','r') s = FH.read() FH.close() r(s) compare_factors = r['compare.factors'] plot_eset = r['plot.eset'] #-------------------------------------------- targets = [ ['ALL1/AF4','NEG'], ['ALL1/AF4','BCR/ABL'], ['E2A/PBX1','BCR/ABL'] ] def get_ofn(t): #print t, len(t) u,v = t s = u + '_' + v ofn = s.replace('/','-') + '.pdf' return ofn def plot(eset_sub,t): ofn = get_ofn(t) gd.pdf(ofn) plot_eset(eset_sub,ofn) gd.dev_off() for t in targets: eset_sub = compare_factors(eset,t) plot(eset_sub,t)

The first plot is at the top, and the second and third are below. It's pretty clear that there are significant commonalities according to mol.biol. BTW, I adjusted the target p-value to keep the number of hits down.

I also added code to substitute the gene names in the plots. That's pretty cool.

There are problems with the alternatives. That's for another time.

Learning to use RPy2 (5)

Here is the result of my experiment with the rpy2-bioconductor-extensions-0.2.


import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
import bioc.biobase as biobase
r = robjects.r
g = robjects.globalenv

importr('ALL')
r('data("ALL")')
data = g['ALL']
eset = biobase.ExpressionSet(data)
print eset

Run as a script.

> python script.py

First (above) section output:

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

Part 2:

expr_data = eset.exprs
print len(expr_data)
L = tuple(expr_data[:5])
print [round(t,3) for t in L]

output:

1616000
[7.597, 5.046, 3.9, 5.904, 5.925]

Part 3:

print list(eset.slotnames())
phenoData = eset.do_slot('phenoData')
print 'phenoData', id(phenoData)
print phenoData

output:

['experimentData', 'assayData', 'phenoData', 'featureData', 'annotation', 'protocolData', '.__classVersion__']
phenoData 0x10bbc35a8
An object of class "AnnotatedDataFrame"
sampleNames: 01005 01010 ... LAL4 (128 total)
varLabels: cod diagnosis ... date last seen (21 total)
varMetadata: labelDescription

Part 4:

pd = eset.get_phenodata()
print 'pd', id(pd)
print pd
print list(pd.slotnames())
print

output:

pd 0x10bbc3488
An object of class "AnnotatedDataFrame"
sampleNames: 01005 01010 ... LAL4 (128 total)
varLabels: cod diagnosis ... date last seen (21 total)
varMetadata: labelDescription

['varMetadata', 'data', 'dimLabels', '.__classVersion__']

Part 5:

data = pd.do_slot('data')
print 'data'
print type(data)
mb = pd.pdata.rx2('mol.biol')
print list(mb)[:10]
levels = r['levels']
print levels(mb)
print list(data.rx2('mol.biol'))[:10]

output:

data
<class 'rpy2.robjects.vectors.DataFrame'>
[2, 4, 2, 1, 4, 4, 4, 4, 4, 2]
[1] "ALL1/AF4" "BCR/ABL" "E2A/PBX1" "NEG" "NUP-98" "p15/p16"

[2, 4, 2, 1, 4, 4, 4, 4, 4, 2]

Part 6:

fd = eset.get_featuredata
print fd
# fd() gives AttributeError, no _featureData
#
featureNames = robjects.r['featureNames']
fn = featureNames(eset)
print len(fn), fn[0]

output:

<bound method ExpressionSet.get_featuredata of >
12625 1000_at

In the last part, the "new" way doesn't seem to work yet. So I did as before.
Well, that's it for Bioconductor and Rpy2 for now. The presentation could be a lot better organized, but at least we figured out how to do most things.

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.

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.