Principal Component Analysis with published CLL gene expression data (Herishanu et al, Blood 2011)

Over the last few years, my understanding of chronic lymphocytic leukaemia has developed. Leukaemic cell proliferation has been shown to be a key element of the disease. The site of CLL cell proliferation was an important question. A gene expression analysis published in Blood by Herishanu et al in 2011 showed a proliferative gene expression pattern in the leukeamic cells from the lymph nodes as distinct from other sites of the body (peripheral blood or bone marrow).

I've been exploring the gene expression data using R. As a first step, I've explored the principal component analysis similar to that shown in Figure 1A of the paper.

The PCA is convincing because the three groups of samples cluster among themselves.  The 3D plot is useful because the three components separate the groups while two dimensions do not.

There are other ways of visualising PCA plots...  Sources are indicated as part of the script. 

Here is the script that I used to analyse this:

SCRIPT START

# install.packages("RCurl", "scatterplot3d")

library(RCurl)

library(scatterplot3d)

# A few year's ago a gene expression profile paper was published in Blood by Herishanu et al

# that gave us a new insight into chronic lymphocytic leukamia. 

# http://www.bloodjournal.org/content/117/2/563.long

# I have written this script to explore the data 

# and to reproduce the principal component analysis. 

# I have downloaded the original CEL files from GEO

# http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21029

# they were placed in a folder and the expression values extracted 

# what's produced is Affymetrix IDs and log2 expression values 

# I have created two data sets

# we can reproduce the PCA workflow with a subset of the data is required (5000 genes selected randomly). 

# this will be faster to download and visualise. 

# get the data - needs connection to the internet

# to use remove the hashtag "#"

# x <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/herishanuMicroArrayDataSubset.tsv")

# here is a link for all the data (>50000 probes) 

x <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/herishanuMicroArrayData.tsv")

# turn the data into a useful data.frame

data <- read.table(text = x, header = TRUE, sep = "\t")

# read in a file with names of samples

# created this from cut and paste of web page combined with manipulation in Word & Excel

sampID <- read.csv("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/Herishanu_Samp_ID.csv", header = FALSE)

str(sampID)

sampID$V1 <- as.character(sampID$V1)

# is order of column names same as list in sampID?

data.col.name <- colnames(data)

data.col.name <- gsub(".CEL", "", data.col.name)

sampID.name <- sampID$V1

all.equal(data.col.name, sampID.name)

# TRUE so that is useful. 

# change column names from chip to PB, BM and LN

# for peripheral blood, bone marrow and lymph node respectively

sampID$Name <- paste0(as.character(sampID$V3), as.character(sampID$V5))

colnames(data) <- sampID$Name

# how do the sample cluster?

# I first need to find the "distances" between the arrays

# and show those distances in the hierarchical cluster plot

# I can pile up all the commands in one

# for more of a step by step approach, see these scripts:

# http://rforbiochemists.blogspot.co.uk/2016/02/visualising-some-cll-proteomic-data-for.html

# http://rforbiochemists.blogspot.co.uk/2016/03/gene-expression-analysis-and.html

plot(hclust(dist(t(data))))

# this is interesting because all the patient samples cluter together!

# with bone marrow and peripheral blood closer together than LN in all cases. 

# this is what it says in the paper and is shown in the supplementary data. 

# http://www.bloodjournal.org/content/117/2/563/tab-figures-only

# What's needed is a way to normalise by patient - HOW??

# Quote from the paper:

# "In the 12 patients in whom all 3 compartments had been arrayed, 

# the patient effect on gene expression was subtracted 

# by mean centering the expression value of each gene across the 3 compartments 

# for each patient separately."

# the twelve with all 3 compartments: #1, #2, #3, #4, #8, #9, #10, #11, #12, #13, #25, #26

# so if I understand this correctly, I need to extract a patient sample

# then mean centre all three samples, 

# and repeat for all 12 samples

# let's explore this concept with just two samples to see if it works

# pick out data for patient #26

samp26 <- data[,grep("#26", sampID$Name)]

# pick out data for patient #25

samp25 <- data[,grep("#25", sampID$Name)]

# put them together

twosamp <- cbind(samp26, samp25)

# visualise with a box plot - look normalised

boxplot(twosamp)

plot(hclust(dist(t(twosamp))))

# cluster by patient

# I think I understand what they have done.

# mean centred for each gene.... i.e. by row in the data

# For us this is really each probeset

# http://gastonsanchez.com/how-to/2014/01/15/Center-data-in-R/

# for each of the values substract the rowMeans

samp25.s <- samp25 - rowMeans(samp25)

samp26.s <- samp26 - rowMeans(samp26)

twosamp.s <- cbind(samp26.s, samp25.s)

plot(hclust(dist(t(twosamp.s))))

# This seems to work and gives a different cluster diagram

# with LN coming together nicely.

# apply this to the whole data set.

# extract each patient cohort

# these twelve were: #1, #2, #3, #4, #8, #9, #10, #11, #12, #13, #25, #26

# do #1 and # 2 first then do the other automatically...

# pat # 1

colnames(data)

pat_1 <- cbind(data[,1], data[,27], data[,46])

colnames(pat_1) <- c("PB#1", "BM#1", "LN#1")

pat_1.s <- pat_1 - rowMeans(pat_1)

# pat # 2

pat_2 <- cbind(data[,2], data[,28], data[,47])

colnames(pat_2) <- c("PB#2", "BM#2", "LN#2")

pat_2.s <- pat_2 - rowMeans(pat_2)

data.n <- cbind(pat_1.s, pat_2.s)

data.reqd <- c("#3", "#4", "#8", "#9", "#10", "#11", "#12", "#13", "#25", "#26")

for(i in 1:length(data.reqd)){

  samp <- data[,grep(data.reqd[i], sampID$Name)]

  samp <- samp - rowMeans(samp)

  data.n <- cbind(data.n, samp)

}

# now have a file called dat.exp.n - normalised within patients.

colnames(data.n)

plot(hclust(dist(t(data.n))))

# nice cluster by region with good separation of LN samples

# some overlap within PB and LN

# extract cell site to use for colours in PCA plots

names <- colnames(data.n)

colourby <- gsub("#\\d+", "", colnames(data.n))

colourby <- gsub("\\.", "", colourby)

colourby <- gsub("\\d+", "", colourby)

# this does cluster the LN samples distinctly but there is still some overlap 

# of Peripheral Blood and Bone Marrow. 

# PCA from the paper looks nice and convincing

# http://www.r-bloggers.com/computing-and-visualizing-pca-in-r/

# uses function princomp()

exp.pca <- princomp(data.n)  # this function does the PCA

print(exp.pca)

plot(exp.pca, type = "l")

Plot of variance by component for CLL gene expression patterns

plot(princomp(data.n)$loadings)

A simple 2D plot of component 1 vs component 2

# http://www.r-bloggers.com/visualizing-principal-components/

p <- princomp(data.n)

loadings <- p$loadings[]

p.variance.explained <- p$sdev^2 / sum(p$sdev^2)

# plot percentage of variance explained for each principal component    

barplot(100*p.variance.explained, las=2, xlab='', ylab='% Variance Explained')

Looking at the effects of each component as a bar plot.

#*****************************************************************

# 2-D Plot

#******************************************************************         

x <- loadings[,1]

y <- loadings[,2]

z <- loadings[,3]

cols <- as.factor(colourby)

cols <- gsub("PB", "green", cols)

cols <- gsub("BM", "red", cols)

cols <- gsub("LN", "blue", cols)

# pch = 24: triangle point-up

# pch = 22: square

# pch = 21: circle

symbols <- as.factor(colourby)

symbols <- gsub("PB", 24, symbols)

symbols <- gsub("BM", 21, symbols)

symbols <- gsub("LN", 22, symbols)

symbols <- as.numeric(symbols)

# plot loadings on the first and second principal components 

# identify sample by body location

plot(x, y, type='p', 

     pch=symbols, 

     xlab='Comp.1', ylab='Comp.2', 

     col = cols, main = "Principal Component Analysis - CLL cell gene expression")

A simple 2D plot of component 1 vs component 2
 Each sample is coloured by location

# add a legend to the top of the plot

legend("top",      # location

       bty="n",              # suppress legend box, shrink text 50%

       title="Body Location of sample",

       c("PB", "BM", "LN"), 

       pch=symbols, col = cols, horiz=TRUE)

# label up the points

text(x, y, colnames(data.n), col=cols, cex=.5, pos=4)

# Comment: lymph node samples separate nicely from other samples

# some overlap between bone marrow and peripheral blood 

#*****************************************************************

# 3-D Plot, for good examples of 3D plots

# http://statmethods.wordpress.com/2012/01/30/getting-fancy-with-3-d-scatterplots/

#******************************************************************                 

# plot all companies loadings on the first, second, and third principal components and highlight points according to the sector they belong

s3d = scatterplot3d(x, y, z, 

                    xlab='Comp.1', ylab='Comp.2', zlab='Comp.3', 

                    color=cols, pch = symbols,

                    main = "Principal Component Analysis in 3D \n CLL cell gene expression")

A 3D graph separates the different groups of samples. 

s3d.coords = s3d$xyz.convert(x, y, z)

text(s3d.coords$x, s3d.coords$y, 

     labels=colnames(data.n), col=cols, cex=.8, pos=4)

Labels can be added but I'm not sure they help in this example. 

# change the order and the angle to make it look a bit more like the figure

s3d = scatterplot3d(y, x, z, 

                    xlab='PC2', ylab='PC1', zlab='PC3',             

                    color=cols, pch = symbols,

                    angle=-25,

                    grid = FALSE,

                    main = "Principal Component Analysis in 3D \n CLL cell gene expression")

par(xpd=TRUE)    # allows legend outside the graph

# add legend

legend(10, -4.5,     # location

       bty="n",              # suppress legend box

       title="Site of sample",

       c("PB", "BM", "LN"), 

       pch=symbols, col = cols, horiz=TRUE)

# add source as text

text(4, -6.5, cex =0.7,

     "Source: Herishanu et al, Blood 2011 117:563-574; doi:10.1182/blood-2010-05-284984")

SCRIPT END

Gene Expression Analysis and Visualization for VizBi 2016 (Pt 1)

UPDATE: During the VizBi2016 tutorial session, the participants noticed a couple of errors in the script. I think these have now been corrected.
Also, Christian Hauer suggested that we visualise the distance matrix using the heatmap packages. He suggested pheatmap, in particular. Thanks Christian.

I have spent some time creating two quite long scripts that generate a selection of visualisation of some gene expression data generated recently by Mel Boyd at the Cardiff CLL Research Group. The script illustrates a workflow that results in a selection of interesting and useful visualisations that allow us to compare the different drug treatments used.

The data was generated from an Affymetrix HTA 2.0 chip. The data has been normalised and mapped using the oligo package. I have focussed on the protein coding transcripts and I have generated a random selection of 5000 of these genes.

This part of the experiment uses the data from four different experimental 'drugs' which target different pathways and have different cellular effects. This analysis and the associated visualisations are used to allow us to compare the drug treatments and look at the data as a whole set.

I am happiest with the visualisation that is generated by the loadings of the principal component analysis but there are some other good visualisations generated by the scripts too.

The cluster at the top right is from the drug that causes the most changes. The other drugs cause few changes and so map closer to the untreated samples.

Next script shows volcano plots and a heatmap.


Here is the script:
# START 
# download the data from github
library(RCurl)
x <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/microArrayData.tsv")
data <- read.table(text = x, header = TRUE, sep = "\t")

# basic commands for looking at the data
head(data)      # view first six rows
str(data)       # shows it's a dataframe
colnames(data)  # column names
View(data)      # this works in R-Studio

# for first part of analysis
# subset of the data with just the numbers.
new.data <- data[2:16]
# these are log2 values 
# The data has been normalized


## FIRST visualisation of the data: a simple box plot
boxplot(new.data, 
        las =2, # las = 2 turns the text around to show sample names
        ylab = expression(bold("expression")),
        main="Boxplot of expression data")  

## SECOND visualisation(s): calculate and visualise distances
# turn into a matrix 
data.m <- as.matrix(new.data)
# transpose the data because a distance matrix works in rows
data.m.t <- t(data.m)
# calculate the distances and put the calculations into an object called distances
distances <- dist(data.m.t)
# convert this distances object into a matrix. 
distances.m <- data.matrix(distances)
# you can look at this object.
View(distances.m)
# we can extract the size of the object and the titles
dim <- ncol(distances.m)
names <- row.names(distances.m)
# now to create the visualisation of the difference matrix. 
# first the coloured boxes
image(1:dim, 1:dim, distances.m, axes = FALSE, xlab = "", ylab = "")
# now label the axis
axis(3, 1:dim, names, cex.axis = 0.8, las=3)
axis(2, 1:dim, names, cex.axis = 0.8, las=1)
# add the values of the differences
text(expand.grid(1:dim, 1:dim), sprintf("%0.1f", distances.m), cex=1)

# heatmaps ccan also be used to visualise the distance matrix. 
# samples can be clustered or you can stop this. Thanks Christian Hauer
heatmap(distances.m)

heatmap3(distances.m)

pheatmap(distances.m)

## THIRD visualisation: use the distance matrix to do some clustering
plot(hclust(distances))

# three clusters - based on the distance matrix and shows the same thing really.

## FOURTH visualisation: principal component analysis
# Do a Principal Component Analysis  
# and widely-used technique for viewing patterns of gross variation
# in datasets
# We can see whether an array is substantially different to the others
pca <- princomp(data.m)  # function that does the PCA
summary(pca)  # one component accounts for 99% of the variance
plot(pca, type = "l")  

names <- factor(gsub("\\.\\d", "", names)) # change names into factor
plot(pca$loadings, col = names, pch = 19)
text(pca$loadings, cex = 0.7, label = colnames(new.data), pos =3)

# some of the treatment cluster well together (e.g. Drug D) others not so much


# can be interesting to look at a subset of the data...
pca <- princomp(data.m[,1:12])
plot(pca$loadings, col = names, pch = 19, 
     main = "Just three drug treatments")
text(pca$loadings, cex = 0.7, label = colnames(new.data)[1:12], pos =3)

# basically happy with our data generally. 
# next step in next script is to look at differentially expressed transcripts.

The chemical logic of glycolysis

Dr Sergio Martínez Cuesta from the EMBL-EBI has prepared the script below which he will demonstrate and discuss during our Training Day tomorrow. 

The script analyses the functional differences between glycolytic enzymes using principal component analysis (PCA) and hierarchical clustering. It generates some very nice heat maps. 

Here's one of the plots, produced with the pheatmap package

Analysing enzymes in terms of bond cleavage

SCRIPT
### Title: The chemical logic of glycolysis
### Author: Sergio Martínez Cuesta, EMBL-EBI, Hinxton, UK
### Aim: To visualise chemical similarities between enzymes reactions in glycolysis using latent variable models (PCA and CA) and clustering techniques

## Data input: load bond change data into R
data <- read.csv("http://science2therapy.com/data/bc.csv", header=TRUE, row.names=1, sep="\t", check.names=FALSE)

# bc.csv contains the frequency counts of chemical bond changes of the 10 glycolytic reactions as obtained from EC-BLAST
# header=TRUE and row.names=1 allow the first row and column to be considered the column and row names, sep="\t" indicates that bc.csv is tab-separated and check.names=FALSE disables the check of variable names


## Explore data
data
class(data)
dim(data)
colnames(data) # types of bond changes (formed and cleaved)
rownames(data) # enzyme ids, see data/glycolysis_enzymes.csv


## Barplots
# Which glycolytic enzyme catalyses more bond changes?
freqs_bc<-sort(rowSums(data), decreasing=TRUE)
freqs_bc
barplot(freqs_bc, col="white", ylab="Frequency of bond changes", xlab="Enzyme name")

# Which bond changes are more common in glycolysis?
freqs_enz<-sort(colSums(data), decreasing=TRUE)
freqs_enz
barplot(freqs_enz, col="white", ylab="Frequency of bond changes", xlab="Type of bond change")

# Average number of bond changes by glycolytic reaction
mean(rowSums(data))


## Principal component analysis (PCA)
fit <- princomp(data)
summary(fit)
biplot(fit, main="PCA")
loadings(fit)
fit$scores

# Additional packages: FactoMineR
install.packages("FactoMineR") # this only needs to be done once
library(FactoMineR)
fit <- PCA(data, scale.unit = FALSE)


## Correspondence analysis (CA)
install.packages("ca")  # this only needs to be done once
library(ca)
fit <- ca(data)
summary(fit)
plot(fit)

# using library(FactoMineR)
# fit<-CA(data)


## Calculate bond change distances between glycolytic reactions
?dist
d <- dist(data, method = "manhattan")


## Hierarchical clustering
fit <- hclust(d, method="ward")
plot(fit, main="Hierarchical clustering")
groups <- cutree(fit, k=4)
rect.hclust(fit, k=4, border="red") 

## Heatmaps
heatmap(as.matrix(data))

heatmap(as.matrix(d))


## Pretty heatmaps
install.packages("pheatmap")  # this only needs to be done once
library(pheatmap)
pheatmap(data, 
         fontsize=13, 
         cellwidth = 15, 
         cellheight = 15, 
         color = colorRampPalette(c("white", "red"))(4),
         clustering_distance_rows=d,
         legend_breaks=c(0,1,2,3))


## How to save pdf images:
pdf("mygraph.pdf")
pheatmap(data)
dev.off()


# References:
# 1- Quick-R: http://www.statmethods.net/

# 2- Crawley, M. J. (2007). The R Book. Chichester, UK: John Wiley & Sons, Ltd. doi:10.1002/9780470515075

# 3- Rahman, S. A., Martínez Cuesta, S., Furnham, N., Holliday, G. L., & Thornton, J. M. (2014). EC-BLAST: a tool to automatically search and compare enzyme reactions. Nature methods, 11(2), 171–174. doi:10.1038/nmeth.2803. url: http://www.ebi.ac.uk/thornton-srv/software/rbl/

# 4- Husson, F., Lê, S., and Pagès, J. (2011) Exploratory multivariate analysis by example using R, [online] http://books.google.com/books?hl=en&lr=&id=EDnNBQAAQBAJ&oi=fnd&pg=PP1&dq=Exploratory+Multivariate+Analysis+by+Example+Using+R&ots=I4bnnmgPJw&sig=ENW39pnO00fOOSO_F01cQ_Iv-UA (Accessed February 6, 2015).

# 5- Keller, M. A., Turchyn, A. V., & Ralser, M. (2014). Non-enzymatic glycolysis and pentose phosphate pathway-like reactions in a plausible Archean ocean. Molecular systems biology, 10, 725. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/24771084