Overview

The made4 function overview() provides a quick way to get an overview or feel for data. overview() will draw a boxplot, histogram and dendrogram of a hierarchical analysis. Hierarchical clustering is produced using average linkage clustering with a Pearson correlation measure of similarity (Eisen et al. 1998) This gives a quick first glance at the data.

overview(k.data)

Often its useful to label the samples using a class vector or covariate of interest, in this case, the tumour type (EWS, BL, NB or RMS).

overview(k.data, labels=k.class)

Overview of Khan data.

Figure 1: A) dendrogram showing results of average linkage clustering, B) boxplot and C) histrogram.

Often one will known classes in the data (eg Normal v Treatment, or different tumor types). We can insert a class colourbar under the dendrogram, and colour the boxplot.

overview(k.data, classvec=k.class, labels=k.class)

Overview of Khan data.

1. dendrogram showing results of average linkage clustering, B) boxplot and C) histrogram. In this case we have added a vector of class (classvec) to color the overview by class membership.

Visualising Results

There are many functions in ade4 and made4 for visualising results from ordination analysis. The simplest way to view the results produced by ord is to use plot. plot(k.ord) will draw a plot of the eigenvalues, along with plots of the variables (genes) and a plot of the cases (microarray samples). In this example Microarray samples are colour-coded using the classvec khan$train.classes which is saved as k.class.

k.class

##  [1] EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS   
## [11] EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS    EWS   
## [21] EWS    EWS    EWS    BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL BL-NHL
## [31] BL-NHL NB     NB     NB     NB     NB     NB     NB     NB     NB    
## [41] NB     NB     NB     RMS    RMS    RMS    RMS    RMS    RMS    RMS   
## [51] RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS    RMS   
## [61] RMS    RMS    RMS    RMS   
## Levels: EWS BL-NHL NB RMS

plot(k.coa, classvec=k.class, genecol="grey3")

Figure 2 Correspondence analysis of Khan dataset. A. plot of the eigenvalues, B. projection of microarray samples from patient with tumour types EWS (red), BL(blue), NB (green) or RMS (brown), C. projection of genes (gray filled circles) and D. biplot showing both genes and samples. Samples and genes with a strong associated are projected in the same direction from the origin. The greater the distance from the origin the stronger the association.

Genes and array projections can also be plotted using plotgenes and plotarrays. The function s.groups required a character vector or factor that indicates the groupings or classes (classvec) and allowed groups to be coloured in different colours. For example, to plot microarray samples (cases),

plotgenes(k.coa)  

To plot microarray samples, colour by group (tumour type) as specified by khan$train.classes

plotarrays(k.coa, arraylabels=k.class)

## [1] "Need to specify groups"

Alternative you can run these analysis and give a class vector to ord and it will automatically colour samples by this class vector

k.coa2<-ord(k.data, classvec=k.class)
plot(k.coa2)

Plot gene projections without labels (clab=0). Typically there are a large number of genes, thus it is not feasible to label all of these. The function plotgenes is more useful to use if you wish to add labels when there are lots of variables (genes)

The gene projections can be also visualised with plotgenes. The number of genes that are labelled at the end of the axis can be defined. The default is 10.

plotgenes(k.coa, n=5, col="red")

By default the variables (genes) are labelled with the rownames of the matrix. Typically these are spot IDs or Affymetrix accession numbers which are not very easy to interpret. But these can be easily labeled by your own labels. For example its often useful to labels using HUGO gene symbols. We find the Bioconductor annotate or biomaRt annotation packages are useful for this.

In this example we provide annotation from the Source database in khan$annotation. The gene symbol are in the vector khan$annotation$Symbol. (In this case its a factor)

gene.symbs<- khan$annotation$Symbol 
gene.symbs[1:4]

## [1] FDFT1 LYN   CKS2  MT1X 
## 1944 Levels:  101F6 182-FIP 7h3 A2LP AADAC AAMP ABCC3 ABCE1 ABCF1 ABCF2 ... ZWINT

plotgenes(k.coa, n=10, col="red", genelabels=gene.symbs)

Figure 3 Projection of genes (filled circles) in Correspondence analysis of Khan dataset. The genes at the ends of each of the axes are labelled with HUGO gene symbols.

To get a list of variables at the end of an axes, use topgenes. For example, to get a list of the 5 genes at the negative and postive end of axes 1.

topgenes(k.coa, axis = 1, n=5)   

To only the a list of the genes (default 10 genes) at the negative end of the first axes

topgenes(k.coa, labels=gene.symbs, end="neg") 

##  [1] "PTPN13"  "OLFM1"   "TNFAIP6" "GYG2"    "CAV1"    "MYC"     "FVT1"   
##  [8] "FCGRT"   "TUBB5"   "MYC"

Two character vectors can be compared using (the very poorly named function) comparelists

To visualise the arrays (or genes) in 3D either use do3d or html3d. do3d is a wrapper for scatterplot3d, but is modified so that groups can be coloured.

do3d(k.coa$ord$co, classvec=k.class, cex.symbols=3)

html3d produces a “pdb” output which can be visualised using rasmol or chime, a free interface for colour, rotating, zooming 3D graphs. The output from Output from html3D can be rotated and visualised on web browsers that can support chime (or a pdb viewer)

html3D(k.coa$ord$co, k.class, writehtml=TRUE)

Output from html3D

explor

It is also worth exploring the package explor which provides a nice R Shiny interface for browsing results from made4 and ade4. It is available from CRAN or from juba.github.io

Classification and Class Prediction using Between Group Analysis

Between Group Analysis (BGA) is a supervised PCA method (Dolédec and Chessel 1987) for classification/prediction. The basis of BGA is to ordinate the groups rather than the individual samples.

In tests on two microarray gene expression datasets, BGA performed comparably to supervised classification methods, including support vector machines and artifical neural networks (Culhane et al. 2002). To train a dataset, use bga, the projection of test data can be assessed using suppl. One leave out cross validation can be performed using bga.jackknife. See the BGA vignette for more details on this method.

k.bga<-bga(k.data, type="coa", classvec=k.class)

plot(k.bga, genelabels=gene.symbs) # Use the gene symbols earlier

Figure 5 Between group analysis of Khan dataset. A. Between.graph of the microarray samples, showing their separation on the discriminating BGA axes, B. Scatterplot of the first 2 axes of microarray samples, coloured by their class, C. graph of positions of genes on the same axis. Genes at the ends of the axis are most discriminating for that group

Sometimes its useful to visualise 1 axes of an analysis. To do this use graph1D or between.graph. The latter function is specifically for visualising results from a bga as it shows the separation of classes achieved.

between.graph(k.bga, ax=1)  # Show the separation on the first axes(ax)

Meta-analysis of microarray gene expression

Coinertia analysis cia (Dolédec S and Chessel D 1994) has been successfully applied to the cross-platform comparison of microarray gene expression datasets Meng et al. (2016), . CIA is a multivariate method that identifies trends or co-structure in the variance of multiple datasets which contain the same samples. That is either the rows or the columns of a matrix must be “matchable”. CIA can be applied to datasets where the number of variables (genes) far exceeds the number of samples (arrays) such is the case with microarray analyses. cia calls coinertia in the ade4 package. See the CIA vignette for more details on this method.

There are extended version of CIA in the Bioconducor packages mogsa and omicade4 Meng et al. (2016)

# Example data are "G1_Ross_1375.txt" and "G5_Affy_1517.txt"
data(NCI60)
coin <- cia(NCI60$Ross, NCI60$Affy)
names(coin)

## [1] "call"      "coinertia" "coa1"      "coa2"

coin$coinertia$RV

## [1] 0.7859656

The RV coefficient $RV which is 0.786 in this instance, is a measure of global similarity between the datasets.
The greater (scale 0-1) the better.

plot(coin, classvec=NCI60$classes[,2], clab=0, cpoint=3)

Figure 6 Coinertia analysis of NCI 60 cell line Spotted and Affymetrix gene expression dataset. The same 60 cell lines were analysed by two different labs on a spotted cDNA array (Ross) and an affymetrix array (Affy). The Ross dataset contains 1375 genes, and the affy dataset contains 1517. There is little overlap betwen the genes represented on these platforms. CIA allows visualisation of genes with similar expression patterns across platforms. A) shows a plot of the 60 microarray samples projected onto the one space. The 60 circles represent dataset 1 (Ross) and the 60 arrows represent dataset 2 (affy). Each circle and arrow are joined by a line, the length of which is proportional to the divergence between that samples in the two datasets. The samples are coloured by cell type. B) The gene projections from datasets 1 (Ross), C) the gene projections from dataset 2 (Affy). Genes and samples projected in the same direction from the origin show genes that are expressed in those samples. See vingette for more help on interpreting these plots.

Data Input

Example datasets provides with made4

Classification and class prediction using Between Group Analysis

Meta analysis of two or more datasets using Coinertia Analysis

Graphical Visualisation of results: 1D Visualisation

Graphical Visualisation of results: 2D Visualisation

Graphical Visualisation of results: 3D Visualisation

Interpretation of results