Introduction

Pharmacogenomics holds much potential to aid in discovering drug response biomarkers and developing novel targeted therapies, leading to development of precision medicine and working towards the goal of personalized therapy. Several large experiments have been conducted, both to molecularly characterize drug dose response across many cell lines, and to examine the molecular response to drug administration. However, the experiments lack a standardization of protocols and annotations, hindering meta-analysis across several experiments.

PharmacoGx was developed to address these challenges, by providing a unified framework for downloading and analyzing large pharmacogenomic datasets which are extensively curated to ensure maximum overlap and consistency. PharmacoGx is based on a level of abstraction from the raw experimental data, and allows bioinformaticians and biologists to work with data at the level of genes, drugs and cell lines. This provides a more intuitive interface and, in combination with unified curation, simplifies analyses between multiple datasets.

To organize the data released by each experiment, we developed the PharmacoSet class. This class efficiently stores different types of data and facilitates interrogating the data by drug or cell line. The PharmacoSet is also versatile in its ability to deal with two distinct types of pharmacogenomic datasets. The first type, known as sensitivity datasets, are datasets where cell lines were profiled on the molecular level, and then tested for drug dose response. The second type of dataset is the perturbation dataset. These types of datasets profile a cell line on the molecular level before and after administration of a compound, to characterize the action of the compound on the molecular level.

With the first release of PharmacoGx we have fully curated and created PharmacoSet objects for three publicly available large pharmacogenomic datasets. Two of these datasets are of the sensitivity type. These are the Genomics of Drug Sensitivity in Cancer Project (GDSC) (Garnett et al. 2012) and the Cancer Cell Line Encyclopedia (CCLE) (Barretina et al. 2012). The third dataset is of the perturbation type, the Connectivity Map (CMAP) project (Lamb et al. 2006).

Furthermore, PharmacoGx provides a suite of parallelized functions which facilitate drug response biomarker discovery, and molecular drug characterization. This vignette will provide two example analysis case studies. The first will be comparing gene expression and drug sensitivity measures across the CCLE and GDSC projects. The second case study will interrogate the CMAP database with a known signature of up and down regulated genes for HDAC inhibitors as published in [Glaser et al. (2003)}. Using the Connectivity Score as defined in (Lamb et al. 2006), it will be seen that known HDAC inhibitors have a high numerical score and high significance.

For the purpose of this vignette, an extremely minuscule subset of all three PharmacoSet objects are included with the package as example data. They are included for illustrative purposes only, and the results obtained with them will likely be meaningless.

Installation and Settings

PharmacoGx requires that several packages are installed. However, all dependencies are available from CRAN or Bioconductor.

if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install('PharmacoGx')

Load PharamacoGx into your current workspace:

library(PharmacoGx)

Requirements

PharmacoGx has been tested on Windows, Mac and Cent OS platforms. The packages uses the core R package parallel to preform parallel computations, and therefore if parallelization is desired, the dependencies for the parallel package must be met.

Downloading PharmacoSet Objects

We have made the PharmacoSet objects of the curated datasets available for download using functions provided in the package. A table of available PharmacoSet objects can be obtained by using the availablePSets() function. Any of the PharmacoSets in the table can then be downloaded by calling downloadPSet, which saves the datasets into a directory of the users choice, and returns the data into the R session.

availablePSets()
GDSC <- downloadPSet("GDSC_2020(v2-8.2)")

Downlading Drug Signatures

The package also provides tools to compute drug perturbation and sensitivity signatures, as explained below. However, the computation of the perturbation signatures is very lengthy, so for users’ convenience we have precomputed the signatures for our perturbation PharmacoSet objects and made them available for download using the function downloadPertSig().

Case Study

(In)Consistency across large pharmacogenomic studies

Our first case study illustrates the functions for analysis of the sensitivity type of dataset. The case study will investigate the consistency between the GDSC and CCLE datasets, recreating the analysis similar to our Inconsistency in Large Pharmacogenomic Studies paper (Haibe-Kains et al. 2013). In both CCLE and GDSC, the transcriptome of cells was profiled using an Affymatrix microarray chip. Cells were also tested for their response to increasing concentrations of various compounds, and form this the IC50 and AUC were computed. However, the cell and drugs names used between the two datasets were not consistent. Furthermore, two different microarray platforms were used. However, PharmacoGx allows us to overcome these differences to do a comparative study between these two datasets.

GDSC was profiled using the hgu133a platform, while CCLE was profiled with the expanded hgu133plus2 platform. While in this case the hgu133a is almost a strict subset of hgu133plus2 platform, the expression information in PharmacoSet objects is summarized by Ensemble Gene Ids, allowing datasets with different platforms to be directly compared. The probe to gene mapping is done using the BrainArray customCDF for each platform (Sabatti, Karsten, and Geschwind 2002).

To begin, you would load the datasets from disk or download them using the downloadPSet function above. In the following example, we use the toy datasets provided with the package to illustrate the process, but to recreate the full analysis the full PharmacoSets have to be downloaded.

We want to investigate the consistency of the data between the two datasets. The common intersection between the datasets can then be found using intersectPSet. We create a summary of the gene expression and drug sensitivity measures for both datasets, so we are left with one gene expression profile and one sensitivity profile per cell line within each dataset. We can then compare the gene expression and sensitivity measures between the datasets using a standard correlation coefficient.

  library(Biobase)
  library(SummarizedExperiment)
  library(S4Vectors)
  library(PharmacoGx)
  data("GDSCsmall")
  data("CCLEsmall")
  commonGenes <- intersect(fNames(GDSCsmall, "rna"),
                           fNames(CCLEsmall,"rna"))
  common <- intersectPSet(list('CCLE'=CCLEsmall,
                               'GDSC'=GDSCsmall),
                          intersectOn=c("cell.lines", "drugs"),
                          strictIntersect=TRUE)


  GDSC.auc <- summarizeSensitivityProfiles(
                common$GDSC,
                sensitivity.measure='auc_published',
                summary.stat="median",
                verbose=FALSE)
  CCLE.auc <- summarizeSensitivityProfiles(
                common$CCLE,
                sensitivity.measure='auc_published',
                summary.stat="median",
                verbose=FALSE)

  GDSC.ic50 <- summarizeSensitivityProfiles(
                common$GDSC,
                sensitivity.measure='ic50_published',
                summary.stat="median",
                verbose=FALSE)
  CCLE.ic50 <- summarizeSensitivityProfiles(
                common$CCLE,
                sensitivity.measure='ic50_published',
                summary.stat="median",
                verbose=FALSE)

  GDSCexpression <- summarizeMolecularProfiles(common$GDSC,
                                        cellNames(common$GDSC),
                                        mDataType="rna",
                                        features=commonGenes,
                                        verbose=FALSE)
  CCLEexpression <- summarizeMolecularProfiles(common$CCLE,
                                         cellNames(common$CCLE),
                                         mDataType="rna",
                                         features=commonGenes,
                                         verbose=FALSE)
  gg <- fNames(common[[1]], 'rna')
  cc <- cellNames(common[[1]])

  ge.cor <- sapply(cc, function (x, d1, d2) {
    stats::cor(d1[ , x], d2[ , x], method="spearman",
                use="pairwise.complete.obs")
  ## TO DO:: Ensure all assays are name so we can call by name instead of index
  }, d1=assay(GDSCexpression, 1), d2=assay(CCLEexpression, 1))
  ic50.cor <- sapply(cc, function (x, d1, d2) {
    stats::cor(d1[, x], d2[ , x], method="spearman",
                use="pairwise.complete.obs")
  }, d1=GDSC.ic50, d2=CCLE.ic50)
  auc.cor <- sapply(cc, function (x, d1, d2) {
    stats::cor(d1[ , x], d2[ , x], method="spearman",
                use="pairwise.complete.obs")
  }, d1=GDSC.auc, d2=CCLE.auc)

  w1 <- stats::wilcox.test(x=ge.cor, y=auc.cor,
                           conf.int=TRUE, exact=FALSE)
  w2 <- stats::wilcox.test(x=ge.cor, y=ic50.cor,
                           conf.int=TRUE, exact=FALSE)
  yylim <- c(-1, 1)
  ss <- sprintf("GE vs. AUC = %.1E\nGE vs. IC50 = %.1E",
                w1$p.value, w2$p.value)

  boxplot(list("GE"=ge.cor,
               "AUC"=auc.cor,
               "IC50"=ic50.cor),
          main="Concordance between cell lines",
          ylab=expression(R[s]),
          sub=ss,
          ylim=yylim,
          col="lightgrey",
          pch=20,
          border="black")

Query the Connectivity Map

The second case study illustrates the analysis of a perturbation type datasets, where the changes in cellular molecular profiles are compared before and after administering a compound to the cell line. Of these datasets, we have currently curated and made available for download the Connectivity Map (CMAP) dataset (Lamb et al. 2006).

For this case study, we will recreate an analysis from the paper by Lamb et al., in which a known signature for HDAC inhibitors (Glaser et al. 2003) is used to recover drugs in the CMAP dataset that are also known HDAC inhibitors. For this example, the package includes this signature, already mapped to the gene level, and it can be loaded by calling data(HDAC\_genes).

Once again, we load the dataset, downloading it if needed using downloadPSet . We then recreate drug signatures for each drug using the function drugPerturbationSig to preform statistical modelling of the transcriptomic response to the application of each drug. We then compare the observed up-regulated and down-regulated genes to a the known HDAC signature, using the GSEA connectivity score to determine the correlation between the two signatures.

  library(PharmacoGx)
  library(pander)
  data(CMAPsmall)
  drug.perturbation <- drugPerturbationSig(CMAPsmall,
                                           mDataType="rna",
                                           verbose=FALSE)
  data(HDAC_genes)

  res <- apply(drug.perturbation[,,c("tstat", "fdr")],
               2, function(x, HDAC){
        return(connectivityScore(x=x,
                                 y=HDAC[,2,drop=FALSE],
                               method="fgsea", nperm=100))
    }, HDAC=HDAC_genes)

  rownames(res) <- c("Connectivity", "P Value")
  res <- t(res)
  res <- res[order(res[,1], decreasing=TRUE),]
  pander::pandoc.table(res,
    caption='Connectivity Score results for HDAC inhibitor gene signature.',
    style = 'rmarkdown')

Connectivity Score results for HDAC inhibitor gene signature.
	Connectivity	P Value
vorinostat	0.9449	0.006998
alvespimycin	0.8831	0.008929
acetylsalicylic acid	0	1
rosiglitazone	-0.714	0.03372
pioglitazone	-0.7335	0.07952

As we can see, the known HDAC inhibitor Varinostat has a very strong connectivity score, as well as a very high significance as determined by permutation testing, in comparison to the other drugs, which score poorly.

This example serves as a validation of the method, and demonstrates the ease with which drug perturbation analysis can be done using PharmacoGx. While in this case we were matching a drug signature with a drug class signature, this method can also be used in the discovery of drugs that are anti-correlated with known disease signatures, to look for potential new treatments and drug repurposing.

Estimating Drug Sensitivity Measures

PharmacoGx includes functions to calculate estimated AUC (Area Under drug response Curve) and IC50 values from drugs dose response experiments that measure cell viability after applications of varying concentrations of drug. Additionally, these measures are recomputed for every sensitivity PharmacoSet we create and included alongside any measures published with the original data. The AUC measures originally published are labelled as auc_published, while the recomputed measures are labelled as auc_recomputed, and likewise for the IC50.

While the PharmacoSets already contain the recomputed data, the AUC and IC50 can be calculated for arbitrary data using the computeIC50 and computeAUC functions. The AUC can be calculated using either the area under the curve defined by the actual points recorded, or the area under the curve fitted to the data.

Curve Fitting

While AUC can be numerically calculated without curve fitting, to estimate the IC50 a drug dose response curve must be fit to the data.The dose-response curves are fitted to the equation $$ y = E_{\infty} + \frac{ 1 - E_{\infty} }{1+(\frac{x}{IC50})^{HS}} $$ where the maximum viability is normalized to be y = y(0) = 1, E_∞ denotes the minimum possible viability achieved by administering any amount of the drug, IC50 is the concentration at which viability is reduced to half of the viability observed in the presence of an arbitrarily large concentration of drug, and HS is a parameter describing the cooperativity of binding. HS < 1 denotes negative binding cooperativity, HS = 1 denotes noncooperative binding, and HS > 1 denotes positive binding cooperativity. The parameters of the curves are fitted using the least squares optimization framework. The fitting of these curves to arbitrary points is implemented by the logLogisticRegression function.

Plotting Drug-Dose Response Data

Drug-Dose response data included in the PharmacoSet objects can be conviniently plotted using the drugDoseResponseCurve function. Given a list of PharmacoSets, a drug name and a cell name, it will plot the drug dose response curves for the given cell-drug combination in each dataset, allowing direct comparisons of data between datasets.

Gene-Drug Association Modelling

PharmacoGx provides methods to model the association between drugs and molecular data such as transcriptomics, genomics and proteomics.

Sensitivity studies allow the discovery of molecular features that improve or inhibit the sensitivity of cell lines to various compounds, by looking at the association between the expression of the feature and the response towards each compound. This allows the selection of drugs to be tailored to each specific patient based on the expressed known sensitivity biomarkers. The function drugSensitivitySig models these associations.

Perturbation studies on the other hand look at the molecular profiles of a cell before and after application of a drug, and therefore can characterize the action of a drug on the molecular level. It is hypothesized that drugs which act to down-regulate expression of known disease biomarkers would be effective in reversing the cell from a diseased to healthy state. The function drugPerturbationSig models the molecular profiles of drugs tested in a perturbation dataset.

Sensitivity Modelling

The association between molecular features and response to a given drug is modelled using a linear regression model adjusted for tissue source:

Y = β₀ + β_iG_i + β_tT + β_bB

where Y denotes the drug sensitivity variable, G_i, T and B denote the expression of gene i, the tissue source and the experimental batch respectively, and βs are the regression coefficients. The strength of gene-drug association is quantified by β_i, above and beyond the relationship between drug sensitivity and tissue source. The variables Y and G are scaled (standard deviation equals to 1) to estimate standardized coefficients from the linear model. Significance of the gene-drug association is estimated by the statistical significance of β_i (two-sided t test). P-values are then corrected for multiple testing using the false discovery rate (FDR) approach.

As an example of the efficacy of the modelling approach, we can model the significance of the association between two drugs and their known biomarkers in CCLE. We examine the association between drug 17-AAG and gene NQO1, as well as drug PD-0325901 and gene BRAF:

  library(pander)
  data(CCLEsmall)
  features <- fNames(CCLEsmall, "rna")[
                          which(featureInfo(CCLEsmall,
                                            "rna")$Symbol == "NQO1")]
  sig.rna <- drugSensitivitySig(object=CCLEsmall,
                            mDataType="rna",
                            drugs=c("17-AAG"),
                            features=features,
                            sensitivity.measure="auc_published",
                            molecular.summary.stat="median",
                            sensitivity.summary.stat="median",
                            verbose=FALSE)
  sig.mut <- drugSensitivitySig(object=CCLEsmall,
                            mDataType="mutation",
                            drugs=c("PD-0325901"),
                            features="BRAF",
                            sensitivity.measure="auc_published",
                            molecular.summary.stat="and",
                            sensitivity.summary.stat="median",
                            verbose=FALSE)
  sig <- rbind(sig.rna, sig.mut)
  rownames(sig) <- c("17-AAG + NQO1","PD-0325901 + BRAF")
  colnames(sig) <- dimnames(sig.mut)[[3]]
  pander::pandoc.table(t(sig), style = "rmarkdown", caption='P Value of Gene-Drug Association' )

P Value of Gene-Drug Association
	17-AAG + NQO1	PD-0325901 + BRAF
estimate	0.601	0.8252
se	0.05368	0.1339
n	492	472
tstat	11.2	6.161
fstat	125.3	37.95
pvalue	6.034e-26	1.613e-09
df	469	449
fdr	6.034e-26	1.613e-09

Perturbation Modelling

The molecular response profile of a given drug is modelled as a linear regression model adjusted experimental batch effects, cell specific differences, and duration of experiment to isolate the effect of the concentration of the drug applied:

G = β₀ + β_iC_i + β_tT + β_dD + β_bB

where G denotes the molecular feature expression (gene), C_i denotes the concentration of the compound applied, T the cell line identity, D denotes the duration of the experiment, B denotes the experimental batch, and βs are the regression coefficients. The strength of feature response is quantified by β_i. Unlike the sensitivity signatures, the G and C variables are unscaled. Significance of the gene-drug association is estimated by the statistical significance of β_i, calculated using an F-test on the improvement in fit after the inclusion of the term. P-values are then corrected for multiple testing using the false discovery rate (FDR) approach.

Connectivity Scoring

The package also provides two methods for quantifying the similarity between two molecular signatures of the form returned by drugPerturbationSig and drugSensitivitySig, or a set of up and down regulated genes as part of a disease signature. The two methods are the GSEA method as introduced by Subramanian et at (Subramanian et al. 2005), and GWC, a method based on a weighted Spearman correlation coefficient. Both methods are implemented by the connectivityScore function.

GSEA

The GSEA method is implemented to compare a signature returned by drugPerturbationSig with a known set of up and down regulated genes in a disease state. For the disease signature, the function expects a vector of features with a value, either binary (1, -1) or continuous, where the sign signifies if the gene is up or down regulated in the disease. The names of the vector are expected to be the feature names, matching the feature names used in the drug signature. The function then returns a GSEA score measuring the concordance of the disease signature to the drug signature, as well as an estimated P-Value for the significance of the connectivity determined by permutation testing.

GWC

The GWC method (Genome Wide Correlation) is implemented to compare two signatures of the same length, such as two drug signatures returned by drugPerturbationSig. The score is a Spearman correlation weighted by the normalized negative logarithm significance of each value. The normalization is done so that datasets of different size can be compared without concern for lower p-values in one dataset due to larger sample sizes.

More precisely, take X_i and Y_i to be the ranks of the first and second set of data respectively, and Px_i, Py_i to be the p-values of those observations. The weight for each pair of observations is:

$$ Wx_i = \frac{-\log_{10}(Px_i)}{\sum_{i}-\log_{10}(Px_i)} $$

$$ Wy_i = \frac{-\log_{10}(Py_i)}{\sum_{i}-\log_{10}(Py_i)} $$

W_i = Wx_i + Wy_i

If we further define the weighted mean as follows: $$ m(X;W) = \frac{\sum_i W_i X_i}{\sum_i{W_i}} $$

Then the weighted correlation is given by:

$$ cor(X,Y,W) = \frac{\frac{\sum_i W_i (X_i - m(X;W))(Y_i - m(Y,W))}{\sum_i W_i}}{\sqrt{(\frac{\sum_i W_i (X_i - m(X;W))^2}{\sum_i W_i}) (\frac{\sum_i W_i (Y_i - m(Y;W))^2}{\sum_i W_i})}} $$

This correlation therefore takes into account not only the ranking of each feature in a signature, but also of the significance of each rank.

Acknowledgements

The authors of the package would like to thank the investigators of the Genomics of Drug Sensitivity in Cancer Project, the Cancer Cell Line Encyclopedia and the Connectivity Map Project who have made their invaluable data available to the scientific community. We would also like to thank Mark Freeman for contributing the code for MLE fitting drug-dose response curves towards this package. We are indebted to Nehme El-Hachem, Donald Wang and Adrian She for their contributions towards the accurate curation of the datasets. Finally, it is important to acknowledge all the members of the Benjamin Haibe-Kains lab for their contribution towards testing and feedback during the design of the package.

Session Info

# set eval = FALSE if you don't want this info (useful for reproducibility) to appear
sessionInfo()

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] pander_0.6.5                ggplot2_3.5.1              
 [3] data.table_1.16.4           PharmacoGx_3.11.0          
 [5] CoreGx_2.11.0               SummarizedExperiment_1.37.0
 [7] Biobase_2.67.0              GenomicRanges_1.59.1       
 [9] GenomeInfoDb_1.43.2         IRanges_2.41.2             
[11] S4Vectors_0.45.2            MatrixGenerics_1.19.0      
[13] matrixStats_1.4.1           BiocGenerics_0.53.3        
[15] generics_0.1.3              knitcitations_1.0.12       
[17] knitr_1.49                  BiocStyle_2.35.0           

loaded via a namespace (and not attached):
  [1] bitops_1.0-9                rlang_1.1.4                
  [3] magrittr_2.0.3              shinydashboard_0.7.2       
  [5] compiler_4.4.2              reshape2_1.4.4             
  [7] vctrs_0.6.5                 relations_0.6-14           
  [9] stringr_1.5.1               pkgconfig_2.0.3            
 [11] crayon_1.5.3                fastmap_1.2.0              
 [13] backports_1.5.0             XVector_0.47.0             
 [15] caTools_1.18.3              utf8_1.2.4                 
 [17] promises_1.3.2              rmarkdown_2.29             
 [19] UCSC.utils_1.3.0            coop_0.6-3                 
 [21] xfun_0.49                   MultiAssayExperiment_1.33.1
 [23] zlibbioc_1.52.0             cachem_1.1.0               
 [25] jsonlite_1.8.9              RefManageR_1.4.0           
 [27] SnowballC_0.7.1             later_1.4.1                
 [29] DelayedArray_0.33.3         BiocParallel_1.41.0        
 [31] parallel_4.4.2              sets_1.0-25                
 [33] cluster_2.1.8               R6_2.5.1                   
 [35] RColorBrewer_1.1-3          bslib_0.8.0                
 [37] stringi_1.8.4               limma_3.63.2               
 [39] boot_1.3-31                 lubridate_1.9.4            
 [41] jquerylib_0.1.4             Rcpp_1.0.13-1              
 [43] R.utils_2.12.3              downloader_0.4             
 [45] httpuv_1.6.15               Matrix_1.7-1               
 [47] igraph_2.1.2                timechange_0.3.0           
 [49] tidyselect_1.2.1            abind_1.4-8                
 [51] yaml_2.3.10                 gplots_3.2.0               
 [53] codetools_0.2-20            lattice_0.22-6             
 [55] tibble_3.2.1                plyr_1.8.9                 
 [57] withr_3.0.2                 shiny_1.10.0               
 [59] BumpyMatrix_1.15.0          evaluate_1.0.1             
 [61] bench_1.1.3                 xml2_1.3.6                 
 [63] pillar_1.9.0                lsa_0.73.3                 
 [65] BiocManager_1.30.25         KernSmooth_2.23-24         
 [67] DT_0.33                     checkmate_2.3.2            
 [69] shinyjs_2.1.0               piano_2.23.0               
 [71] munsell_0.5.1               scales_1.3.0               
 [73] gtools_3.9.5                xtable_1.8-4               
 [75] marray_1.85.0               glue_1.8.0                 
 [77] slam_0.1-55                 maketools_1.3.1            
 [79] tools_4.4.2                 sys_3.4.3                  
 [81] fgsea_1.33.0                visNetwork_2.1.2           
 [83] buildtools_1.0.0            fastmatch_1.1-4            
 [85] cowplot_1.1.3               grid_4.4.2                 
 [87] bibtex_0.5.1                colorspace_2.1-1           
 [89] GenomeInfoDbData_1.2.13     cli_3.6.3                  
 [91] fansi_1.0.6                 S4Arrays_1.7.1             
 [93] dplyr_1.1.4                 gtable_0.3.6               
 [95] R.methodsS3_1.8.2           sass_0.4.9                 
 [97] digest_0.6.37               SparseArray_1.7.2          
 [99] htmlwidgets_1.6.4           R.oo_1.27.0                
[101] htmltools_0.5.8.1           lifecycle_1.0.4            
[103] httr_1.4.7                  statmod_1.5.0              
[105] mime_0.12

References

Barretina, Jordi, Giordano Caponigro, Nicolas Stransky, Kavitha Venkatesan, Adam A Margolin, Sungjoon Kim, Christopher J Wilson, et al. 2012. “The Cancer Cell Line Encyclopedia Enables Predictive Modelling of Anticancer Drug Sensitivity.” Nature 483 (7391): 603–7.

Garnett, Mathew J, Elena J Edelman, Sonja J Heidorn, Chris D Greenman, Anahita Dastur, King Wai Lau, Patricia Greninger, et al. 2012. “Systematic Identification of Genomic Markers of Drug Sensitivity in Cancer Cells.” Nature 483 (7391): 570–75.

Glaser, Keith B, Michael J Staver, Jeffrey F Waring, Joshua Stender, Roger G Ulrich, and Steven K Davidsen. 2003. “Gene Expression Profiling of Multiple Histone Deacetylase (HDAC) Inhibitors: Defining a Common Gene Set Produced by HDAC Inhibition in T24 and MDA Carcinoma Cell Lines.” Molecular Cancer Therapeutics 2 (2): 151–63.

Haibe-Kains, Benjamin, Nehme El-Hachem, Nicolai Juul Birkbak, Andrew C Jin, Andrew H Beck, Hugo J W L Aerts, and John Quackenbush. 2013. “Inconsistency in Large Pharmacogenomic Studies.” Nature 504 (7480): 389–93.

Lamb, Justin, Emily D. Crawford, David Peck, Joshua W. Modell, Irene C. Blat, Matthew J. Wrobel, Jim Lerner, et al. 2006. “The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease.” Science 313 (5795): 1929–35. https://doi.org/10.1126/science.1132939.

Sabatti, Chiara, Stanislav L Karsten, and Daniel H Geschwind. 2002. “Thresholding Rules for Recovering a Sparse Signal from Microarray Experiments.” Mathematical Biosciences 176 (1): 17–34.

Subramanian, Aravind, Pablo Tamayo, Vamsi K. Mootha, Sayan Mukherjee, Benjamin L. Ebert, Michael A. Gillette, Amanda Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proceedings of the National Academy of Sciences 102 (43): 15545–50. https://doi.org/10.1073/pnas.0506580102.

PharmacoGx: An R Package for Analysis of Large Pharmacogenomic Datasets