Preparing the gene expression matrix

library(airway)
library(DESeq2)
data(airway)

# import data to DESeq2 and variance stabilize
dset = DESeqDataSetFromMatrix(assay(airway),
    colData=as.data.frame(colData(airway)), design=~dex)
dset = estimateSizeFactors(dset)
dset = estimateDispersions(dset)
gene_expr = getVarianceStabilizedData(dset)

# annotate matrix with HGNC symbols
library(biomaRt)
mart = useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes = getBM(attributes = c("ensembl_gene_id","hgnc_symbol"),
              values=rownames(gene_expr), mart=mart)
matched = match(rownames(gene_expr), genes$ensembl_gene_id)
rownames(gene_expr) = genes$hgnc_symbol[matched]

Obtaining pathway scores

We can then use the progeny function to score the expression matrix. Note that we are scaling the pathway scores with respect to the controls only.

library(progeny)
pathways = progeny(gene_expr, scale=FALSE)
controls = airway$dex == "untrt"
ctl_mean = apply(pathways[controls,], 2, mean)
ctl_sd = apply(pathways[controls,], 2, sd)
pathways = t(apply(pathways, 1, function(x) x - ctl_mean))
pathways = apply(pathways, 1, function(x) x / ctl_sd)

Checking for differences between the groups

So now we might be interested how the treatment with dexamethasone affects signaling pathways. To do this, we check if the control is different to the perturbed condition using a linear model:

library(dplyr)
result = apply(pathways, 1, function(x) {
    broom::tidy(lm(x ~ !controls)) %>%
        filter(term == "!controlsTRUE") %>%
        select(-term)
})
mutate(bind_rows(result), pathway=names(result))

## # A tibble: 11 x 5
##    estimate std.error statistic p.value pathway 
##       <dbl>     <dbl>     <dbl>   <dbl> <chr>   
##  1    3.36      5.74      0.586 0.579   EGFR    
##  2    0.966     1.92      0.504 0.633   Hypoxia 
##  3   -1.02      0.651    -1.57  0.167   JAK.STAT
##  4    2.25      1.45      1.56  0.170   MAPK    
##  5    0.354     0.792     0.447 0.670   NFkB    
##  6   -3.61      3.54     -1.02  0.346   PI3K    
##  7    1.41      0.980     1.44  0.199   TGFb    
##  8    0.841     0.685     1.23  0.265   TNFa    
##  9    1.56      0.897     1.74  0.132   Trail   
## 10    0.584     0.835     0.700 0.510   VEGF    
## 11   -4.34      0.756    -5.74  0.00122 p53

What we see is that indeed the p53/DNA damage response pathway is less active after treatment than before.

Getting the data

This example shows how to use the GDSC gene expression data of multiple cell lines together with PROGENy to calculate pathway activity and then to check for associations with drug sensitivity.

First, we need the GDSC data for both gene expression and drug response. They are available on the GDSC1000 web site:

# set up a file cache so we download only once
library(BiocFileCache)
bfc = BiocFileCache(".")
# gene expression and drug response
base = "http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources/Data/"
paths = bfcrpath(bfc, paste0(base, c("suppData/TableS4A.xlsx",
            "preprocessed/Cell_line_RMA_proc_basalExp.txt.zip")))

You can also download the files manually (adjust the file names when loading):

• Processed gene expression matrix
• Drug sensitivities

Creating the right objects to work with

First, we need to load the files we just downloaded into R to be able to perform the analysis:

# load the downloaded files
drug_table = readxl::read_excel(paths[1], skip=5)
gene_table = readr::read_tsv(paths[2])

# we need drug response with COSMIC IDs
drug_response = data.matrix(drug_table[,3:ncol(drug_table)])
rownames(drug_response) = drug_table[[1]]

# we need genes in rows and samples in columns
gene_expr = data.matrix(gene_table[,3:ncol(gene_table)])
colnames(gene_expr) = sub("DATA.", "", colnames(gene_expr), fixed=TRUE)
rownames(gene_expr) = gene_table$GENE_SYMBOLS

Running PROGENy to get pathway activity scores

Activity inference is done using a weighted sum of the model genes. We can run this without worrying about the order of genes in the expression matrix using:

library(progeny)
pathways = progeny(gene_expr)

We now have the pathway activity scores for the pathways defined in PROGENy:

head(pathways)

##                EGFR     Hypoxia    JAK.STAT       MAPK       NFkB
## 906826   0.03030286 -0.09136142 -0.36490995 -0.1758001 -0.5793367
## 687983  -0.99125434 -1.32673898 -0.93152060 -0.4946866 -1.3799417
## 910927  -0.10673190 -0.78816420 -1.06002081  0.1370551 -0.5497209
## 1240138 -0.05592591 -0.74266270 -0.07989446 -0.8259452  0.3418629
## 1240139 -0.15157011  0.11136425 -0.58596025 -0.2583581 -0.7256043
## 906792   0.71386069  0.39667896 -0.50001888  1.1967197 -0.4005830
##               PI3K       TGFb       TNFa      Trail        VEGF        p53
## 906826  -0.1999210 -0.6198524 -0.4724567 -0.5891909  0.18688452 -1.1725585
## 687983   0.3824370 -0.6696468 -1.0229424 -0.6113840 -0.06262960 -1.0818725
## 910927  -0.2155790  0.6214328 -0.1737935 -0.9185408  0.24335159  0.8249120
## 1240138  0.5883394  1.8891349  1.0191163  0.1214765 -0.15953605  2.1774919
## 1240139  1.0191110  0.9312615 -0.4347272 -0.2985134  0.36720972  0.8348820
## 906792  -2.1897400 -0.3093659 -0.1523604 -0.1621503  0.08751554 -0.9558531

Testing if MAPK activity is significantly associated with Trametinib

Trametinib is a MEK inhibitor, so we would assume that cell lines that have a higher MAPK activity are more sensitive to MEK inhibition.

We can test this the following way:

cell_lines = intersect(rownames(pathways), rownames(drug_response))
trametinib = drug_response[cell_lines, "Trametinib"]
mapk = pathways[cell_lines, "MAPK"]

associations = lm(trametinib ~ mapk)
summary(associations)

## 
## Call:
## lm(formula = trametinib ~ mapk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9965 -1.5286  0.3535  1.5446  6.8271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.03670    0.07155  -14.49   <2e-16 ***
## mapk        -1.31733    0.07095  -18.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.125 on 880 degrees of freedom
##   (80 observations deleted due to missingness)
## Multiple R-squared:  0.2815, Adjusted R-squared:  0.2806 
## F-statistic: 344.7 on 1 and 880 DF,  p-value: < 2.2e-16

And indeed we find that MAPK activity is strongly associated with sensitivity to Trametinib: the Pr(>|t|) is much smaller than the conventional threshold of 0.05.

The intercept is significant as well, but we’re not really interested if the mean drug response is above or below 0 in this case.

Note, however, that we tested all cell lines at once and did not adjust for the effect different tissues may have.