Using rcellminer

In rcellminer: rcellminer: Molecular Profiles, Drug Response, and Chemical Structures for the NCI-60 Cell Lines

library(knitr)
opts_chunk$set(out.extra='style="display:block; margin: auto"', 
                             fig.align="center", fig.width=8, fig.height=8, tidy=FALSE)

Overview

The NCI-60 cancer cell line panel has been used over the course of several decades as an anti-cancer drug screen. This panel was developed as part of the Developmental Therapeutics Program (DTP, http://dtp.nci.nih.gov/) of the U.S. National Cancer Institute (NCI). Thousands of compounds have been tested on the NCI-60, which have been extensively characterized by many platforms for gene and protein expression, copy number, mutation, and others (Reinhold, et al., 2012). The purpose of the CellMiner project (http://discover.nci.nih.gov/cellminer) has been to integrate data from multiple platforms used to analyze the NCI-60, and to provide a powerful suite of tools for exploration of NCI-60 data. While CellMiner is an unmatched resource for online exploration of the NCI-60 data, consideration of more specialized scientific questions often requires custom programming. The rcellminer R package complements the functionality of CellMiner, providing programmatic data access, together with functions for data visualization and analysis. These functions are approachable for even beginning R users, as illustrated by the initial examples below. The subsequent case studies, inspired by CellMiner-related publications, show how modest amounts of code can script specialized analyses, integrating multiple types of data to yield new scientific insights. rcellminer functions also provide robust building blocks for more extensive tools, as exemplifed by the package's interactive Shiny applications.

Basics

Installation

if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("rcellminer")
BiocManager::install("rcellminerData")

Getting Started

Load rcellminer and rcellminerData packages:

library(rcellminer)
library(rcellminerData)

A list of all accessible vignettes and methods is available with the following command.

help.search("rcellminer")

Searching for Compounds

The NSC number is a numeric identifier for substances submitted to the National Cancer Institute (NCI) for testing and evaluation. It is a registration number for the Developmental Therapeutics Program (DTP) repository, and it is used as the unique identifier for compounds in the CellMiner database. NSC stands for National Service Center.

rcellminer allows users to quickly search for NSC IDs by compound name or partial name. For example, many kinase inhibitors end with the suffix "nib". Users can quickly search NSCs for compound names with this suffix; queries are case insensitive and are treated as regular expressions.

searchForNscs("nib$")  

Profile Visualization

Often, it is useful for researchers to plot multiple data profiles next to each other in order to visually identify patterns. Below are examples for the visualization of various profiles: single drugs and multiple drugs, as well as molecular profiles and combinations of drug and molecular profiles.

# Get Cellminer data
drugAct <- exprs(getAct(rcellminerData::drugData))
molData <- getMolDataMatrices()

# One drug
nsc <- "94600"
plots <- c("drug") 
plotCellMiner(drugAct, molData, plots, nsc, NULL)

# One expression
gene <- "TP53"
plots <- c("exp") 
plotCellMiner(drugAct, molData, plots, NULL, gene)

# Two drugs: Irinotecan and topotecan 
nsc <- c("616348", "609699")
plots <- c("drug", "drug") 
plotCellMiner(drugAct, molData, plots, nsc, NULL)

# Two genes 
# NOTE: subscript out of bounds Errors likely mean the gene is not present for that data type
gene <- c("TP53", "MDM2")
plots <- c("exp", "mut", "exp") 
plotCellMiner(drugAct, molData, plots, NULL, gene)

# Gene and drug to plot 
nsc <- "609699"
gene <- "SLFN11"
plots <- c("exp", "cop", "drug") 
plotCellMiner(drugAct, molData, plots, nsc, gene)

Visualizing Drug Sets

For related drugs, it is often useful to visualize a summary of drug activity. rcellminer allows you to easily plot a drug set's average activity z-score pattern over the NCI-60, together with 1 standard deviation error bars to gauge variation in response behavior.

# Get CellMiner data
drugAct <- exprs(getAct(rcellminerData::drugData))

drugs <- searchForNscs("camptothecin")
drugAct <- drugAct[drugs,]
mainLabel <- paste("Drug Set: Camptothecin Derivatives, Drugs:", length(drugs), sep=" ")

plotDrugSets(drugAct, drugs, mainLabel) 

Structure Visualization

Drug compound structures are provided as annotations for drug compounds. The structures of CellMiner compounds can be visualized and compared using the rcdk package.

drugAnnot <- getFeatureAnnot(rcellminerData::drugData)[["drug"]]
drugAnnot[1:5, "SMILES"]

Working with Additional Drug Information

Mechanism of action (MOA)

rcellminer provides the mechanism of action (MOA) category for compounds for which this information is available. The MOA specifies the biochemical targets or processes that a compound affects. MOAs are indicated by CellMiner abbreviations, though more detailed descriptions can be obtained as indicated below.

Get MOA information

Find known MOA drugs and organize their essential information in a table.

# getFeatureAnnot() returns a named list of data frames with annotation data 
# for drugs ("drug") and drug activity repeat experiments ("drugRepeat").
drugAnnot <- getFeatureAnnot(rcellminerData::drugData)[["drug"]]

# Get the names and MOA for compounds with MOA entries
knownMoaDrugs <- unique(c(getMoaToCompounds(), recursive = TRUE))
knownMoaDrugInfo <- data.frame(NSC=knownMoaDrugs, stringsAsFactors = FALSE)
knownMoaDrugInfo$Name <- drugAnnot[knownMoaDrugInfo$NSC, "NAME"]

# Process all NSCs 
knownMoaDrugInfo$MOA <- sapply(knownMoaDrugInfo$NSC, getMoaStr)

# Order drugs by mechanism of action.
knownMoaDrugInfo <- knownMoaDrugInfo[order(knownMoaDrugInfo$MOA), ]

# Drug_MOA_Key data frame provides further details on MOA abbrevations.
Drug_MOA_Key[c("A2", "A7"), ]

Drug Activity

Additionally, rcellminer provides GI50 (50% growth inhibition) values for the compounds in the database. GI50 values are similar to IC50 values, which are the concentrations that cause 50% growth inhibition, but have been renamed to emphasize the correction for the cell count at time zero. Further details on the assay used are available on the DTP website.

Get GI50 values

Compute GI50 data matrix for known MOA drugs.

negLogGI50Data <- getDrugActivityData(nscSet = knownMoaDrugInfo$NSC)

gi50Data <- 10^(-negLogGI50Data)

Construct integrated data table (drug information and NCI-60 GI50 activity).

knownMoaDrugAct <- as.data.frame(cbind(knownMoaDrugInfo, gi50Data), stringsAsFactors = FALSE)

# This table can be written out to a file
#write.table(knownMoaDrugAct, file="knownMoaDrugAct.txt", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE, na="NA")

Use Cases

Pattern Comparison

In this example, we construct and visualize a composite gene inactivation pattern for the widely studied PTEN tumor suppressor gene, integrating gene expression, copy, and mutation data. The composite pattern is then used to identify compounds with correlated activity patterns.

# Gene of interest
gene <- "PTEN"

# Get Data
drugAct <- exprs(getAct(rcellminerData::drugData))
molData <- getMolDataMatrices()

# Get the cell lines names for cell lines meeting particular thresholds
copKnockdown <- names(which(molData[["cop"]][paste0("cop", gene), ] < -1))
expKnockdown <- names(which(molData[["exp"]][paste0("exp", gene), ] < -1.5))
mutKnockdown <- names(which(molData[["mut"]][paste0("mut", gene), ] == 1))

# Make composite pattern
pattern <- rep(0, length(molData[["cop"]][paste0("cop", gene), ]))
names(pattern) <- names(molData[["cop"]][paste0("cop", gene), ])
tmp <- Reduce(union, list(copKnockdown, expKnockdown, mutKnockdown))
pattern[tmp] <- 1

# Composite plot data
extraPlot <- list()
extraPlot[["title"]] <- "Composite Pattern"
extraPlot[["label"]] <- "Knockdown Composite (Binary)"
extraPlot[["values"]] <- pattern

# Plot data
plotCellMiner(molData=molData, plots=c("cop", "exp", "mut"), gene=gene, extraPlot=extraPlot)

# Significant drug correlations to FDA-approved/Clinical trial drugs
# getFeatureAnnot() returns a named list of data frames with annotation data 
# for drugs ("drug") and drug activity repeat experiments ("drugRepeat").
drugAnnot <- getFeatureAnnot(rcellminerData::drugData)[["drug"]]
tmpDA <- drugAnnot[drugAnnot$FDA_STATUS != "-", c("NSC", "FDA_STATUS")]

tmpDrugAct <- drugAct[rownames(drugAct) %in% tmpDA$NSC,]

results <- patternComparison(pattern, tmpDrugAct)
sigResults <- results[results$PVAL < 0.01, ]

# Get names of drugs that are FDA approved or in clinical trials
getDrugName(rownames(sigResults))

Correlating DNA Copy Number Alteration and Gene Expression

We can assess the extent to which a gene's expression may be regulated by copy gain or loss by computing the correlation between its NCI-60 transcript expression and DNA copy number alteration pattern. The overall distribution of such correlations can guide consideration of their significance. In this setting, it is insightful to consider the extent to which particular classes of genes may be regulated by copy alteration. This example illustrates the computation and visualization of the gene copy number vs. expression correlation distribution, while situating the median correlation for a set of oncogenes within the observed distribution.

molData <- getMolDataMatrices()

# Get Data
copGenes <- removeMolDataType(rownames(molData[["cop"]]))
expGenes <- removeMolDataType(rownames(molData[["exp"]]))
genes <- intersect(copGenes, expGenes)

# Generate the appropriate rownames 
expGeneLabels <- paste0("exp", genes)
copGeneLabels <- paste0("cop", genes)

a <- molData[["exp"]][expGeneLabels,]
b <- molData[["cop"]][copGeneLabels,]
allGenes <- rowCors(a, b)

#selectedOncogenes <- c("FZD1", "JAK2", "ALK", "PIK3CG", "RET", "CDC25A", "PDGFB", "PIK3CB", "WNT3")
selectedOncogenes <- c("ABL1", "ALK", "BRAF", "CCND1", "CCND3", "CCNE1", "CCNE2", 
                                             "CDC25A", "EGFR", "ERBB2", "EZH2", "FOS", "FOXL2", "HRAS", 
                                             "IDH1", "IDH2", "JAK2", "KIT", "KRAS", "MET", "MOS", "MYC", 
                                             "NRAS", "PDGFB", "PDGFRA", "PIK3CA", "PIK3CB", "PIK3CD", 
                                             "PIK3CG", "PIM1", "PML", "RAF1", "REL", "RET", "SRC", "STK11", 
                                             "TP63", "WNT10B", "WNT4", "WNT2B", "WNT9A", "WNT3", "WNT5A", 
                                             "WNT5B", "WNT10A", "WNT11", "WNT2", "WNT1", "WNT7B", "WISP1", 
                                             "WNT8B", "WNT7A", "WNT16", "WISP2", "WISP3", "FZD5", "FZD1")
selectedOncogenes <- intersect(selectedOncogenes, genes)

# Generate the appropriate rownames 
expGeneLabels <- paste0("exp", selectedOncogenes)
copGeneLabels <- paste0("cop", selectedOncogenes)

a <- molData[["exp"]][expGeneLabels,]
b <- molData[["cop"]][copGeneLabels,]
selectedOncogenesCor <- rowCors(a, b)

hist(allGenes$cor, main="", xlab="Pearson correlation between expression and copy number", breaks=200, col="lightgray", border="lightgray")

segments(x0=median(allGenes$cor), y0=0, x1=median(allGenes$cor), y1=175, lty=2, lwd=2)
text(median(allGenes$cor)+0.02, y=175, adj=0, labels="Median Correlation\nAll Genes", cex=0.75)

segments(x0=median(selectedOncogenesCor$cor), y0=0, x1=median(selectedOncogenesCor$cor), y1=140, lty=2, lwd=2, col="red")
text(median(selectedOncogenesCor$cor)+0.02, y=140, adj=0, labels="Median Correlation\nOncogenes", cex=0.75)

rug(selectedOncogenesCor$cor, col="red") 

Assessing Correlative Evidence for a Drug MOA Hypothesis

rcellminer makes it easy to consider whether the available pharmacogenomic data support a hypothesis relating to drug mechanism of action. For example, higher expression of anti-apoptotic genes could plausibly promote cell survival in the face of DNA damage induced by topoisomerase inhibitors, thereby contributing to drug resistance. To explore this further, we can check if the expression of known anti-apoptotic genes is significantly negatively correlated with the activity of the topoisomerase 1 inhibitor camptothecin.

# Get normalized (Z-score) NCI-60 gene expression and drug activity data.
nci60DrugActZ <- exprs(getAct(rcellminerData::drugData))
nci60GeneExpZ <- getAllFeatureData(rcellminerData::molData)[["exp"]]

antiApoptosisGenes <- c("BAG1", "BAG3", "BAG4", "BCL10", "BCL2", 
                                                "BCL2A1", "BCL2L1", "BCL2L10", "BCL2L2", 
                                                "BFAR", "BIRC3", "BIRC6", "BNIP1", "BNIP2", 
                                                "BNIP3", "BNIP3L", "BRAF", "CASP3", "CD27", 
                                                "CD40LG", "CFLAR", "CIDEA", "DAPK1", "DFFA",
                                                "FAS", "IGF1R", "MCL1", "NOL3", "TMBIM1", 
                                                "TP53", "TP73", "XIAP")
camptothecinNsc <- "94600"

# Compute table of correlations between camptothecin activity and anti-apoptosis gene expression.
pcResults <- patternComparison(nci60DrugActZ[camptothecinNsc, ], nci60GeneExpZ[antiApoptosisGenes, ])

# Adjust p-values for multiple comparisons, sort with respect to q-values.
pcResults$QVAL <- p.adjust(pcResults$PVAL, method = "fdr")
pcResults <- pcResults[order(pcResults$QVAL), ]

# Identify genes with significant negative correlations (FDR < 0.1)
pcResults[((pcResults$COR < 0) & (pcResults$QVAL < 0.10)), ]

Two genes, the BAX-antagonists TMBIM1 and BCL2L2, have NCI-60 expression patterns that are significantly negatively correlated with camptothecin activity. Plots illustrating these relationships are constructed below.

colorTab <- loadNciColorSet(returnDf = TRUE)
tissueColorTab <- unique(colorTab[, c("tissues", "colors")])

plotData <- as.data.frame(t(rbind(nci60DrugActZ[camptothecinNsc, , drop=FALSE], 
                                                                    nci60GeneExpZ[c("TMBIM1", "BCL2L2"), ])))
colnames(plotData) <- c("Camptothecin", "TMBIM1", "BCL2L2")

plot(x=plotData$TMBIM1, y=plotData$Camptothecin, col=colorTab$colors, pch=16,
         xlab="TMBIM1 mRNA exp (Z-score)", ylab="Camptothecin Activity (Z-score)",
         main=paste0("Camptothecin Activity vs. TMBIM Expression (r = ", 
                                round(pcResults["TMBIM1", "COR"], 2), ")"))
abline(lm(formula("Camptothecin ~ TMBIM1"), plotData), col="red")
legend("bottomleft", legend=tissueColorTab$tissues, col=tissueColorTab$colors, 
             cex=0.7, pch = 16)

plot(x=plotData$BCL2L2, y=plotData$Camptothecin, col=colorTab$colors, pch=16,
         xlab="BCL2L2 mRNA exp (Z-score)", ylab="Camptothecin Activity (Z-score)",
         main=paste0("Camptothecin Activity vs. BCL2L2 Expression (r = ", 
                                round(pcResults["BCL2L2", "COR"], 2), ")"))
abline(lm(formula("Camptothecin ~ BCL2L2"), plotData), col="red")
legend("bottomleft", legend=tissueColorTab$tissues, col=tissueColorTab$colors, 
             cex=0.7, pch = 16)

Relating Gene Alteration Patterns to Drug Responses

Alterations of selected DNA repair genes may have predictive value for anticancer drug activity. Using rcellminer, we can flexibly script analyses to explore such gene-drug associations. The analyses below follow elements of the approach presented in Sousa et al., though the detailed results may differ somewhat due to simplifications made for illustrative purposes. Our aim is to catalogue impaired DNA repair gene function in the NCI-60, due either to deleterious mutation or lack of expression. These alteration patterns can then be related to drug response behavior by considering differences in drug activity between altered and non-altered cell lines.

We start by obtaining NCI-60 drug activity data (-logGI50) for 158 FDA-approved drugs, together with corresponding gene expression data (log2 intensity), and binary muation data (where a value of one indicates the presence of a deleterious mution, in either the heterozygous or homozygous state).

# Get Cellminer data
# getFeatureAnnot() returns a named list of data frames with annotation data 
# for drugs ("drug") and drug activity repeat experiments ("drugRepeat").
drugAnnot <- getFeatureAnnot(rcellminerData::drugData)[["drug"]]
fdaDrugAnnot <- drugAnnot[which(drugAnnot$FDA_STATUS == "FDA approved"), ]

nci60FdaDrugAct <- getDrugActivityData(nscSet = fdaDrugAnnot$NSC)
nci60GeneExp <- getAllFeatureData(rcellminerData::molData)[["xai"]]
nci60GeneMut <- getAllFeatureData(rcellminerData::molData)[["mut"]]

We then visualize the NCI-60 mutation patterns for a subset of DNA repair genes.

dnarGeneSet <- c("APC", "APLF", "ATM", "CLSPN", "ERCC6", "FANCI", 
                 "FANCM", "GEN1", "HLTF", "MLH1", "POLD1", "POLE", "POLG",
                 "POLQ", "RAD54L", "REV3L", "RMI1", "SLX4", "SMARCA4", "SMC2",
                 "TP53", "WRN")
dnarGeneSet <- intersect(dnarGeneSet, rownames(nci60GeneMut))

dnarGeneMut <- nci60GeneMut[dnarGeneSet, ]

# Identify most frequently mutated genes.
numMutLines <- apply(dnarGeneMut, MARGIN = 1, FUN = sum)
sort(numMutLines, decreasing = TRUE)

dnarGeneMutPlotData <- dnarGeneMut[order(numMutLines), ]
heatmap(dnarGeneMutPlotData, scale="none", Rowv = NA, Colv = NA, col = c("grey", "red"), 
                main="Deleterious mutations")

A similar plot can be developed to indicate DNA repair genes with little or no expression in particular NCI-60 cell lines.

dnarGeneExp <- nci60GeneExp[dnarGeneSet, ]

hist(dnarGeneExp, xlab="mRNA expression (log2 intensity)",
         main="Distribution of DNA repair gene expression values")

# Set low expression threshold to 1st percentile of expression values.
lowExpThreshold <- quantile(dnarGeneExp, na.rm = TRUE, probs = c(0.01))
lowExpThreshold

# Construct binary (potential) expression knockout matrix.
dnarGeneExpKo <- matrix(0, nrow=nrow(dnarGeneExp), ncol=ncol(dnarGeneExp))
rownames(dnarGeneExpKo) <- rownames(dnarGeneExp)
colnames(dnarGeneExpKo) <- colnames(dnarGeneExp)
dnarGeneExpKo[which(dnarGeneExp < lowExpThreshold)] <- 1

# Restrict to wild type expression knockouts.
dnarGeneExpKo[which(dnarGeneMut == 1)] <- 0

# Identify genes with most frequent loss of expression.
numExpKoLines <- apply(dnarGeneExpKo, MARGIN = 1, FUN = sum)
sort(numExpKoLines, decreasing = TRUE)

dnarGeneExpKoPlotData <- dnarGeneExpKo[order(numExpKoLines), ]
heatmap(dnarGeneExpKoPlotData, scale="none", Rowv = NA, Colv = NA, 
                col = c("grey", "blue"), main="Potential expression knockouts")

Finally, we merge the data from the above plots to show DNA repair gene alterations by mutation or lack of expression.

dnarGeneAlt <- matrix(0, nrow=nrow(dnarGeneExp), ncol=ncol(dnarGeneExp))
rownames(dnarGeneAlt) <- rownames(dnarGeneExp)
colnames(dnarGeneAlt) <- colnames(dnarGeneExp)
dnarGeneAlt[which(dnarGeneMut == 1)] <- 1
dnarGeneAlt[which(dnarGeneExpKo == 1)] <- 2

# Identify genes with most frequent alterations (by mutation or lack of expression).
numAltLines <- apply(dnarGeneAlt, MARGIN = 1, FUN = sum)
sort(numAltLines, decreasing = TRUE)

dnarGeneAltPlotData <- dnarGeneAlt[order(numAltLines), ]
heatmap(dnarGeneAltPlotData, scale="none", Rowv = NA, Colv = NA, 
                col = c("grey", "red", "blue"), 
                main="Altered genes by mutation (red) or low expression (blue)")

For any altered gene, we can identify drugs with significantly different activity between altered and non-altered cell lines.

geneName <- "APLF"
altLineIndices <- which(dnarGeneAlt[geneName, ] != 0)
nonAltLineIndices <- which(dnarGeneAlt[geneName, ] == 0)

drugAssocTab <- data.frame(GENE=geneName, DRUG_NAME=fdaDrugAnnot$NAME, 
                                                     DRUG_NSC=fdaDrugAnnot$NSC, TTEST_PVAL=NA, 
                                                     ADJ_PVAL=NA, MEAN_ACT_DIFF=NA, 
                                                     stringsAsFactors = FALSE)
rownames(drugAssocTab) <- drugAssocTab$DRUG_NSC
for (drugNsc in drugAssocTab$DRUG_NSC){
    ttestResult <- t.test(x=nci60FdaDrugAct[drugNsc, altLineIndices], 
                                                y=nci60FdaDrugAct[drugNsc, nonAltLineIndices])
    drugAssocTab[drugNsc, "TTEST_PVAL"] <- ttestResult$p.value
    drugAssocTab[drugNsc, "MEAN_ACT_DIFF"] <- 
        ttestResult$estimate[1] - ttestResult$estimate[2]
}

drugAssocTab$ADJ_PVAL <- p.adjust(drugAssocTab$TTEST_PVAL, method = "bonferroni")
drugAssocTab <- drugAssocTab[order(drugAssocTab$ADJ_PVAL), ]

drugAssocTab[drugAssocTab$ADJ_PVAL < 0.05, 
                         c("GENE", "DRUG_NAME", "DRUG_NSC", "ADJ_PVAL", "MEAN_ACT_DIFF")]

meanActDiff <- drugAssocTab$MEAN_ACT_DIFF
negLog10Pval <- -log10(drugAssocTab$TTEST_PVAL)
plot(x=meanActDiff, y=negLog10Pval, xlim = c((min(meanActDiff)-1), (max(meanActDiff)+1)),
         xlab="Difference between -logGI50 means (altered vs. non-altered lines)",
         ylab="-log10(p-value)", main=(paste0(geneName, " Drug Response Associations")))

ptLabels <- character(length(drugAssocTab$DRUG_NAME))
numLabeledPts <- 7 # label the k most significant drug associations
ptLabels[1:numLabeledPts] <- drugAssocTab$DRUG_NAME[1:numLabeledPts]
text(x=meanActDiff, negLog10Pval, ptLabels, cex=0.9, pos=4, col="red")

Session Information

sessionInfo()

References

• Reinhold, W.C., et al. (2012) CellMiner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set. Cancer Research, 72, 3499-3511.

• Sousa, F.G., et al. (2015) Alterations of DNA repair genes in the NCi-60 cell lines and their predictive value for anticancer drug activity. DNA Repair 28 (2015) 107-115.

• Varma, S., et al. (2014) High Resolution Copy Number Variation Data in the NCI-60 Cancer Cell Lines from Whole Genome Microarrays Accessible through CellMiner. PLoS ONE 9(3): e92047.

rcellminer documentation built on Nov. 26, 2020, 2:02 a.m.