We will show here a typical pipeline for analysis of gene expression data. This tutorial is based on the use of Monocytes and Macrophage data from the following paper:

van de Laar L, Saelens W, De Prijck S, Martens L et al.Ā Yolk Sac Macrophages, Fetal Liver, and Adult Monocytes Can Colonize an Empty Niche and Develop into Functional Tissue-Resident Macrophages. Immunity 2016 Apr 19;44(4):755-68. PMID: 26992565

This tutorial contains the following steps:

  1. Download the data from GEO
  2. Process the expression table for analyses
  3. DE Analysis by Limma
  4. PCA (Principal component analysis)
  5. Clustering by Heatmap (only top 250 genes)
  6. GO (Gene Ontology) Enrichment Analysis
  7. GSEA (Gene Set Enrichment Analysis)
  8. Combining other dataset for PCA and clustering

You can find 3 files given to you under handouts directory:

File Description
handout.pdf Step-by-step explanation of the tasks.
Imun_gen.ex.Rdump Preprocessed object from GEO for loading in the analysis.
pheno.table A table storing the phenotype data for comparison.
mouse_hallmark_genesets.rdata Hallmark mouse gene sets for GSEA

Set the working directory

Before you start to work on this handout, please set your working directory to the place where you save the directory. For example,

setwd("/home/grasshoff/Documents/BIR_course")
# Please use your path instead of the above example

Loading the libraries

rm(list=ls())
if("mogene10sttranscriptcluster.db" %in% rownames(installed.packages()) == FALSE) {
  BiocManager::install(c("mogene10sttranscriptcluster.db"))
}
if("msigdbr" %in% rownames(installed.packages()) == FALSE) {
  BiocManager::install("msigdbr")
}

library(Biobase)
library(GEOquery)
library(limma)
library(mogene10sttranscriptcluster.db)
library(RColorBrewer)
library(gProfileR)
library(genefilter)
library(openxlsx)
library(piano)
library(gplots)
library(sva)
library(msigdbr)

Biobase contains base functions for Bioconductor; GEOquery is the bridge between GEO and BioConductor; limma is a library for the analysis of gene expression microarray data, especially the use of linear models for analysing designed experiments and the assessment of differential expression; mogene10sttranscriptcluster.db contains the mouse gene annotation of mapping between different systems; RColorBrewer can make colorful palettes; gProfileR is the package to perform gene enrichment analysis; openxlsx is for creating excel table in R; piano is the package for GSEA; gplots contains many plotting tools; and sva (ComBat) is used here to remove known batch effects from microarray data.

1 Download the data from GEO

After you find a dataset, which you want to explore, you can search in GEO using their GEO number. All the details of the protocal, platform and processing steps can be found in the GEO webpage.

Please take a look of this example: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76999

First, you need to download the data from GEO by the GEO number in R. We use the function getGEO (from package GEOquery) to download a list of all the available data.

# Load the data from GEO for GSE76999
GEO_List <- getGEO(filename="GSE76999_series_matrix.txt.gz", GSEMatrix = FALSE, getGPL=FALSE)
# Check the type of this variable 
class(GEO_List)
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
# Check the length of this list
length(GEO_List)
## [1] 1
# An element of a list always comes together with a name, right?
names(GEO_List)
## NULL
# Since there is only one ExpressionSet, let's assign it into another variable
eset <- GEO_List

ExpressionSet (eset here) is a Container for high-throughput assays and experimental metadata.

In order to understand ExpressionSet, type:

eset # Get an overview
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 35556 features, 36 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: GSM2042244 GSM2042245 ... GSM2042279 (36 total)
##   varLabels: title geo_accession ... tissue:ch1 (40 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
##   pubMedIds: 26992565 
## Annotation: GPL6246
varLabels(eset) # To see what kinds of labels inside
##  [1] "title"                   "geo_accession"          
##  [3] "status"                  "submission_date"        
##  [5] "last_update_date"        "type"                   
##  [7] "channel_count"           "source_name_ch1"        
##  [9] "organism_ch1"            "characteristics_ch1"    
## [11] "characteristics_ch1.1"   "characteristics_ch1.2"  
## [13] "treatment_protocol_ch1"  "growth_protocol_ch1"    
## [15] "molecule_ch1"            "extract_protocol_ch1"   
## [17] "label_ch1"               "label_protocol_ch1"     
## [19] "taxid_ch1"               "hyb_protocol"           
## [21] "scan_protocol"           "description"            
## [23] "data_processing"         "data_processing.1"      
## [25] "data_processing.2"       "platform_id"            
## [27] "contact_name"            "contact_email"          
## [29] "contact_department"      "contact_institute"      
## [31] "contact_address"         "contact_city"           
## [33] "contact_zip/postal_code" "contact_country"        
## [35] "supplementary_file"      "data_row_count"         
## [37] "relation"                "age:ch1"                
## [39] "strain:ch1"              "tissue:ch1"
eset$title # To see the names of the experiments
##  [1] "Monocyte extracted from adult (wk6-12) Bone Marrow, biological replicate 1"                                     
##  [2] "Monocyte extracted from adult (wk6-12) Bone Marrow, biological replicate 2"                                     
##  [3] "Monocyte extracted from adult (wk6-12) Bone Marrow, biological replicate 3"                                     
##  [4] "Monocyte extracted from adult (wk6-12) Bone Marrow, biological replicate 4"                                     
##  [5] "Monocyte extracted from E15.5 Fetal Liver, biological replicate 1"                                              
##  [6] "Monocyte extracted from E15.5 Fetal Liver, biological replicate 2"                                              
##  [7] "Monocyte extracted from E15.5 Fetal Liver, biological replicate 3"                                              
##  [8] "Monocyte extracted from E15.5 Fetal Liver, biological replicate 4"                                              
##  [9] "Macrophage extracted from E12.5 Yolk Sac, biological replicate 1"                                               
## [10] "Macrophage extracted from E12.5 Yolk Sac, biological replicate 2"                                               
## [11] "Macrophage extracted from E12.5 Yolk Sac, biological replicate 3"                                               
## [12] "Macrophage extracted from E12.5 Yolk Sac, biological replicate 4"                                               
## [13] "Alveolar Macrophage extracted from adult (wk6-12) mice via Bronchoalveolar lavage, biological replicate 1"      
## [14] "Alveolar Macrophage extracted from adult (wk6-12) mice via Bronchoalveolar lavage, biological replicate 2"      
## [15] "Alveolar Macrophage extracted from adult (wk6-12) mice via Bronchoalveolar lavage, biological replicate 3"      
## [16] "Alveolar Macrophage extracted from adult (wk6-12) mice via Bronchoalveolar lavage, biological replicate 4"      
## [17] "Alveolar Macrophage derived from transferred Bone Marrow Monocytes, 6 wks post-transfer, biological replicate 1"
## [18] "Alveolar Macrophage derived from transferred Bone Marrow Monocytes, 6 wks post-transfer, biological replicate 2"
## [19] "Alveolar Macrophage derived from transferred Bone Marrow Monocytes, 6 wks post-transfer, biological replicate 3"
## [20] "Alveolar Macrophage derived from transferred Bone Marrow Monocytes, 6 wks post-transfer, biological replicate 4"
## [21] "Alveolar Macrophage derived from transferred Fetal Liver Monocytes, 6 wks post-transfer, biological replicate 1"
## [22] "Alveolar Macrophage derived from transferred Fetal Liver Monocytes, 6 wks post-transfer, biological replicate 2"
## [23] "Alveolar Macrophage derived from transferred Fetal Liver Monocytes, 6 wks post-transfer, biological replicate 3"
## [24] "Alveolar Macrophage derived from transferred Fetal Liver Monocytes, 6 wks post-transfer, biological replicate 4"
## [25] "Alveolar Macrophage derived from transferred Yolk Sac Macrophages, 6 wks post-transfer, biological replicate 1" 
## [26] "Alveolar Macrophage derived from transferred Yolk Sac Macrophages, 6 wks post-transfer, biological replicate 2" 
## [27] "Alveolar Macrophage derived from transferred Yolk Sac Macrophages, 6 wks post-transfer, biological replicate 3" 
## [28] "Alveolar Macrophage derived from transferred Yolk Sac Macrophages, 6 wks post-transfer, biological replicate 4" 
## [29] "Alveolar Macrophage derived from transferred Alveolar macrophages, 6 wks post-transfer, biological replicate 1" 
## [30] "Alveolar Macrophage derived from transferred Alveolar macrophages, 6 wks post-transfer, biological replicate 2" 
## [31] "Alveolar Macrophage derived from transferred Alveolar macrophages, 6 wks post-transfer, biological replicate 3" 
## [32] "Alveolar Macrophage derived from transferred Alveolar macrophages, 6 wks post-transfer, biological replicate 4" 
## [33] "Alveolar Macrophage extracted from 6 wk old wild type mice via Bronchoalveolar lavage, biological replicate 1"  
## [34] "Alveolar Macrophage extracted from 6 wk old wild type mice via Bronchoalveolar lavage, biological replicate 2"  
## [35] "Alveolar Macrophage extracted from 6 wk old wild type mice via Bronchoalveolar lavage, biological replicate 3"  
## [36] "Alveolar Macrophage extracted from 6 wk old wild type mice via Bronchoalveolar lavage, biological replicate 4"
eset$organism_ch1 # Show the organism of the experiments
##  [1] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
##  [6] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [11] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [16] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [21] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [26] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [31] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [36] "Mus musculus"

2 Process the expression table for analysis

Here we select only the first 12 samples for downstream analysis and eliminate others. These 12 smaples are BM-MOs (bone marrow monocytes), FL-MOs (fetal liver monocytes), and YS-Macs (yolk sac-derived macrophages), with 4 replicates for each.

# Select the samples we need
sel <- 1:12
eset <- eset[ ,sel]
# We also need to modify its expression matrix
exprs(eset) <- exprs(eset)[ ,sel]
# Create labels for later steps
labels <- c("BM_MOs","FL_MOs","YS_Macs")
sel_groups <- factor(rep(labels, each=4))

# Lots of colors you can choose: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
# But let's try the most simple ones
# palette function is used to set the customized colors
palette(c("red","blue", "green"))

Before we do any further analysis, we need to do quality control of the data. We use boxplot to check whether we need further normalization on the data.

boxplot(exprs(eset), col=sel_groups)

How is it? Are you satisfied? I am not. In order to learn all the parameters you can adjust for this boxplot, you can read the help message by typing:

?boxplot

I will shown a few parameters that can be used to improve this plot. Letā€™s rotate the labels on x axis!!

boxplot(exprs(eset), col=sel_groups, las=2)

A little bit better, but not complete. Letā€™s make the label smaller.

boxplot(exprs(eset), col=sel_groups, las=2, cex.axis = 0.7)

I am a little bit picky and donā€™t like the fat box. I want to change its ratio and remove the outliers.

boxplot(exprs(eset), col=sel_groups, las=2, cex.axis = 0.7, boxwex=0.6, outline=FALSE)

The boxplot indicates all samples have the similar distribution and are normalized. This indicates there is no need to perfrom normalization, as normalization was already performed inthe data deposited in GEO. In case you need to normalize the data, this can be done with the following commands:

# normalizeBetweenArrays is a function in LIMMA
exprs_matrix.norm <- normalizeBetweenArrays(exprs(eset))
boxplot(exprs_matrix.norm, col=sel_groups, las=2, cex.axis = 0.7, boxwex=0.6, outline=FALSE)

3 Differential Expression Analysis by Limma

Next, we want to find genes, which are differentially expressed between BM-MOs, FL-MOs and YS-Macs using Limma. First, we perform a log transformation of the data.

# transform into log2 scale
ex <- exprs(eset)
exprs(eset) <- log2(ex)

Next, we improve the class labels from the 3 distinct studies and define a design matrices, which indicates conditions to be compared in the DE analysis.

# set up the data and proceed with analysis
eset$description <- sel_groups # Replace the long description by short labels
# Construct design matrices by model.matrix
# model.matrix loads the columns of eset and matches to its description. 
design <- model.matrix(~ description + 0, eset) # 0 defines the way we compare the samples
colnames(design) <- levels(sel_groups)

This is how the design table and array labels will look like.

eset$description
##  [1] BM_MOs  BM_MOs  BM_MOs  BM_MOs  FL_MOs  FL_MOs  FL_MOs  FL_MOs  YS_Macs
## [10] YS_Macs YS_Macs YS_Macs
## Levels: BM_MOs FL_MOs YS_Macs
design
##            BM_MOs FL_MOs YS_Macs
## GSM2042244      1      0       0
## GSM2042245      1      0       0
## GSM2042246      1      0       0
## GSM2042247      1      0       0
## GSM2042248      0      1       0
## GSM2042249      0      1       0
## GSM2042250      0      1       0
## GSM2042251      0      1       0
## GSM2042252      0      0       1
## GSM2042253      0      0       1
## GSM2042254      0      0       1
## GSM2042255      0      0       1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$description
## [1] "contr.treatment"

Finally, we will fit a liner model, perform differential analysis and estimate p-values.

# Fit linear model for each gene given a series of arrays
fit <- lmFit(eset, design)
# Build comparison and compute the satistics
cont.matrix <- makeContrasts(BM_MOs-FL_MOs, BM_MOs-YS_Macs, FL_MOs-YS_Macs, levels=design)
# Apply the contrast matrix to the linear model
fit2 <- contrasts.fit(fit, cont.matrix)
# Generate the statistics
fit2 <- eBayes(fit2, 0.05)

The makecontrast function indicates which pairs of conditions should be compared (BM_MOs vs FL_MOs, BM_MOs vs YS_Macs and FL_MOs vs.Ā YS_Macs). We can then use the function top table to obtain all DE genes for a given p-value (after FDR multiple test correction) for each of these 3 comparirsons.

# Extract the sig. DE genes
# coef number represents the order of column in cont.matrix:
# coef=1  BM_MOs-FL_MOs
# coef=2  BM_MOs-YS_Macs
# coef=3  FL_MOs-YS_Macs
de_genes_BM_MOs_FL_MOs <- topTable(fit2, coef=1, adjust="fdr", p.value=0.05, number=Inf)
de_genes_BM_MOs_YS_Macs <- topTable(fit2, coef=2, adjust="fdr", p.value=0.05, number=Inf)
de_genes_FL_MOs_YS_Macs <- topTable(fit2, coef=3, adjust="fdr", p.value=0.05, number=Inf)

head(de_genes_BM_MOs_FL_MOs) # show top de genes in BM_MO vs FL_MO comparirson. 

# Get all the DE genes regardless the comparisons
all_de_genes <- topTable(fit2, adjust="fdr", p.value=0.05, number=Inf)
all_genes_BM_MOs_FL_MOs <- topTable(fit2, coef=1, adjust="fdr", p.value=1, number=Inf)
# These tables will be used later

Please take a look into the content of these DE gene lists. Because the table is too big to be shown here, you can open them in Rstudio or output them to excel tables:

# Output in an excel file
wb <- createWorkbook()
# Create each sheet
addWorksheet(wb, "de_genes_BM_MOs_FL_MOs")
addWorksheet(wb, "de_genes_BM_MOs_YS_Macs")
addWorksheet(wb, "de_genes_FL_MOs_YS_Macs")
# Write the data into each sheet
writeData(wb, "de_genes_BM_MOs_FL_MOs", de_genes_BM_MOs_FL_MOs, rowNames = FALSE)
writeData(wb, "de_genes_BM_MOs_YS_Macs", de_genes_BM_MOs_YS_Macs, rowNames = FALSE)
writeData(wb, "de_genes_FL_MOs_YS_Macs", de_genes_FL_MOs_YS_Macs, rowNames = FALSE)
saveWorkbook(wb, "DE_genes_GSE76999.xlsx", overwrite = TRUE)

# Output CSV files if necessary
# write.table(de_genes_BM_MOs_FL_MOs, file="de_genes_BM_MOs_FL_MOs.csv", row.names=F, sep="\t")
# write.table(de_genes_BM_MOs_YS_Macs, file="de_genes_BM_MOs_YS_Macs.csv", row.names=F, sep="\t")
# write.table(de_genes_FL_MOs_YS_Macs, file="de_genes_FL_MOs_YS_Macs.csv", row.names=F, sep="\t")

Please open the output table DE_genes_GSE76999.xlsx and try to understand each column.

Get gene annotation

You probalby find that result is not very useful because you cannot recognize any genes by their probe ID (array-specific). What we prefer is the the gene symbol, but it is not included in the data. You may think, ā€œcan we convert the probe ID to gene symbol?ā€ Yes, this is doable.

Almost every micro array has a particular R package for this mapping task. Here we use package mogene10sttranscriptcluster.db which corresponds the Affymetrix mouse genes 1.0 described in GSE76999.

# mogene10sttranscriptcluster.db package contains the annotation of Affymetrix mogene10
# Let's see what is inside this package
columns(mogene10sttranscriptcluster.db)
# Get all the probe ID in our ExpressionSet
ID <- featureNames(eset)
# Select the symbols according to the IDs we just got 
tmp <- select(mogene10sttranscriptcluster.db, ID, c("SYMBOL", "ENTREZID"))
tmp <- na.omit(tmp)
# Then we can insert a column for symbols which match the ID for each row
all_genes_BM_MOs_FL_MOs$ENTREZID <- tmp$ENTREZID[match(rownames(all_genes_BM_MOs_FL_MOs), tmp$PROBEID)]
all_de_genes$symbol <- tmp$SYMBOL[match(rownames(all_de_genes), tmp$PROBEID)]
all_de_genes <- na.omit(all_de_genes)
de_genes_BM_MOs_FL_MOs$ENTREZID <- tmp$ENTREZID[match(rownames(de_genes_BM_MOs_FL_MOs), tmp$PROBEID)]
de_genes_BM_MOs_YS_Macs$ENTREZID <- tmp$ENTREZID[match(rownames(de_genes_BM_MOs_YS_Macs), tmp$PROBEID)]
de_genes_FL_MOs_YS_Macs$ENTREZID <- tmp$ENTREZID[match(rownames(de_genes_FL_MOs_YS_Macs), tmp$PROBEID)]

What if I have data from a distnct array? For example, you find the array name as Affymetrix Human Gene 2.0 ST Array [transcript (gene) version] in the given GEO description. What you can do is google this name with R bioconductor to see what package is available for annotation conversion. Or, you can visit the webpage of Bioconductor: https://bioconductor.org/packages/release/BiocViews.html#___ChipManufacturer where the list of Packages for all ChipManufacturer is found.

MA plot

The MA plot visualizes the differences between measurements taken in two samples, by transforming the data into M (log ratio) and A (mean average) scales, then plotting these values. Here we take BM-MOs v.s. FL_MOs as an example.

# By Limma build-in function
ma <- fit2[,"BM_MOs - FL_MOs"]
plotMA(ma)

Volcano plot

A volcano plot is a type of scatter-plot that is used to quickly identify the significance of changes in large data sets. It plots p-values versus fold-change.

{volcanoplot(fit2, coeff=1, xlim=c(-2,2), col="gray")
title("BM_MOs - FL_MOs")}

You can also save any figure into PDF file by the codes below:

pdf("volcanoplot.pdf", width=4, height=4) # Save the plot into PDF format
{volcanoplot(fit2, coeff=1, xlim=c(-2,2), col="gray")
title("BM_MOs - FL_MOs")}
dev.off() # close the graphical device

4 PCA (Principal component analysis)

Principal component analysis (PCA) is a mathematical algorithm that reduces the dimensionality of the data while retaining most of the variation in the data set1. It accomplishes this reduction by identifying directions, called principal components, along which the variation in the data is maximal. By using a few components, each sample can be represented by relatively few numbers instead of by values for thousands of variables. Samples can then be plotted, making it possible to visually assess similarities and differences between samples and determine whether samples can be grouped (see Nature Biotechnology 26, 303 - 304 (2008) ā€œWhat is principal component analysis?ā€).

# Performs a principal components analysis on the given data matrix
pca <- prcomp(t(ex))
# Take PC1 and PC2 for the plot
{plot(pca$x[,1:2],col=sel_groups, pch=19)
# include a legend for points 
legend("topright", inset=.05, labels, pch=19, col=1:3, horiz=TRUE)}

5 Clustering by Heatmap

Another way to visualize the similarity and the difference of the samples is to perform hierarchical clustering and plot results as a heatmap.

First, we will create a color maps for generating values from blue (low expression) to red (high expression). Try out this command (display.brewer.all()) to check other pallets. We also restrict the matrix to only show differentially expressed genes.

# ex is the matrix for all the expression values
# all_de_genes is all the DE genes
# Let's filter ex by the list in all_de_genes
sel <- match(rownames(all_de_genes), rownames(ex)) # Get the index of DE genes in ex
sel_ex <- ex[sel,] # Select those rows only

Next, we create a few functions defining the use of Pearson correlation as a distance matrix and complete linkage for cluster.

clust <- function(x) hclust(x, method="complete") # perfrom complete linkage clustering
dist <- function(x) as.dist((1-cor(t(x)))/2)  # use the inverse of correlation as distance. 
hmcol <- rev(brewer.pal(11, "RdBu"))

# heatmap.2 is a powerful and complicated function, please read its help message
hp <- heatmap.2(sel_ex, # input matrix 
                col=hmcol, # define the customized colormap
                scale="row", # scale the values for each gene (row)
                labRow=F, # Not showing row labels
                density.info="none", # Not showing density info
                trace="none", # Not showing trace line
                margins = c(8,3), # modify the margin to fit the labels
                labCol = rep(labels,each=4), # add the column labels
                # define hierarchical clustering method
                hclust = clust,
                # Using correlation coefficient for distance function
                distfun = dist
                )

This is a pretty heamap, or? As you noticed, this involves calling a complex function with several parameters.Play around with parameters: changing the colormap, linkage function of the clustering and see what you get. Also notice that computing this takes a few minutes and obtaninng a nice looking heatmap can involve a lot of work on parameters.

Get genes from clustering

If we want to get the clusters from the dendrogram on the left hand side, we can use the script below:

# Extract the hierarchical cluster from heatmap to class "hclust"
hc <- as.hclust(hp$rowDendrogram)
# Cut the tree by the dedired number of groups
clusters <- cutree(hc, k=5)
# Converting the ID into gene symbol into a list. Please try to understand this code step by step...
clustered_genes <- list(all_de_genes$symbol[match(names(clusters)[clusters==1],all_de_genes$ID)],
                        all_de_genes$symbol[match(names(clusters)[clusters==2],all_de_genes$ID)],
                        all_de_genes$symbol[match(names(clusters)[clusters==3],all_de_genes$ID)],
                        all_de_genes$symbol[match(names(clusters)[clusters==4],all_de_genes$ID)],
                        all_de_genes$symbol[match(names(clusters)[clusters==5],all_de_genes$ID)])

# Alternative way to do it
clustered_genes <- list(all_de_genes[names(clusters[clusters==1]),]$symbol,
                        all_de_genes[names(clusters[clusters==2]),]$symbol,
                        all_de_genes[names(clusters[clusters==3]),]$symbol,
                        all_de_genes[names(clusters[clusters==4]),]$symbol,
                        all_de_genes[names(clusters[clusters==5]),]$symbol)
# Output in an excel file
wb <- createWorkbook()
addWorksheet(wb, "cluster 1")
addWorksheet(wb, "cluster 2")
addWorksheet(wb, "cluster 3")
addWorksheet(wb, "cluster 4")
addWorksheet(wb, "cluster 5")
writeData(wb, "cluster 1", clustered_genes[[1]], rowNames = FALSE)
writeData(wb, "cluster 2", clustered_genes[[2]], rowNames = FALSE)
writeData(wb, "cluster 3", clustered_genes[[3]], rowNames = FALSE)
writeData(wb, "cluster 4", clustered_genes[[4]], rowNames = FALSE)
writeData(wb, "cluster 5", clustered_genes[[5]], rowNames = FALSE)
saveWorkbook(wb, "gene_clusters.xlsx", overwrite = TRUE)

Exercise Idea: try to reimplement the previous code with a loop so that you can also evaluate results with distinct number of clusters.

6 GO (Gene Ontology) Enrichment Analysis

Next, we will do a GO analysis to check GO terms associated to each of the DE gene sets. Here we will look at genes up- or down-regulated in the comparirsosn BM-MOs (G0) v.s. YS-Macs (G2).

go_enrich_up <- gprofiler(as.vector(all_de_genes$symbol[all_de_genes$BM_MOs...YS_Macs > 0]),
                          organism="mmusculus")
go_enrich_down <- gprofiler(as.vector(all_de_genes$symbol[all_de_genes$BM_MOs...YS_Macs < 0]),
                            organism="mmusculus")

head(go_enrich_up)
##   query.number significant p.value term.size query.size overlap.size precision
## 1            1        TRUE 0.00736        34       1908           13     0.007
## 2            1        TRUE 0.00588        78       1908           21     0.011
## 3            1        TRUE 0.01550      2346       1908          263     0.138
## 4            1        TRUE 0.02830       433       1908           65     0.034
## 5            1        TRUE 0.03480       135       1908           28     0.015
## 6            1        TRUE 0.00861       126       1908           28     0.015
##   recall    term.id domain subgraph.number
## 1  0.382 GO:0032731     BP              22
## 2  0.269 GO:0050764     BP              79
## 3  0.112 GO:0051128     BP              65
## 4  0.150 GO:0071900     BP             129
## 5  0.207 GO:0032640     BP             149
## 6  0.222 GO:0001906     BP              39
##                                                term.name relative.depth
## 1   positive regulation of interleukin-1 beta production              1
## 2                             regulation of phagocytosis              1
## 3          regulation of cellular component organization              1
## 4 regulation of protein serine/threonine kinase activity              1
## 5                       tumor necrosis factor production              1
## 6                                           cell killing              1
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             intersection
## 1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            GSDMD,JAK2,CASP1,PYCARD,PANX1,CASP4,HK1,NOD1,TLR8,RIPK2,TRIM16,NLRP1A,NAIP5
## 2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 HCK,CNN2,PTEN,LBP,ABR,ATG3,C3,FGR,STAP1,ATG7,PYCARD,RAB27A,CD300A,ABCA7,SCARB1,ATG5,RAP1GAP,CD300LF,CD47,PRTN3,SIRPB1A
## 3 BCL11A,DEF8,UNC119,TIAM1,AXL,PNKP,APOC2,CDKN1B,KLF4,HCK,PIH1D1,ARHGEF18,CNN2,CHMP2B,MAP2K1,CD44,CRBN,TNFSF14,INPP5K,TINF2,RHOA,ARPC5,RALA,TRPM2,MET,FAF1,PTEN,MYL2,DGUOK,DLL1,PRUNE1,DSTN,LBP,LAMP2,CRLF3,ABR,RETREG3,SGK2,BAIAP2L2,TRPV2,FIS1,STYXL1,GLE1,PKIB,LIMS1,SGK1,APAF1,COBL,IRAK3,SREBF1,FKBP1B,HNF1B,ITGB3,NIN,VRK1,RPS6KA5,SEMA4D,SMN1,HEXB,ANXA7,PRKCD,RGCC,BNIP3L,BIN3,FBXO4,NUBP1,CLDN1,MEFV,ATG3,FGD4,APP,IGFBP6,PARP3,TWF2,RPS6KA2,TBC1D5,CRIPT,C3,BAMBI,RIOK3,DIAPH1,SMAD4,NEDD4L,ANXA1,FAS,SLK,HDGFL3,MOB2,LIMD1,SGK3,PIKFYVE,FN1,ITM2C,ACTR3,RGS2,RAB29,ADIPOR1,MYOG,CAPN2,VANGL2,RABGAP1L,VIM,ARPC5L,LCN2,OLFM1,TOR1A,GSN,NUSAP1,PTPN1,P2RY1,GPSM2,ZBTB7B,SEMA4A,LMO4,GTF2B,PDLIM5,MIER1,TNFRSF1B,TNFRSF8,PPT1,PINK1,ADGRB2,HGF,RSPO1,FGR,PADI2,NPPB,BST1,TBC1D1,STAP1,CXCL9,SCARB2,TMEM106B,OCM,ATG7,CLEC2I,DMPK,CORO1A,PYCARD,TRIM30A,CASK,MSN,FLNA,MTM1,ARHGAP6,G6PDX,IKBKB,NECTIN1,LDLR,RAB27A,ANXA2,TMEM30A,CCDC88A,RAC2,CDC14A,CHMP2A,PLCG2,MTMR3,INPP5J,CD300A,CNOT6L,ANKRD27,MYL12B,CEBPA,DPP4,VPS13C,BRSK1,ABCA7,PALM,SYT1,SLAIN2,LRRK2,CDC42EP3,CYLD,WNT4,INF2,DDHD1,SSH2,SCARB1,ATG5,F11R,FIG4,RAP1GAP2,SS18L1,TBC1D2,PTGER4,RHOU,STAT2,XAF1,APOBEC1,CD53,EML2,SH3KBP1,RAP1GAP,IRF8,NOP53,CDC42EP4,TSPO,AGRN,S100A10,ARSB,INPP5F,RELN,STN1,ID1,DOCK5,PRKCE,INSIG1,PLAUR,STXBP6,CDC42SE1,IRGM1,CD24A,ABHD17B,CD300LF,DCUN1D3,SORL1,FZD9,ARHGAP15,FZD4,LGALS3,PLCB1,CD177,ATAD2B,RAB4B,ADAM10,SERTAD3,CD47,FMNL1,S100A9,CAPG,ARHGAP24,PRTN3,BAK1,UNC13D,GDA,BCAS3,PTK2B,FDPS,H3F3A,GMFG,H2-K1,ITGB1BP1,NRG1,MAPK3,STK24,VRK2,SIVA1,MIR195A,BUB3,CTCFL,CSF2RB,TRIM40,H2-D1,SVIP,CEACAM1,ARFIP1,TMSB15B2,SIRPB1A,CDKN2D,CHMP1B,PDE2A,VPS28
## 4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            DAZAP2,TIAM1,PRKAR2B,CDKN1B,PIH1D1,PTPN6,MAP2K1,INPP5K,RHOA,PTEN,CD40,ACSL1,PKIB,IRAK3,SHC2,ERN1,PIK3R5,FGF10,PRKCD,LATS2,RGCC,FGD4,STK38,CD74,JAK2,PDCD4,RGS2,PTPRC,ATP2B4,VANGL2,PTPN1,PDCD10,PTPN22,MAP3K6,HGF,PRKAG2,FLT1,PYCARD,STK26,SYAP1,LTF,MST1,CNPPD1,MAP3K13,CCND3,CD300A,UVRAG,LRRK2,GADD45A,MAP4K1,THBS1,MDFIC,ITGA1,PIK3R6,CD24A,SORL1,FZD4,TLR6,RGS14,PTK2B,NRG1,MAPK3,STK24,CEACAM1,CDKN2D
## 5                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      AXL,CYBB,LBP,CD34,VSIR,IRAK3,APP,IFIH1,PTPN22,TNFRSF8,CLEC4A2,PYCARD,TRIM30A,LTF,CD300LD,LRRK2,RARA,NOD1,PTGER4,THBS1,RIPK2,TSPO,TICAM1,CD24A,CCR2,GPR18,SPN,BCL3
## 6                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                PTPN6,FCER2A,H2-M3,LYST,STX7,ELANE,STAT5B,TAP2,PTPRC,STAP1,KLRK1,CLEC2D,PGLYRP1,CTSC,CORO1A,RAB27A,CTSH,LTF,TAP1,CAMP,CTSG,PIK3R6,CEBPG,UNC13D,H2-K1,H2-D1,CEACAM1,H60B

Then, of course, we can save these tables into excel files.

# Output in an excel file
wb <- createWorkbook()
addWorksheet(wb, "GO_BM_MOs_YS_Macs_UP")
addWorksheet(wb, "GO_BM_MOs_YS_Macs_DOWN")
writeData(wb, "GO_BM_MOs_YS_Macs_UP", go_enrich_up, rowNames = FALSE)
writeData(wb, "GO_BM_MOs_YS_Macs_DOWN", go_enrich_down, rowNames = FALSE)
saveWorkbook(wb, "GO_analysis.xlsx", overwrite = TRUE)

7 GSEA

What should we do if we wish to test a specific (maybe self-defined) gene set, to see whether they are up-regulated or down-regulated compared to other genes in terms of differential expression?

Then here is the Gene set enrichment GSEA.

We assume that the first 20 genes are a set of genes interested. We use camera from limma to see whether it is going up or down.

# assume the set of first 20 genes is the interested gene set
interested_gene_set <- rownames(exprs(eset))[1:20]

# Do the GSEA on limma
res.gsea <- camera(exprs(eset), list(set1 = interested_gene_set), design)
print(res.gsea)
##      NGenes Direction    PValue
## set1     20        Up 0.2180146

We can also take use of previous found gene set from Database MSigDB. (see MSigDB homepageā€). Now we import the NABA geneset from MSigDB.

m_df = msigdbr(species = "Mus musculus", category = "C2", subcategory = "CP")

naba <- m_df[grep("NABA", m_df$gs_name), ]

naba.l <- list()
for(i in unique(naba$gs_name)){
  sub <- naba[naba$gs_name == i, ]
  sub <- as.data.frame(sub)
  naba.l[[i]] <- as.character(sub$gene_symbol)
}

# Now we show what genes are inside the NABA_CORE_MATRISOME geneset
naba.l$NABA_CORE_MATRISOME
##   [1] "Abi3bp"   "Acan"     "Adipoq"   "Aebp1"    "Agrn"     "Ambn"    
##   [7] "Amelx"    "Amelx"    "Aspn"     "Bcan"     "Bglap"    "Bglap2"  
##  [13] "Bglap3"   "Bgn"      "Bmper"    "Bsph1"    "Ccn1"     "Ccn2"    
##  [19] "Ccn3"     "Ccn4"     "Ccn5"     "Ccn6"     "Cdcp2"    "Chad"    
##  [25] "Chadl"    "Cilp"     "Cilp2"    "Coch"     "Col10a1"  "Col11a1" 
##  [31] "Col11a2"  "Col12a1"  "Col13a1"  "Col14a1"  "Col15a1"  "Col16a1" 
##  [37] "Col17a1"  "Col18a1"  "Col19a1"  "Col1a1"   "Col1a2"   "Col20a1" 
##  [43] "Col22a1"  "Col23a1"  "Col24a1"  "Col25a1"  "Col26a1"  "Col27a1" 
##  [49] "Col28a1"  "Col2a1"   "Col3a1"   "Col4a1"   "Col4a2"   "Col4a3"  
##  [55] "Col4a4"   "Col4a5"   "Col4a6"   "Col5a1"   "Col5a2"   "Col5a3"  
##  [61] "Col6a1"   "Col6a2"   "Col6a3"   "Col6a5"   "Col6a6"   "Col7a1"  
##  [67] "Col8a1"   "Col8a2"   "Col9a1"   "Col9a2"   "Col9a3"   "Colq"    
##  [73] "Comp"     "Creld1"   "Creld2"   "Crim1"    "Crispld1" "Crispld2"
##  [79] "Cthrc1"   "Dcn"      "Dmbt1"    "Dmp1"     "Dpt"      "Dspp"    
##  [85] "Ecm1"     "Ecm2"     "Edil3"    "Efemp1"   "Efemp2"   "Egflam"  
##  [91] "Eln"      "Emid1"    "Emilin1"  "Emilin2"  "Emilin3"  "Epyc"    
##  [97] "Esm1"     "Fbln1"    "Fbln2"    "Fbln5"    "Fbln7"    "Fbn1"    
## [103] "Fbn2"     "Fga"      "Fgb"      "Fgg"      "Fgl1"     "Fgl2"    
## [109] "Fmod"     "Fn1"      "Fndc1"    "Fndc7"    "Fndc8"    "Fras1"   
## [115] "Gas6"     "Gldn"     "Hapln1"   "Hapln2"   "Hapln3"   "Hapln4"  
## [121] "Hmcn1"    "Hmcn2"    "Hspg2"    "Ibsp"     "Igfals"   "Igfbp1"  
## [127] "Igfbp2"   "Igfbp3"   "Igfbp4"   "Igfbp5"   "Igfbp6"   "Igfbp7"  
## [133] "Igfbpl1"  "Igsf10"   "Impg1"    "Impg2"    "Ints14"   "Ints6l"  
## [139] "Kcp"      "Kera"     "Lama1"    "Lama2"    "Lama3"    "Lama4"   
## [145] "Lama5"    "Lamb1"    "Lamb2"    "Lamb3"    "Lamc1"    "Lamc2"   
## [151] "Lamc3"    "Lgi1"     "Lgi2"     "Lgi3"     "Lgi4"     "Lrg1"    
## [157] "Ltbp1"    "Ltbp2"    "Ltbp3"    "Ltbp4"    "Lum"      "Matn1"   
## [163] "Matn2"    "Matn3"    "Matn4"    "Mepe"     "Mfap1a"   "Mfap2"   
## [169] "Mfap3"    "Mfap4"    "Mfap5"    "Mfge8"    "Mgp"      "Mmrn1"   
## [175] "Mmrn2"    "Igsf10"   "Ncan"     "Ndnf"     "Nell1"    "Nell2"   
## [181] "Nid1"     "Nid2"     "Npnt"     "Ntn1"     "Ntn3"     "Ntn4"    
## [187] "Ntn5"     "Ntng1"    "Ntng2"    "Nyx"      "Ogn"      "Oit3"    
## [193] "Omd"      "Optc"     "Otog"     "Otol1"    "Papln"    "Pcolce"  
## [199] "Pcolce2"  "Podn"     "Podnl1"   "Zp3"      "Postn"    "Prelp"   
## [205] "Prg2"     "Prg3"     "Prg4"     "Pxdn"     "Pxdn"     "Reln"    
## [211] "Rspo1"    "Rspo2"    "Rspo3"    "Rspo4"    "Sbspon"   "Slit1"   
## [217] "Slit2"    "Slit3"    "Smoc1"    "Smoc2"    "Sned1"    "Sparc"   
## [223] "Sparcl1"  "Spock1"   "Spock2"   "Spock3"   "Spon1"    "Spon2"   
## [229] "Spp1"     "Srgn"     "Srpx"     "Srpx2"    "Sspo"     "Svep1"   
## [235] "Tecta"    "Tectb"    "Tgfbi"    "Thbs1"    "Thbs2"    "Thbs3"   
## [241] "Thbs4"    "Thsd4"    "Tinag"    "Tinagl1"  "Tnc"      "Tnfaip6" 
## [247] "Tnn"      "Tnr"      "Tnxb"     "Tsku"     "Tspear"   "Ush2a"   
## [253] "Vcan"     "Vit"      "Vtn"      "Vwa1"     "Vwa2"     "Vwa3a"   
## [259] "Vwa3b"    "Vwa5a"    "Vwa5b1"   "Vwa5b2"   "Vwa7"     "Vwce"    
## [265] "Vwde"     "Vwf"      "Zp1"      "Zp2"      "Zp3"      "Zpld1"

One quiz: >How to get the full geneset of HALLMARK_GLYCOLYSIS? >How to test all gene sets from MSigDB

8 Combining other dataset for PCA and clustering

Finally, as an advanced example of data combination, we have selected Monocyte and Macrophage populations from the ImmGenn data (https://www.immgen.org/).

For simplicity, we will load the pre-processed data and a table containg details on the experiments (pheno data). To make these data comparable, we need to remove the batch effects caused by different labs, by using a tool called combat. Here, we are particularly interested in contrasting the expression patterns of monocytes and macrophages from distinct orgarms.

# get imm gene data from saved object
ex_imun_gen <- dget("Imun_gen.ex.Rdump")
# get all pheno data from tab sep. file
phenoData_table <- read.table(file="pheno.table.csv", header=T, sep=",", quote="")

phenoData_table # see the content of the phenotype table.

# combine two expression matrices 
ex_combined.raw <- merge(ex_imun_gen, ex, by="row.names")  
rownames(ex_combined.raw) <- ex_combined.raw$"Row.names"
ex_combined.raw$"Row.names" <- NULL
ex_combined.raw <- as.matrix(ex_combined.raw)

PCA on merged data

Once the data has been combined, we can perform PCA and clustering analysis in both data sets. Here tissues are represented by distinct colors and cell types by distinct symbols.

pca_im <- prcomp(t(ex_combined.raw))  # perform PCA
exp_names <- paste0(phenoData_table$cell, "_", phenoData_table$Tissue) # combine cell and tissue data. 
# Take PC1 and PC2 for the plot
#par(mai = c(1,  1, .25, 2)) 
{plot(pca_im$x[,1:2], 
     pch = ifelse(phenoData_table$cell == "Macrophage", 16, 17),#c(16, 17)[as.numeric(phenoData_table$cell)],  # cell type with distinct symbols
     col = brewer.pal(length(unique(phenoData_table$Tissue)), "Paired"))# brewer.pal(10,"Set1")[as.numeric(phenoData_table$Tissue)])  # cell tissue with distinct colors
# Legend outside the plot on right
legend("top",             # Location of legend 
  xpd = TRUE,                       # Allow drawing outside plot area
  xjust = 0,                        # Left justify legend box on x
  yjust = 1,                        # Center legend box on y
  legend = unique(phenoData_table$Tissue),
  col = brewer.pal(length(unique(phenoData_table$Tissue)), "Paired"),      # Legend Element Colors
  pch = c(16,17),                      # Legend Element Styles
  title = "Tissue")                 # Legend Title

# Legend outside the plot on bottom
legend("bottom",            # Location of legend
  xpd = TRUE,                       # Allow drawing outside plot area
  xjust = 0,                        # Left justify legend box on x
  yjust = 0,                        # Center legend box on y
  legend = unique(phenoData_table$cell),
  col = c("black", gray(.2)),       # Legend Element Colors
  pch = c(16, 17),                  # Symbol styles for legend elements
  title = "Cell")                   # Legend Title
}

This PCA plot are clearly separatedy into 2 parts, each representing the two data sets. This is knows as batch effect. Because we combine two different datasets, their different background, condition, or precedure may result in this separation. In order to remove it, we use ComBat function from package sva to remove it.

modcombat = model.matrix(~1, data=phenoData_table)
ex_combined.batch_removed = ComBat(dat=ex_combined.raw, 
                                   batch=phenoData_table$batch, mod=modcombat, par.prior=TRUE)

Lets check how the PCA looks like after this.

pca_im <- prcomp(t(ex_combined.batch_removed)) # pca analysis of batch removed matrix

# Take PC1 and PC2 for the plot
{plot(pca_im$x[,1:2], 
     pch = ifelse(phenoData_table$cell == "Macrophage", 16, 17),#c(16, 17)[as.numeric(phenoData_table$cell)],  # cell type with distinct symbols
     col= brewer.pal(length(unique(phenoData_table$Tissue)), "Paired"))# brewer.pal(10,"Set1")[as.numeric(phenoData_table$Tissue)])  # cell tissue with distinct colors
# Legend outside the plot on right
legend("topleft",             # Location of legend 
  xpd = TRUE,                       # Allow drawing outside plot area
  xjust = 0,                        # Left justify legend box on x
  yjust = 1,                        # Center legend box on y
  legend = unique(phenoData_table$Tissue),
  col = brewer.pal(length(unique(phenoData_table$Tissue)), "Paired"),      # Legend Element Colors
  pch = c(16,17),                      # Legend Element Styles
  title = "Tissue")                 # Legend Title

# Legend outside the plot on bottom
legend("bottomright",            # Location of legend
  xpd = TRUE,                       # Allow drawing outside plot area
  xjust = 0,                        # Left justify legend box on x
  yjust = 0,                        # Center legend box on y
  legend = unique(phenoData_table$cell),
  col = c("black", gray(.2)),       # Legend Element Colors
  pch = c(16, 17),                  # Symbol styles for legend elements
  title = "Cell")                   # Legend Title
}

Now, we see that similar cells (and tissue of origin) cluster appart.

Clustering of merged data

Alternativelly, we can also cluster these gene expresison profiles.

# Select the de genes from the expression matrix after batch effect removal
sel <- match(rownames(all_de_genes), rownames(ex_combined.batch_removed))
sel <- sel[!is.na(sel)] # remove the NAs

clust <- function(x) hclust(x, method="complete") # perfrom complete linkage clustering
dist <- function(x) as.dist((1-cor(t(x)))/2)  # use the inverse of correlation as distance. 
hmcol <- rev(brewer.pal(11, "RdBu"))

# heatmap.2 is a powerful and complicated function, please read its help message
hp <- heatmap.2(ex_combined.batch_removed[sel,], # input matrix 
                col=hmcol, # define the customized colormap
                scale="row", # scale the values for each gene (row)
                labRow=F, # Not showing row labels
                density.info="none", # Not showing density info
                trace="none", # Not showing trace line
                margins = c(11,3), # modify the margin to fit the labels
                labCol=paste0(phenoData_table$cell, " / ", phenoData_table$Tissue), # add the column labels
                # define hierarchical clustering method
                hclust = clust,
                # Using correlation coefficient for distance function
                distfun = dist
                )