Overview of the approach. We present a computational framework that maps the functional genomic landscape in MM for a given drug by relying on MM-patient RNA sequencing (RNA-seq), whole exome sequencing (WES), and cytogenetics data that are matched with ex vivo drug sensitivity measures of CD138+ cells isolated from patients. An overview of the flow of information across various stages of this framework is depicted in Fig. 5a, which begins with a TCC-consented (Total Cancer Care) patient donating their bone marrow (BM) specimen. The sorted CD138+ cells are characterized by performing WES, RNA-seq, and ex vivo drug sensitivity screening of several MM drugs. Further, the same patient’s clinical data is abstracted to transcribe their treatment history and fluorescence in-situ hybridization (FISH) data into a well annotated database that can be used to programmatically query this information. The patient’s cytogenetic abnormalities and mutated genes identified from FISH and WES, respectively, characterize an MM patient’s tumor genetically (presence or absence of a genetic event). RNA-seq data alone characterizes the transcriptome of the patient’s tumor, leading to an MM-specific transcriptomic landscape. Finally, ex vivo drug sensitivity metrics like LD50 (lethal dose for median effect) or LD20, and AUC (area under the dose response curve) characterize an MM patient’s ex vivo drug sensitivity to each drug tested. We identify all functional (associated with a phenotype) genetic events that lead to sensitivity, or resistance to a given drug across the entire MM cohort, resulting in genetic biomarkers that can putatively predict ex vivo drug sensitivity in MM. We then rely on paired RNA-seq and ex vivo drug sensitivity data from patient samples to infer transcriptomic footprints, or biomarkers for resistance, or sensitivity to each drug. These inferred biomarkers are examined to develop novel therapeutic strategies in MM.
Primary cancer cells: Moffitt cohort & MCC17814. We investigated the ex vivo responses of cancer cells from 399 MM patients. Out of those, 265 had molecular data (WES). Demographic information can be found in Table 1. Investigators obtained signed informed consent from all patients who were enrolled in the clinical trials/protocols MCC17814, MCC14745, MCC14690, and MCC18608 conducted at the H. Lee Moffitt Cancer Center and Research Institute, as approved by the institutional review board. To this end, patient samples were used in accordance with the Declaration of Helsinki, International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS), Belmont Report, and U.S. Common Rule. The medical records were deidentified and only the following clinically relevant information was reviewed: (A) the treatment administered (therapeutic agents, doses, and schedule) prior to biopsy, (B) cytogenetics, and (C) serum and urine electrophoresis results.
MM patient data repository: Moffitt cohort. At Moffitt Cancer Center, we collected over 1,000 MM patient samples in total and the tumor cells were characterized for ex vivo drug sensitivity, molecular features, and clinical traits in partnership with the ORIEN/AVATAR consortium. Fig. 5b presents a Venn diagram featuring 399 samples screened for ex vivo drug sensitivity, 918 samples with WES, 930 samples with RNA-seq, and 724 samples with cytogenetics data abstracted. Out of all these samples 716 of them have WES data, RNA-seq data, and cytogenetics; 199 of them also have ex vivo drug sensitivity data. The Venn diagram depicts not only the depth of the database (number of samples) but also the breadth (types of data) within it.
Extended Data Fig. 1 presents an oncoplot for patients that have ex vivo drug sensitivity, RNA-seq, and WES data (239 patients), the top ten most common mutations by frequency as rows and patients as columns, which are ordered by their disease states from smoldering multiple myeloma, NDMM, ERMM, to LRMM. The patients’ cytogenetic abnormalities are also highlighted.
Ex vivo drug sensitivity characterization: Moffitt cohort. An ex vivo assay was used to quantify the chemosensitivity of primary MM cells. Fresh BM aspirate cells were enriched for CD138+ expression using Miltenyi (Bergisch Gladbach, Germany) 130-051-301 antibody-conjugated magnetic beads. MM cells (CD138+) were seeded in Corning (Corning, NY) CellBIND 384 well plates with collagen I and previously established human-derived stroma, containing approximately 4000 MM cells and 1000 stromal cells. Each well was filled with 80 μL of Roswell Park Memorial Institute (RPMI) 1640 media supplemented with fetal bovine serum (FBS, heat inactivated), penicillin/streptomycin, and patient-derived plasma (10%, freshly obtained from patient's own aspirate, filtered) and left overnight for adhesion of stroma. The next day, drugs were added using a robotic plate handler so that every drug/combination was tested at 5 (fixed concentration ratio, for combinations) concentrations (1:3 serial dilution) in two replicates. Negative controls (supplemented growth media with and without the vehicle control dimethyl sulfoxide [DMSO]) were included, as well as positive controls for each drug (cell line MM1.S at highest drug concentration). Plates were placed in a motorized stage microscope (EVOS Auto FL, Life Technologies, Carlsbad, CA) equipped with an incubator and maintained at 5% CO2 and 37 °C. Each well was imaged every 30 min for a total duration of up to 6 days.
Digital image analysis algorithm. A digital image analysis algorithm60 was implemented to determine changes in viability of each well longitudinally across the 96 or 144 h intervals. This algorithm computes differences in sequential images and identifies live cells with continuous membrane deformations resulting from their interaction with the surrounding extracellular matrix. These interactions cease upon cell death. By applying this operation to all 288 images acquired for each well, we quantified non-destructively, and without the need to separate the stroma and myeloma, the effect of drugs as a function of concentration and exposure time.
Estimating ex vivo drug sensitivity metrics: AUC, LD50, and LD20. Digital image analysis computes percent viability of MM cells for each time point and experimental condition (drug and concentration). For each patient-drug, we have a dose-time-response surface, which is abstracted into two metrics AUC (area under the curve) and LD50/LD20 (lethal dose for 50%/20% cell kill). AUC is an area/integral measure of ex vivo response to therapy computed by taking an average of all ex vivo responses across all time and concentration, while LD50 is a point measure of response characterized by the concentration that leads to 50% (20% in case of LD20) cell kill at a given point in time (96 h unless otherwise stated).
Sample preparation: Fresh BM aspirate cells were enriched for CD138 expression using Miltenyi (Bergisch Gladbach, Germany) 130-051-301 antibody-conjugated magnetic beads. 1.0 x 106 viably frozen CD138+ cells were shipped for molecular analysis in the context of the ORIEN Avatar program.
Nucleic Acid Extraction: For frozen tissue DNA extraction, Qiagen QIASymphony DNA purification was performed, generating 213 bp average insert size. For frozen tissue RNA extraction, Qiagen RNAeasy plus mini kit was performed, generating 216 bp average insert size.
DNA Sequencing: Preparation of Whole Exome Sequencing (WES) libraries involved hybrid capture using an enhanced IDT WES and Nimblegen SeqCap EZ kits (38.7 Mb) with additional custom designed probes for double coverage of 440 cancer genes. Library hybridization was performed at either single or 8-plex and sequenced on an Illumina NovaSeq 6000 instrument, generating 100 bp paired reads. WES was performed on tumor/normal matched samples with the normal covered at 100X and the tumor covered at 300X (additional 440 cancer genes covered at 600X) depth. Both tumor/normal concordance and gender identity QC checks are performed. Minimum threshold for hybrid selection is >80% of bases with >20X fold coverage; ORIEN Avatar WES libraries typically meet or exceed 90% of bases with >50X fold coverage for tumor and 90% of bases with >30X fold coverage for normal samples.
RNA Sequencing (Moffitt cohort & MCC17814) was performed using the Illumina TruSeq RNA Exome with single library hybridization, cDNA synthesis, library preparation, sequencing (at either 100 or 150 bp paired reads) to a coverage of 100M total reads / 50M paired reads.
RNA Sequencing (BOSTON): was performed on CD138+ bone marrow cells. RNA was extracted using the Qiagen Allprep RNA mini kit, and library preparation was performed with either the TruSeq Stranded mRNA (non-FFPE compatible), or the Smart-Seq V4 Ultra Low Input Nextera XT kit. Total RNA sequencing was performed with 100 bp reads using an Illumina HiSeq 2500 instrument to a coverage of at least 24M total reads / 12M paired reads.
Mutation calling: Moffitt cohort. Individual VCF files (one per sample) were converted to tab separated format using the software vcf2tsv (https://github.com/sigven/vcf2tsv version=0.3.4). Only gene mutation records with column value “PASS” for field “FILTER” and “exonic status” as “exonic” and type “protein coding” were considered. All files were merged and formatted according to minimum requirements and processed using the R package maftools (bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html version 2.4.12). Mutational summaries for non-synonymous mutations were created using the maftools functions oncoplot and plotmafSummary, while “lollipop” visualizations of individual gene mutation sites were generated with maftools function lollipopPlot. Over- and under-mutated genes were assessed by ratio of number of mutated samples by the length of the protein, in terms of number of amino acids.
Association of ex vivo drug sensitivity with disease state, cytogenetic abnormalities, and driver mutations. For this purpose, we split the ex vivo patient cohort into several groups, where each ex vivo drug response measure (AUC, LD50, and LD20) is compared between those patients that belong to that subgroup versus those who don’t via an unpaired t-test that yields a p-value for the comparison and a difference in median for AUCs and a log2 fold-change in median LD50/LD20. The p-values are adjusted for multi-test correction using Benjamin-Hochberg method and controlled for a false discovery rate less than 5% (FDR ≤ 0.05).
RNA sequencing (RNA-seq) analysis: Moffitt cohort & MCC17814. RNAseq Tumor Pipeline Analysis is processed according to the workflow outlined below using GRCh38/hg38 human genome reference sequencing and GenCode build version 32.
Adapter trimming: Adapter sequences are trimmed from the raw tumor sequencing FASTQ file. Adapter-trimming via k-mer matching is performed along with quality-trimming and filtering, contaminant-filtering, sequence masking, GC-filtering, length filtering and entropy-filtering. The trimmed FASTQ file is used as input to the read alignment process.
Read Alignment: The tumor adapter-trimmed FASTQ file is aligned to the human genome reference (GRCh38/hg38) and the Gencode genome annotation v32 using the STAR aligner. The STAR aligner generates multiple output files used for Gene Fusion Prediction and Gene Expression Analysis.
RNA expression: RNA expression values are calculated and reported using estimated mapped reads, Fragments Per Kilobase of transcript per Million mapped reads (FPKM), and Transcripts Per Million mapped reads (TPM) at both transcript level and gene level based on transcriptome alignment generated by STAR. Gene expression data was obtained from DNAnexus files containing FPKM and TPM values for 59,368 records. Among these, 19,933 were protein coding genes, which were further analyzed; the remainder were discarded. For each gene/sample, we have calculated log2(FPKM+10-3) and removed any genes whose values for quartile 1 and quartile 3 were the same (i.e., any gene must be expressed in at least 25% of samples to be considered in this analysis). The remaining 16,738 genes were Z-normalized across all samples using MATLAB’s function normalize.
RNA sequencing (RNA-seq) analysis: BOSTON
RNAseq Analysis for the BOSTON study was performed as follows. Trimming: Adapter and quality trimming were performed on the FASTQ files using trimmomatic61. Read Alignment: The adapter-trimmed FASTQ files were aligned to the human genome reference (GRCh38/hg38) using the STAR aligner. RNA Expression: FeatureCounts was used to obtain gene-level counts. After removing immunoglobulin genes, ribosomal transcripts, genes with zero counts, and genes with zero variance across all samples, the log2 counts per million mapped reads (CPM) values were normalized with voom62. Covariates were identified using variancePartition63, then corrected for using from sva64, including ComBat to remove batch effects and library preparation bias effects.
Identifying transcriptomic footprints associated with ex vivo drug resistance or sensitivity. To unravel the transcriptomic topology of a complex disease like MM, we employed RNA-seq data from 844 patients to identify modules of co-expressing genes using a robust dimensionality reduction technique and an efficient clustering method. Z-normalized expression of 16,738 genes across 844 MM patients is used to identify groups of co-expressing genes that are likely to play disease-specific functional roles. Although it is typical to consider genes as variables (dimensions) and patients as observations (typically used with single-cell sequencing data to identify clones, or cell types), we instead perceive patients as variables uniquely contributing to a high-dimensional MM heterogeneity space and genes as observations that govern MM transcriptomic topology. This leaves us with 16,738 genes spread across a massively high-dimensional (844) MM patient space. We project this high-dimensional data onto a two-dimensional space65 using t-Student Neighbor Embedding (tSNE)41, a well-known dimensionality reduction technique that specializes in extracting features (co-expression of genes) that lie on various low-dimensional embedded manifolds41, thereby serving as an excellent visualization tool depicting a disease-specific two-dimensional transcriptomic map. Locations of genes on this 2D map are used to identify functional modules by employing an efficient clustering algorithm called fuzzy C-means42, which results in 500 distinct gene clusters (gene sets) of varying sizes as shown in Extended Data Fig. 2.
We use the ex vivo drug sensitivity data for each drug from EMMA to identify clusters that are enriched for resistance and sensitivity using gene set enrichment analysis (GSEA)66. GSEA estimates an enrichment score for each cluster (gene set) using a running-sum statistic along the ranked-list of all genes (16,738) based on the correlation between their expression and the continuous phenotypic variable (AUC). Whenever GSEA encounters a gene that belongs to the cluster, it increases the running-sum statistic and decreases it if it doesn’t encounter a gene from that cluster. The maximum value of this running-sum statistic is the enrichment score for that cluster associated with positive correlation to the continuous phenotypic variable. GSEA estimates the statistical significance of such an enrichment by randomly scrambling the phenotypic variable several times, and for each case it generates a ranked gene list and the corresponding enrichment score for the cluster of interest. All these enrichment scores form a null distribution, which is compared to the enrichment score for the cluster using the actual input data to estimate the nominal p-value of enrichment. This approach is repeated for all 500 clusters and their nominal p-values are corrected for multiple hypothesis testing. The clusters that are enriched for resistance and sensitivity are identified by a family-wise error rate that is less than 5%. These enriched clusters collectively form the transcriptomic footprint of the drug in MM.
Transcription factors and epigenetic alterations enriched for a transcriptomic footprint. The transcriptomic footprint identified from paired RNA-seq and ex vivo drug response data for a drug using MM-specific transcriptomic landscape is a collection of co-expressing gene clusters in MM that are either enriched for sensitivity or resistance to a given drug. From this transcriptomic footprint we derive two gene sets: genes overexpressed in resistant patients and genes overexpressed in sensitive patients. We subject each of these two gene sets to enrichment analysis using enrichR44-46 by estimating the adjusted p-value for a one-sided Fisher-exact test quantifying the significance of overlap between genes implicated in ex vivo resistance/sensitivity and each of ChEA 201667 (for transcription factors that bind to genes implicated in ex vivo resistance/sensitivity) and epigenomics roadmap68 (epigenetic alterations that impact the expression of genes implicated in ex vivo resistance/sensitivity) supervised gene set databases. If the adjusted p-value is less than 10-10 (Moffitt cohort) or 0.05 (BOSTON and MCC17814) for a given test, that transcription factor or epigenetic alteration is considered statistically significantly enriched for either resistance or sensitivity to a given drug. This significance is denoted by the size of the bubble in the bubble plots featured in Figs 2d,e and 3c,d, while the intensity of the color signifies the combined enrichR score. Bubbles shown in blue represent enrichment analyses done on genes implicated in sensitivity to a given drug, while red bubbles signify enrichment associated with genes implicated in resistance phenotype.
Biomarker score estimation from transcriptomic footprint: Regression tree modeling. We relied on a regression tree model to estimate a molecular biomarker score from the median expression of the genes in each enriched cluster of the transcriptomic footprint as inputs to the model. The model was trained using the extensive paired ex vivo AUC response as output and RNA-seq data as input from the Moffitt cohort, where the predicted ex vivo AUC response from gene expression data would serve as a surrogate for clinical response in MM. The regression tree modeling was implemented in Matlab computational environment using the regression tree class available in the Statistics and Machine Learning Toolbox.