In greenelab/ADAGEpath: ADAGE-based Signature Analysis

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

This is an example for analyzing a publicly available microarray dataset from the ArrayExpress database.

Data preparation

Load in required libraries.

library("ADAGEpath")
library("DT")
library("dplyr")

Before any analysis, we need to specify the ADAGE model and the data compendium we want to use.

model <- eADAGEmodel
compendium <- PAcompendium
probe_dist <- probedistribution

Here we use the dataset E-GEOD-41926 as an example. Modify the accession number to load other P.a. microarray dataset from ArrayExpress.

accession <- "E-GEOD-41926"
data_raw <- load_dataset(input = accession, isProcessed = FALSE,
                         isRNAseq = FALSE, model = model,
                         compendium = compendium,
                         quantile_ref = probe_dist,
                         download_folder = "./download", norm01 = FALSE)

ADAGE only accepts expression values in the (0,1) range. We linearly transform expression values to be between 0 and 1 using the Pa compendium as the reference.

data_normed <- zeroone_norm(input_data = data_raw, use_ref = TRUE,
                            ref_data = compendium)

Now we query the sample information from ArrayExpress to understand the experimental condition of each sample.

pheno_table <- get_sample_info(accession)

We reorder samples in the pheno_table to have the same order as the expression data.

pheno_table <- pheno_table[match(colnames(data_raw)[-1],
                                 pheno_table$`Array Data File`), ]
DT::datatable(pheno_table, class = 'cell-border stripe')

We set the sample phenotypes according to the genotype information.

data_pheno <- pheno_table$`Characteristics [genotype]`

ADAGE signature analysis

Activity calculation

We calculate the activity of each signature for each sample in the dataset.

data_activity <- calculate_activity(input_data = data_normed, model = model)

The returned data_activity is a data.frame with signature names in the first column and activity values per sample starting from the second column.

Active signature detection

We are interested in identifying differentially active signatures between wildtype and plcH mutant.

We set indices to samples that will be included in the two-group comparison. Here the first four samples are two wildtype and two plcH mutant. Modify the indices to include different samples in the two-group comparison.

indices <- 1:4

We use limma to perform a differential activation test. limma is more robust than a simple t test when sample size is small. A two-group limma analysis is provided in the function build_limma(). You can also build other limma models to test signatures' activities when the experimental design is more complex.

limma_result <- build_limma(input_data = data_activity[, c(1, indices + 1)],
                            phenotypes = data_pheno[indices],
                            use.bonferroni = FALSE)

To take both absolute activity difference and significance into account, we use pareto fronts to pick the most differentially active signatures. We extract differentially active signatures in the first 10 layers of pareto fronts. Modify N_fronts to get more or fewer signatures.

active_sigs <- get_active_signatures(limma_result = limma_result,
                                     pheno_group = "both",
                                     method = "pareto", N_fronts = 10)

Signatures that are differentially active between wildtype and plcH mutant are:

active_sigs

Plot each signature's activity changes and significance in the limma test.

plot_volcano(limma_result = limma_result, highlight_signatures = active_sigs,
             interactive = TRUE)

Look at how the activities of active signature vary across samples.

plot_activity_heatmap(activity = data_activity, signatures = active_sigs)

Based on the volcano plot and the activity heatmap, we can see that Node34pos, Node28neg, Node41pos, Node50pos, Node69pos, Node49neg, and Node65pos are active in wildtype, while the rest signatures become more active in plcH mutant.

Overlapping signature removal

To reduce the number of signatures to look at, we can check whether these active signatures greatly overlap with others.

plot_signature_overlap creates a heatmap of odds ratios. The odds ratio represents the odds that two signatures share a specific number of genes.

signature_similarity <- plot_signature_overlap(selected_signatures = active_sigs,
                                               model = model)

We can see that the four signatures (Node50pos, Node28neg, Node34pos, and Node41pos) that are most strongly active in wildtype share significant number of genes. Node49neg overlap with Node65pos. Node34neg and Node28pos that are active in plcH mutant also overlap.

Next we calculate the marginal activities of active signatures. Marginal activity is defined as the activity of signature A after removing genes that are shared with signature B.

marginal_activity <- calculate_marginal_activity(
  input_data = data_normed,
  selected_signatures = active_sigs, model = model)

Again, we build a limma model to test whether these marginal activities are still strongly different between two conditions.

marginal_limma <- build_limma(input_data = marginal_activity[, c(1, indices + 1)],
                              phenotypes = data_pheno[indices])

Let's visualize the marginal activties in a matrix heatmap. The value in this matrix represents the -log10 transformed adjusted p value in the activation test when the effect of the column signature is removed from the row signature. Values in the diagonal of the heatmap are the activation significance of signatures themselves. Activation significance below the cutoff is marked by a cross sign.

plot_marginal_activation(marginal_limma_result = marginal_limma,
                         signature_order = colnames(signature_similarity),
                         sig_cutoff = 0.05)

Node58neg, Node116pos, Node272pos are masked by Node257pos. Node68neg is masked by Node265neg. Node 34neg and Node28pos can mask each other. We keep Node34neg because it has higher significance. Node50pos, Node34pos, Node41pos, and Node49neg are masked by Node28neg.

unique_active_sigs <- remove_redundant_signatures(marginal_limma,
                                                  sig_cutoff = 0.05)
unique_active_sigs

Active signature interpretation

Associated pathways

We download Pseudomonas aeruginosa KEGG pathway terms from the TRIBE web server.

KEGG <- fetch_geneset(type = "KEGG")
# we only consider KEGG pathways with more than 5 genes and less than 100 genes
# as meaningful pathways
KEGG_subset <- KEGG[lengths(KEGG) >= 5 & lengths(KEGG) <= 100]

We associate active signatures to known KEGG pathways.

pathway_association <- annotate_signatures_with_genesets(
  selected_signatures = unique_active_sigs, model = model, genesets = KEGG_subset)

Calculate the activity of associated pathways inside active signatures. This help differentiate active vs. non-active pathways assoicated with active signatures.

pathway_activity <- signature_geneset_activity(
  signature_geneset_df = pathway_association[, c("signature", "geneset")],
  gene_set_list = KEGG_subset, model = model, input_data = data_normed)
plot_activity_heatmap(pathway_activity)

We run a limma test on pathway activities and find pathways that are truly active in this dataset.

pathway_limma <- build_limma(pathway_activity[, c(1, indices + 1)],
                             phenotypes = data_pheno[indices])
# combine pathway association and pathway activation test results
combined_result <- combine_geneset_outputs(
  signature_geneset_association = pathway_association,
  geneset_limma_result = pathway_limma)
knitr::kable(combined_result, digits = 4, row.names = FALSE, align = "c")

There are also signatures uncharacterized by KEGG.

uncharacterized_sigs <- setdiff(unique_active_sigs, pathway_association$signature)
uncharacterized_sigs

Check the signature similarity of the uncharacterized signatures.

plot_signature_overlap(selected_signatures = uncharacterized_sigs, model = model)

Gene-gene network

We check how genes in the active signatures cluster in the ADAGE gene-gene network.

We calculate an expression fold change for each gene and pass it to the gene-gene network to show as node color. Again, we use limma to test differential expression and get the logFC.

data_raw_limma <- build_limma(input_data = data_raw[, c(1, indices + 1)],
                              phenotypes = data_pheno[indices],
                              use.bonferroni = FALSE)
# build a gene:fold change table from limma result
gene_logFC <- data.frame(geneID = rownames(data_raw_limma),
                         logFC = data_raw_limma$logFC)

Visualize the ADAGE gene-gene network of the active signatures.

visualize_gene_network(selected_signatures = unique_active_sigs,
                       model = model, cor_cutoff = 0.5,
                       gene_color_value = gene_logFC,
                       curated_pathways = KEGG)

Genes in signatures

To review a signature, we can list all genes in it.

gene_annotation <- annotate_genes_in_signatures(selected_signatures = "Node28neg",
                                                model = model,
                                                curated_pathways = KEGG)
# add expression fold change to the gene table
gene_annotation <- dplyr::right_join(gene_logFC, gene_annotation,
                                     by = c("geneID" = "LocusTag"))
DT::datatable(gene_annotation)