TOPS: a versatile software tool for statistical analysis and visualization of combinatorial gene-gene and gene-drug interaction screens
- Markus K Muellner†1Email author,
- Gerhard Duernberger†1, 2,
- Florian Ganglberger†1,
- Claudia Kerzendorfer1,
- Iris Z Uras1,
- Andreas Schoenegger1,
- Klaudia Bagienski1,
- Jacques Colinge1 and
- Sebastian MB Nijman1Email author
© Muellner et al.; licensee BioMed Central Ltd. 2014
Received: 17 May 2013
Accepted: 28 March 2014
Published: 8 April 2014
Measuring the impact of combinations of genetic or chemical perturbations on cellular fitness, sometimes referred to as synthetic lethal screening, is a powerful method for obtaining novel insights into gene function and drug action. Especially when performed at large scales, gene-gene or gene-drug interaction screens can reveal complex genetic interactions or drug mechanism of action or even identify novel therapeutics for the treatment of diseases.
The result of such large-scale screen results can be represented as a matrix with a numeric score indicating the cellular fitness (e.g. viability or doubling time) for each double perturbation. In a typical screen, the majority of combinations do not impact the cellular fitness. Thus, it is critical to first discern true "hits" from noise. Subsequent data exploration and visualization methods can assist to extract meaningful biological information from the data. However, despite the increasing interest in combination perturbation screens, no user friendly open-source program exists that combines statistical analysis, data exploration tools and visualization.
We developed TOPS (Tool for Combination Perturbation Screen Analysis), a Java and R-based software tool with a simple graphical user interface that allows the user to import, analyze, filter and plot data from double perturbation screens as well as other compatible data. TOPS was designed in a modular fashion to allow the user to add alternative importers for data formats or custom analysis scripts not covered by the original release.
We demonstrate the utility of TOPS on two datasets derived from functional genetic screens using different methods. Dataset 1 is a gene-drug interaction screen and is based on Luminex xMAP technology. Dataset 2 is a gene-gene short hairpin (sh)RNAi screen exploring the interactions between deubiquitinating enzymes and a number of prominent oncogenes using massive parallel sequencing (MPS).
TOPS provides the benchtop scientist with a free toolset to analyze, filter and visualize data from functional genomic gene-gene and gene-drug interaction screens with a flexible interface to accommodate different technologies and analysis algorithms in addition to those already provided here. TOPS is freely available for academic and non-academic users and is released as open source.
KeywordsDouble perturbation screens Synthetic lethality Functional genetics Massive parallel sequencing Luminex xMAP Drug screens Synergy Epistasis
Genes operate in complex cellular networks, and elucidating this connectivity is critical for understanding normal physiology and disease. Functional genomic screens that combine two perturbations have been used with great success to uncover such genetic interactions, and also to reveal mechanisms of drug action [1, 2]. This is generally done by perturbing a biological system, most often a cell, in a defined manner i.e. by overexpression or knockdown of a gene or by addition of a drug and subsequently administering a second set of perturbations and analyzing the system for non-additive readout (i.e. synergistic/antagonistic or synthetic lethality/synthetic rescue) .
These perturbation experiments have traditionally been performed in model organisms such as yeast where genetic manipulation is relatively easy compared to mammalian cells . Now, with the advent of technologies including RNA interference (RNAi) and the recent emergence of nuclease-based genome engineering, human cells can also be used to perform functional genomic screens. Cellular fitness, often represented by cell number, division rate or death, is used as a proxy for many of these screens since it is a parameter that can be measured with relative ease. This proxy is also of particular relevance in case genetic interactions are sought as potential anti-cancer therapies. Here, the aim is to kill only cells carrying cancer mutations by targeting an interacting gene which is in synthetic lethal configuration to the mutation. This can be achieved by targeting this second gene product with a small molecule inhibitor or antibody. The non-cancer cells would remain unaffected from this therapy as inhibiting only one interaction partner will not lead to lethality. Indeed, many cancer mutations can be tested systematically against drugs to uncover selective resistance/sensitivity mechanisms, or against small interfering RNAs to find potential new drug targets [4–6].
In practice, the number of data points typically generated is further increased to T x n by the requirement of n replicate measurements. Thus, a screen of, for example, 100 genes against 100 drugs measured in quadruplicate requires 40,000 measurements. A Dataset of this size is usually beyond the limits of what can be conveniently analyzed without the use of dedicated statistics software.
To generate these large datasets, multiplexing strategies greatly improve the throughput of the experiment, by combining genetic perturbations in a single well (multiplexing drug treatments is typically not feasible). A convenient multiplexing strategy is to introduce DNA barcodes into the genome of mutant cell lines such that each genetic perturbation corresponds to a single barcode sequence. In this scenario DNA barcodes act both as unique identifiers for each genetic perturbation and as a proxy for cell number/fitness. Barcodes can be encoded by short hairpin (sh)RNA sequences themselves or introduced into the genome as non-transcribed DNA sequences using lentiviral vectors [7, 8]. Practically, a set of barcoded cell lines (A) can then be pooled together and divided into aliquots that are then screened against the second set (B) of perturbations.
Luminex xMAP – a hybridization based multiplexing technology that allows up to 500 measurements in a single sample. Here, DNA barcoded panels consisting of genetically modified cell lines (i.e. expressing an oncogene or knockdown of a tumor suppressor gene) and an unmodified control from the same genetic background are pooled in a single well of a multiwell plate (genetic perturbation, set A) and treated with either a drug or a second genetic perturbation (drug perturbation, set B). In this case only perturbation set A is multiplexed.
Massive Parallel Sequencing (MPS). In this case, a panel of non-barcoded cell lines carrying perturbation set A from (1) are infected with a library of lentiviral shRNA vectors (genetic perturbation, set B). Here, the shRNA sequence that integrate into the genome serves both as perturbation and as DNA barcode. Further multiplexing can be achieved by indexing the second perturbation in a single sequencing run. For example, a MPS run yielding 108 aligned reads on average allows 100 cell lines to be screened against 100 drugs at a reasonable average sequencing depth of 10,000 reads per DNA barcode.
Programming Language, external software and system requirements
TOPS was developed as a graphical user interface for an analysis pipeline that we originally developed in R. In order to make these R scripts accessible to non-expert users, in particular scientists with little background in scripting languages, we took advantage of the R package rJava (http://www.rforge.net/rJava/) to run the scripts from a Java based interface. TOPS can therefore run its R scripts on all platforms for which rJava is available (currently MacOS, Linux and Windows).
Output from the R interface, including and error messages, is stored in text files for troubleshooting. The overall progress of a session can be monitored and is logged separately.
User skills and installation
Installation of TOPS consists of unpacking the zipped file into a directory. For TOPS to run, Java, R and the package rJava are must also be installed. Upon first launch, definition of the paths in which R and rJava are located may be required. The most common locations for MacOS, Linux and Windows are checked automatically and the user will be prompted only if the libraries cannot be found. For convenience, alternative locations to be checked in this manner can be added in libPaths.txt in the TOPS directory.
TOPS uses a single semicolon separated values (CSV) data input file. This file has to be generated by the user and can for instance be exported from a Microsoft Office Excel spreadsheet. Alternatively the example datasets can be used for testing. To use TOPS no scripting skills are necessary.
Results and discussion
Example datasets provided
TOPS comes with two datasets to test the different procedures. Thus, users can explore the software without having to generate original data.
Example dataset 1: a gene-drug interaction screen in breast cancer
To simulate oncogenic events in breast cancer we created a panel of 70 different isogenic cell lines ectopically expressing oncogenes or knocking down tumor suppressor genes commonly found in breast cancer. Findings based on this dataset have been published previously .
For multiplexing purposes, each cell line was infected with a vector introducing a genetic barcode (xMAP tag) into the genomic DNA. The cell lines were then pooled, distributed among multi-well plates and treated with 87 individual drugs. After a period of exponential growth, genomic DNA was harvested and barcodes were quantified using Luminex xMAP beads. Raw bead data was converted to median measurements per interaction by calculating the truncated mean. The resulting data is in the format ‘Cell line (name)’, ‘Drug (name)’, ‘Barcode score (numerical)’. The barcode score is a measure of the relative fitness as it is an estimate of the total number of cells for each isogenic cell line under the drug treatment condition.
Example dataset 2: a gene-gene interactions screen in breast cancer
We screened 27 isogenic cell lines expressing oncogenes against a library consisting of approximately 400 shRNA vectors covering 80 deubiquitinating enzymes (DUBs). DUBs represent an emerging class of cancer targets involved in many different cellular processes with therapeutic potential. After a period of exponential growth cells were harvested and genomic DNA was extracted in order to PCR-amplify the shRNA barcode cassettes. During the amplification process primers with 60 unique barcode sequences corresponding to the 27 cell lines in 4 replicates were used as described in . The tagged PCR samples were then pooled, a sequencing library was generated and the samples were sequenced on two lanes (50 basepairs) of a Hiseq 2000 instrument (Illumina). The MPS raw data were converted to FASTQ format and split according to the indexing PCR primers with a custom Python script. Each sample’s NGS reads were aligned to the shRNA sequences using shALIGN  allowing for zero, one, and two mismatches. The output CSV file is in the format: ‘Cell line (name)’, ‘shRNA (name)’, ‘number of reads’. Number of reads corresponds to shRNA abundance and this is the readout for relative fitness.
A typical input file consists of five columns with the following information:
perturbation A; perturbation B; fitness score; replicate 1; replicate 2
In the provided example dataset 1, perturbation A corresponds to the gene (i.e. overexpression or knockdown) and perturbation B is a drug treatment. For dataset 2, perturbation B is an shRNA vector targeting a specific DUB and perturbation A is the gene. Fitness score represents the barcode abundance in these examples. Replicates 1 and 2 designate if measurements are replicates of each other. Replicate 1 identifies multiple perturbations towards the same outcome (i.e. multiple different hairpins for one gene, multiple vectors encoding the same cDNA) while replicate 2 identifies the replicate number for a given A-B combination.
TOPS has custom importers covering both types of files from the example datasets. These importers convert the data coming from Luminex or MPS to the input file format.
Entries with non-unique identifiers for perturbation A and perturbation B are automatically interpreted to be replicate measurements by the importers.
TOPS is supplied with two analysis models that are well suited for interaction data of the type that is supplied in the example datasets. However, custom analysis methods can be conveniently plugged into TOPS. The first method is based on a multiple linear regression model and the second one is based on the Mann-Whitney U Test.
Statistical analysis: linear model
Linear models have been used previously to normalize [12, 13] and analyze [14, 15] high-throughput measurements. Accordingly, our analysis is composed of two consecutive linear models: the first one is used for data normalization; the second one identifies outliers that represent "hits" i.e. potentially biologically meaningful interactions that warrant further study.
the systematic effect of perturbation A on fitness
the systematic effect of perturbation B on fitness
the signal that is due to the (synergistic) combination of perturbation A and B consistent across biological replicates, which is the signal of biological interest
noise caused by biological and technical variability between the replicates.
where intensity i represents the measured fitness score, β A , β B are regression coefficients, A, B are categorical variables representing perturbations A and B respectively, and s i is a remaining signal including both noise (biological and technical) and a potential synergistic contribution of A and B perturbations. The model is fit to the experimental data to estimate the first two components from the fitness score corresponding to the systematic, non-synergistic perturbations induced by A and B. The linear model was fit by robust regression using an M estimator (function rlm in R library "MASS") to avoid influence of outliers.
The coefficients of these linear models (βA|B and βB|A) capture the effect of a perturbation within the context of a specific other perturbation. The coefficients combine the signal of replicate experiments and estimate the magnitude of the biological effect thanks to the noise component ε.
The ‘rlm’ function of the R ‘MASS’ library is employed to compute a robust estimate of the coefficients of the linear model. This function also offers the possibility to calculate the statistical significance of each coefficient by performing a one-sided t-test (βA|B ≠ 0 or βB|A ≠ 0). For small sample sizes this test can suffer from false positives due to underestimating variance as described in Axelsson et al. .
As both linear models are applied after log transformation of the input data these additive coefficients reflect modeling a multiplicative effect of the fitness response of the pairwise fitness effects β A and β B [17, 18].
In the case of the sLM split linear model aimed at reducing computation time and memory usage, the dataset is subdivided into groups of 50 perturbations (A and B), resulting in less data points available for the t-test, hence the reduced statistical power.
Statistical analysis: Mann-Whitney U-test
The second method we implemented to analyze fitness score data is based on the Mann-Whitney U-test, a more robust alternative to the Student’s t-test that does not assume a normal distribution. For this analysis a distinct normalization procedure is employed.
Here, i and j represent row and column positions in the matrix, and A j and B i represent sets containing all values for row j and column i, respectively. Thus, each value in a row is divided by the corresponding row median, and each value in a column by the column median across the matrix. Each fitness score in the matrix is subsequently log10 transformed.
Next, the algorithm performs two-sided tests to calculate p-values. The model requires at least 3 replicates per interaction to perform this test. As with the linear method algorithm we calculate two p-values for each interaction, corresponding to perturbation set A or B.
Importing analyzed data into the "results browser" and generating plots
Once data has been analyzed by either method, the results of the analysis are stored in a text file in the following format:
perturbation A; perturbation B; relative magnitude A; relative magnitude B; absolute magnitude A; absolute magnitude B; p-value A; p-value B; combined p-value by Brown’s Method.
Where perturbation A and B are the identifiers for the tested perturbations, relative magnitude A and B gives an indication of the ratio of treatment versus the median of all treatments (control), absolute magnitude A and B are the raw signals from the assay and the p-values are calculated as discussed previously. Given this format, other datasets not generated with our analysis pipeline can also be imported, filtered and plotted using TOPS.
Comparison of performance of the analysis models
In order to compare the performance of all three models (linear, split-linear, and Mann-Whitney U-test) we generated a randomized dataset with a distribution based on dataset 1. In this set we generated true positive hits with an average signal of 0.5, 0.7 (reduced fitness) 1.3 and 1.5 fold (increased fitness) compared to the average of true negatives.
Experimentally validated true and false positive hits from dataset 1
We present a user-friendly software package to analyze and visualize interaction data from functional genomic dual perturbation screens. Although multiple tools for analyzing cell viability screens are available [16, 21–24] these have their limitations by either being based on commercial software (mostly MATLAB) or requiring command line skills. TOPS incorporates statistical models designed for the analysis of pairwise interactions of larger gene/drug sets. Furthermore it is fully based on free software and provides a graphical user interface. The software is easily accessible and offers a powerful analysis tool for the benchtop scientist while being expandable enough to be attractive to users who would like to run their own analysis methods. Importantly, not only Luminex xMAP and Sequencing data can be analyzed with the presented methods but in principle any data from other technologies can be imported as long as the data can be reduced to a "perturbation A", "perturbation B", "score" format and true positives are relatively rare in comparison to true negatives which is a necessity due to the nature of normalization of the in-built analysis. We have included two analysis pipelines based on different methods to demonstrate the versatility of TOPS. We have also included two importers for Luminex xMAP data and for pre-processed screening data that makes these two technologies particularly easy to use with the software.
Availability and requirements
Project name: TOPS.
Project home page: https://sourceforge.net/p/topscemm/wiki/Home/.
Operating system(s): Win32/OSX/Linux.
Programming language: Java and R.
Other requirements: Java 6 or newer. R 2.14 or newer.
License: Creative Commons Attribution ShareAlike License V3.0.
Any restrictions to use by non-academics: Only those imposed by the license.
We would like to thank Luminex Corp. for support with the Luminex assay, Helen Pickersgill for help in the preparation of this manuscript, and colleagues from CeMM who acted as beta-testers for the software.
- Hillenmeyer ME, Ericson E, Davis RW, Nislow C, Koller D, Giaever G: Systematic analysis of genome-wide fitness data in yeast reveals novel gene function and drug action. Genome Biol. 2010, 11 (3): R30-10.1186/gb-2010-11-3-r30.View ArticlePubMed CentralPubMed
- Hillenmeyer ME, Fung E, Wildenhain J, Pierce SE, Hoon S, Lee W, Proctor M, St Onge RP, Tyers M, Koller D, Altman RB, Davis RW, Nislow C, Giaever G: The chemical genomic portrait of yeast: uncovering a phenotype for all genes. Science. 2008, 320 (5874): 362-365. 10.1126/science.1150021.View ArticlePubMed CentralPubMed
- Nijman SM: Synthetic lethality: general principles, utility and detection using genetic screens in human cells. FEBS Lett. 2010, 585 (1): 1-6.View ArticlePubMed
- Kessler JD, Kahle KT, Sun T, Meerbrey KL, Schlabach MR, Schmitt EM, Skinner SO, Xu Q, Li MZ, Hartman ZC, Rao M, Yu P, Dominguez-Vidana R, Liang AC, Solimini NL, Bernardi RJ, Yu B, Hsu T, Golding I, Luo J, Osborne CK, Creighton CJ, Hilsenbeck SG, Schiff R, Shaw CA, Elledge SJ, Westbrook TF: A SUMOylation-dependent transcriptional subprogram is required for Myc-driven tumorigenesis. Science. 2011, 335 (6066): 348-353.View ArticlePubMed CentralPubMed
- Scholl C, Fröhling S, Dunn IF, Schinzel AC, Barbie DA, Kim SY, Silver SJ, Tamayo P, Wadlow RC, Ramaswamy S, Döhner K, Bullinger L, Sandy P, Boehm JS, Root DE, Jacks T, Hahn WC, Gilliland DG: Synthetic lethal interaction between oncogenic KRAS dependency and STK33 suppression in human cancer cells. Cell. 2009, 137 (5): 821-834. 10.1016/j.cell.2009.03.017.View ArticlePubMed
- Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, Fröhling S, Chan EM, Sos ML, Michel K, Mermel C, Silver SJ, Weir BA, Reiling JH, Sheng Q, Gupta PB, Wadlow RC, Le H, Hoersch S, Wittner BS, Ramaswamy S, Livingston DM, Sabatini DM, Meyerson M, Thomas RK, Lander ES, et al: Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009, 462 (7269): 108-112. 10.1038/nature08460.View ArticlePubMed CentralPubMed
- Muellner MK, Uras IZ, Gapp BV, Kerzendorfer C, Smida M, Lechtermann H, Craig-Mueller N, Colinge J, Duernberger G, Nijman SM: A chemical-genetic screen reveals a mechanism of resistance to PI3K inhibitors in cancer. Nat Chem Biol. 2011, 7 (11): 787-793.View ArticlePubMed CentralPubMed
- Brummelkamp TR, Bernards R: New tools for functional mammalian cancer genetics. Nat Rev Cancer. 2003, 3 (10): 781-789. 10.1038/nrc1191.View ArticlePubMed
- Marcotte R, Brown KR, Suarez F, Sayad A, Karamboulas K, Krzyzanowski PM, Sircoulomb F, Medrano M, Fedyshyn Y, Koh JL, van Dyk D, Fedyshyn B, Luhova M, Brito GC, Vizeacoumar FJ, Vizeacoumar FS, Datti A, Kasimer D, Buzina A, Mero P, Misquitta C, Normand J, Haider M, Ketela T, Wrana JL, Rottapel R, Neel BG, Moffat J: Essential gene profiles in breast, pancreatic, and ovarian cancer cells. Cancer Discov. 2012, 2 (2): 172-189. 10.1158/2159-8290.CD-11-0224.View ArticlePubMed
- Sims D, Mendes-Pereira AM, Frankum J, Burgess D, Cerone MA, Lombardelli C, Mitsopoulos C, Hakas J, Murugaesu N, Isacke CM, Fenwick K, Assiotis I, Kozarewa I, Zvelebil M, Ashworth A, Lord CJ: High-throughput RNA interference screening using pooled shRNA libraries and next generation sequencing. Genome Biol. 2011, 12 (10): R104-10.1186/gb-2011-12-10-r104.View ArticlePubMed CentralPubMed
- Tischler J, Lehner B, Fraser AG: Evolutionary plasticity of genetic interaction networks. Nat Genet. 2008, 40 (4): 390-391. 10.1038/ng.114.View ArticlePubMed
- Fang Y, Brass A, Hoyle DC, Hayes A, Bashein A, Oliver SG, Waddington D, Rattray M: A model-based analysis of microarray experimental error and normalisation. Nucleic Acids Res. 2003, 31 (16): e96-10.1093/nar/gng097.View ArticlePubMed CentralPubMed
- Schadt EE, Li C, Ellis B, Wong WH: Feature extraction and normalization algorithms for high-density oligonucleotide gene expression array data. J Cell Biochem Suppl. 2001, Suppl 37: 120-125.View ArticlePubMed
- Huber P: Robust Statistics. 2003, Hoboken, NJ, USA: Wiley-Interscience
- Marazzi A: Algorithms, Routines, and S-Functions for Robust Statistics. 1993, London, UK: Taylor & Francis
- Axelsson E, Sandmann T, Horn T, Boutros M, Huber W, Fischer B: Extracting quantitative genetic interaction phenotypes from matrix combinatorial RNAi. BMC Bioinforma. 2011, 12: 342-10.1186/1471-2105-12-342.View Article
- Mani R, St Onge RP, Hartman JL, Giaever G, Roth FP: Defining genetic interaction. Proc Natl Acad Sci U S A. 2008, 105 (9): 3461-3466. 10.1073/pnas.0712255105.View ArticlePubMed CentralPubMed
- Laufer C, Fischer B, Billmann M, Huber W, Boutros M: Mapping genetic interactions in human cancer cells with RNAi and multiparametric phenotyping. Nat Methods. 2013, 10 (5): 427-431. 10.1038/nmeth.2436.View ArticlePubMed
- Brown M: A method for combining non-independent, one-sided tests of significance. Biometrics. 1975, 31 (4): 987-992. 10.2307/2529826.View Article
- Ilouga PE, Hesterkamp T: On the prediction of statistical parameters in high-throughput screening using resampling techniques. J Biomol Screen. 2012, 17 (6): 705-712. 10.1177/1087057112441623.View ArticlePubMed
- Boutros M, Bras LP, Huber W: Analysis of cell-based RNAi screens. Genome Biol. 2006, 7 (7): R66-10.1186/gb-2006-7-7-r66.View ArticlePubMed CentralPubMed
- Luo B, Cheung HW, Subramanian A, Sharifnia T, Okamoto M, Yang X, Hinkle G, Boehm JS, Beroukhim R, Weir BA, Mermel C, Barbie DA, Awad T, Zhou X, Nguyen T, Piqani B, Li C, Golub TR, Meyerson M, Hacohen N, Hahn WC, Lander ES, Sabatini DM, Root DE: Highly parallel identification of essential genes in cancer cells. Proc Natl Acad Sci U S A. 2008, 105 (51): 20380-20385. 10.1073/pnas.0810485105.View ArticlePubMed CentralPubMed
- Collins SR, Schuldiner M, Krogan NJ, Weissman JS: A strategy for extracting and analyzing large-scale quantitative epistatic interaction data. Genome Biol. 2006, 7 (7): R63-10.1186/gb-2006-7-7-r63.View ArticlePubMed CentralPubMed
- Baryshnikova A, Costanzo M, Kim Y, Ding H, Koh J, Toufighi K, Youn JY, Ou J, San Luis BJ, Bandyopadhyay S, Hibbs M, Hess D, Gingras AC, Bader GD, Troyanskaya OG, Brown GW, Andrews B, Boone C, Myers CL: Quantitative analysis of fitness and genetic interactions in yeast on a genome scale. Nat Methods. 2010, 7 (12): 1017-1024. 10.1038/nmeth.1534.View ArticlePubMed CentralPubMed
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.